Introduction

Neuromorphic electronics aim at the realization of intelligent systems that emulate the immense capability of the nervous system to efficiently cope with a wide diversity of environmental and biological signals1,2,3. By emulating the brain’s basic building blocks such as synapses and neurons, neuromorphic electronics can perform bio-inspired processing and computation, opening opportunities in a broad range of application fields, including edge computing4, wearables5, point-of-care diagnostics6,7, bioelectronics8,9, (bio)robotics and environmental intelligence10. A wide range of materials has been used over the past decade for the realization of neuromorphic electronics including metal oxides11,12,13,14,15,16, ferroelectrics17,18, ferromagnetics19,20, phase change21,22,23, and 2D materials24,25.

Neuromorphic electronics made of organic materials are of particular interest because of their close resemblance with biology26. Neuromorphic electronics with soft electrochemical matter such as organic mixed ionic electronic conductors (OMIECs), can emulate realistically biological phenomena because of their responsiveness to the biological carriers of information (alkaline ions, neurotransmitters, neuromodulators, etc.)26,27,28,29,30. Moreover, their operation in wet biological environments opens entirely new directions for neuromorphic biosensors and bioelectronics. Indeed, over the past few years, the basic building blocks of the nervous system have been realized with OMIECs and their corresponding devices (i.e., organic electrochemical transistors or OECTs) such as low-voltage artificial synapses and synaptic networks31,32.

Organic electrochemical artificial neurons (OANs) are the latest entry of building blocks, with a few different approaches for circuit realization. OANs possess the remarkable capability to realistically mimic biological phenomena by responding to key biological information carriers, including alkaline ions, noise in the electrolyte, and biological conditions. An organic artificial neuron with a cascade-like topology made of OECT inverters has shown basic (regular) firing behavior and firing frequency that is responsive to the concentration of ionic species (Na+, K+) of the host liquid electrolyte33. An organic artificial neuron consisting of a non-linear building block that displays S-shape negative differential resistance (S-NDR) has also been recently demonstrated34. Due to the realization of the non-linear circuit theory with OECTs and the sharp threshold for oscillations, this artificial neuron displays biorealistic firing properties and neuronal excitability that can be found in the biological domain such as input voltage-induced regular and irregular firing, ion and neurotransmitter-induced excitability and ion-specific oscillations. Biohybrid devices comprising artificial neurons and biological membranes have also shown to operate synergistically, with membrane impedance state modulating the firing properties of the biohybrid in situ. More recently, a circuit leveraging the non-linear properties of antiambipolar OMIECs, which exhibit negative differential transconductance, has been realized35. These neurons show biorealistic properties such as various firing modes and responsivity to biologically relevant ions and neurotransmitters. With this neuron, ex-situ electrical stimulation has been shown in a living biological model. Therefore, the class of OANs perfectly complements the broad range of features already demonstrated by solid-state spiking circuits (Supplementary Table 1), offering opportunities for both hybrid interfacing between these technologies and new developments in neuromorphic bioelectronics.

Despite the promising recent realizations of organic artificial neurons, all approaches still remain in the qualitative demonstration domain and a rigorous investigation of circuit operation is still missing. Indeed, quantitative models exist only for inorganic, solid-state artificial neurons without the inclusion of physical soft-matter parameters and the biological wetware (i.e., aqueous electrolytes, alkaline ions, biomembranes)36,37. This gap in knowledge significantly impedes the simulation of larger-scale functional circuits, and therefore the design and development of integrated organic neuromorphic electronics, biohybrids, OAN-based neural networks, and intelligent bioelectronics.

In this work, we unravel the operation of organic artificial neurons that display non-linear phenomena such as S-shape negative differential resistance (S-NDR). By combining experiments, numerical simulations of non-linear iontronic circuits, and newly developed analytical expressions, we investigate, reproduce, rationalize, and design the wide biorealistic repertoire of organic electrochemical artificial neurons including their firing properties, neuronal excitability, wetware operation, and biohybrid formation. The OAN operation is efficiently rationalized to include how neuronal dynamics are probed by biochemical stimuli in the electrolyte medium. The OAN behavior is also extended on the biohybrid formation, with a solid rationale of the in situ interaction of OANs with biomembranes. Non-linear simulations of OANs are rooted in a physics-based framework, considering ion type, ion concentration, organic mixed ionic–electronic parameters, and biomembrane properties. The derived analytical expressions establish a direct link between OAN spiking features and its physical parameters and therefore provide a mapping between neuronal behavior and materials/device parameters. The proposed approach open opportunities for the design and engineering of advanced biorealistic OAN systems, establishing essential knowledge and tools for the development of neuromorphic bioelectronics, in-liquid neural networks, biohybrids, and biorobotics.

Results

The OAN is obtained by connecting an organic electrochemical non-linear device (OEND) with a biasing network comprising a resistor RL, a capacitor CL, and a DC voltage generator VIN (Fig. 1a). The OEND comprises two organic electrochemical transistors (OECTs), named T1 and T2, and two resistors, named R1 and R2. We used the mixed ionic–electronic conductors poly(3,4-ethylenedioxythiophene) doped with poly(styrene sulfonate) (PEDOT:PSS) and poly(2-(3,3-bis(2-(2-(2-methoxyethoxy)ethoxy)ethoxy)-[2,2-bithiophen]−5-yl) thieno [3,2-b] thiophene) (p(g2T-TT)) as channel materials for T1 and T2, respectively. T1 is a normally on, viz. depletion-mode, OECT while T2 is a normally off, viz. accumulation mode, OECT. Typical OECT transfer characteristics (\({I}_{{{{{{\rm{D}}}}}}}-{V}_{{{{{{\rm{G}}}}}}}\)) measured at various drain voltages (\({V}_{{{{{{\rm{D}}}}}}}\)) are displayed in Supplementary Fig. 1 of the Supplementary Information.

Fig. 1: Organic artificial neuron (OAN).
figure 1

a Circuit diagram of the OAN highlighting the organic electrochemical non-linear device (OEND). VIN is the load voltage generator, RL represents the load resistor, CL represents the load capacitance, T1 is the normally on OECT, T2 is the normally off OECT, and R1 and R2 are the two resistors of the OEND. b Measured (symbols) and simulated (full line) electrical characteristic of the OEND accessed in voltage mode: \({V}_{{{{{{\rm{OEND}}}}}}}\) is applied and \({I}_{{{{{{\rm{OEND}}}}}}}\) is measured. Arrows show the direction of \({V}_{{{{{{\rm{OEND}}}}}}}\) sweep. c Measured (symbols) and simulated (full line) electrical characteristic of the OEND accessed in current mode: \({I}_{{{{{{\rm{OEND}}}}}}}\) is applied and \({V}_{{{{{{\rm{OEND}}}}}}}\) is measured. Arrows show the direction of \({I}_{{{{{{\rm{OEND}}}}}}}\) sweep. d Measured (symbols) and simulated (full line) spiking current (\({I}_{{{{{{\rm{spike}}}}}}}\)) of the OAN. e Measured (symbols) and simulated (full line) spiking voltage (\({V}_{{{{{{\rm{spike}}}}}}}\)). T1 is a p-type depletion-mode OECTs based on PEDOT:PSS, and T2 is a p-type accumulation mode OECT based on p(g2T-TT).

To investigate the OAN, we model the OECT electrical characteristics (see Supplementary Note 1) accounting for material and device parameters such as the volumetric capacitance \( C_{{{\rm{V}}}}\)38,39, ion-concentration-dependent threshold voltage \( V_{{{\rm{TH}}}}\)40, energy disorder \( \gamma \)41,42, and channel-length modulation43,44. Then, we implemented the OECT model in a circuit simulator. As displayed in Fig. 1b, c, we accurately reproduced the measured electrical characteristics of the OEND. Voltage vs. current operation mode of the OEND is a typical measurement conducted for characterizing non-linear devices, such as redox-diffusive memristors, Mott memristors, single-transistor latch, and Gaussian heterojunction transistors45. Serving as the non-linear core component of the OAN, the OEND exhibits a behavior similar to other spiking neuron devices based on non-linear elements, displaying an abrupt increase or decrease in either the current-voltage or voltage-current relationship. The current-voltage characteristic of the OEND is crucial for designing and shaping the spiking behavior of the OAN. Specifically, the OEND can be assessed either as \({I}_{{{{{{\rm{OEND}}}}}}}({V}_{{{{{{\rm{OEND}}}}}}})\) where \({V}_{{{{{{\rm{OEND}}}}}}}\) is the independent (input) variable (Fig. 1b), or as \({V}_{{{{{{\rm{OEND}}}}}}}({I}_{{{{{{\rm{OEND}}}}}}})\) where \({I}_{{{{{{\rm{OEND}}}}}}}\) is the independent variable (Fig. 1c). When the OEND is operated in voltage mode, \({V}_{{{{{{\rm{OEND}}}}}}}\) is the independent input variable: \({V}_{{{{{{\rm{OEND}}}}}}}\) is swept forward and backward and the current \({I}_{{{{{{\rm{OEND}}}}}}}\) flowing through the OEND is recorded. As displayed in Fig. 1b, a hysteretic characteristic is obtained under steady-state operation, which is inherently attributed to the non-linear circuit configuration. Conversely, when the OEND is operated in current mode, \({I}_{{{{{{\rm{OEND}}}}}}}\) is the independent input variable: \({I}_{{{{{{\rm{OEND}}}}}}}\) is swept forward and backward and the voltage \({V}_{{{{{{\rm{OEND}}}}}}}\) is recorded. As displayed in Fig. 1c the \({V}_{{{{{{\rm{OEND}}}}}}}({I}_{{{{{{\rm{OEND}}}}}}})\) characteristic is non-hysteretic and shows S-shape negative differential resistance (S-NDR). A comprehensive discussion on OEND operation is provided later in the manuscript (“Organic electrochemical non-linear device operation” section).

Finally, we connected the OEND to the biasing network and we performed non-linear transient simulations of the OAN. As displayed in Fig. 1d, e the simulations predict both the experimental current and voltage output oscillations. The material and device parameters are listed in the “Methods” section unless otherwise stated. As we recently reported34, the artificial action potentials of the OAN and its firing properties, including excitability, spike latency, oscillation frequency, and amplitude, are sensitive to the local electro-bio-chemical signals in the liquid environment, viz. the extracellular space.

Organic electrochemical non-linear device operation

As a first step, we focus on the OEND, which is the core component of the OAN. When operating the OEND in voltage mode, \({V}_{{{{{{\rm{OEND}}}}}}}\) serves as the independent input variable. \({V}_{{{{{{\rm{OEND}}}}}}}\) is swept forward and backward, while recording the current \({I}_{{{{{{\rm{OEND}}}}}}}\) flowing through the OEND. As illustrated in Fig. 2a, under steady-state conditions, a hysteretic \({I}_{{{{{{\rm{OEND}}}}}}}({V}_{{{{{{\rm{OEND}}}}}}})\) characteristic is observed. This hysteresis is inherently due to the non-linear switching of the OETCs. The OEND internal voltages controlling the operation of the OECTs T1 and T2 are marked in Fig. 2b and are named \({V}_{{{{{{\rm{GS}}}}}}1}\), \({V}_{{{{{{\rm{SG}}}}}}2}\), and \({V}_{{{{{{\rm{GD}}}}}}2}\). The OEND internal voltages calculated by means of numerical simulations are displayed in Fig. 2c, d and Supplementary Fig. 2. For the sake of clarity, the operating regions of the OECTs are also highlighted. During the forward sweep of \({V}_{{{{{{\rm{OEND}}}}}}}\) from 0 V to positive voltages lower than \({V}_{{{{{{\rm{ON}}}}}}}\), the current \({I}_{{{{{{\rm{OEND}}}}}}}\) flows through the branch R1-T1-R2 (region 1, Fig. 2a). This occurs because T1 is a depletion-mode p-type OECT (\({V}_{{{{{{\rm{TH}}}}}}1} \, > \, 0 \,{{{\rm{V}}}}\)) and T2 is an accumulation mode p-type OECT (\({V}_{{{{{{\rm{TH}}}}}}2} \, < \, 0 \,{{{\rm{V}}}}\)). Specifically, at small \({V}_{{{{{{\rm{OEND}}}}}}}\), T1 operates in the linear region, resulting in a small channel resistance \({R}_{{{{{{\rm{T}}}}}}1}\) and consequently a small source-drain voltage \({V}_{{{{{{\rm{SD}}}}}}1}={R}_{{{{{{\rm{T}}}}}}1}{I}_{{{{{{\rm{OEND}}}}}}}\) is obtained. The circuit topology dictates that \({V}_{{{{{{\rm{SD}}}}}}1}={V}_{{{{{{\rm{SG}}}}}}2}\), and hence \({V}_{{{{{{\rm{SG}}}}}}2} < |{V}_{{{{{{\rm{TH}}}}}}2}|\) (Fig. 2d). During the forward sweep of \({V}_{{{{{{\rm{OEND}}}}}}}\), as long as \({V}_{{{{{{\rm{OEND}}}}}}} < {V}_{{{{{{\rm{ON}}}}}}}\) (region 1, Fig. 2a), T1 remains ON while T2 remains OFF, causing \({I}_{{{{{{\rm{OEND}}}}}}}\) to flow through the branch R1-T1-R2, resulting in a slope of the \({I}_{{{{{{\rm{OEND}}}}}}}({V}_{{{{{{\rm{OEND}}}}}}})\) characteristic of \({{{{{\rm{d}}}}}}{I}_{{{{{{\rm{OEND}}}}}}}/{{{{{\rm{d}}}}}}{V}_{{{{{{\rm{OEND}}}}}}}=1/({R}_{1}+{R}_{{{{{{\rm{T}}}}}}1}+{R}_{2})\).

Fig. 2: OEND operation.
figure 2

a OEND electrical characteristic accessed in voltage mode. A voltage ramp (\({V}_{{{{{{\rm{OEND}}}}}}}\)) is applied forward (red line) and backward (blue line), and the current \({I}_{{{{{{\rm{OEND}}}}}}}\) is measured. The four operating regions and the relevant circuit components in each region of operation are highlighted. b OEND circuit showing the internal voltages. c Gate-source voltage of transistor T1 (\({V}_{{{{{{\rm{GS}}}}}}1}\)) as a function of \({V}_{{{{{{\rm{OEND}}}}}}}\). T1 threshold voltage (\({V}_{{{{{{\rm{TH}}}}}}1}\)) is displayed. If \({V}_{{{{{{\rm{GS}}}}}}1}\ge {V}_{{{{{{\rm{TH}}}}}}1}\), T1 is turned OFF while if \({V}_{{{{{{\rm{GS}}}}}}1} < {V}_{{{{{{\rm{TH}}}}}}1}\), T1 is operated in a linear regime. Forward sweep (red line) and backward sweep (dashed blue line). The numbers refer to the four operating regions of the OEND, as highlighted in (a). d Source-gate voltage of transistor T2 (\({V}_{{{{{{\rm{SG}}}}}}2}\)) as a function of \({V}_{{{{{{\rm{OEND}}}}}}}\). If \({V}_{{{{{{\rm{SG}}}}}}2}\le |{V}_{{{{{{\rm{TH}}}}}}2}|\), T2 is in the OFF state and, if \({V}_{{{{{{\rm{SG}}}}}}2} > |{V}_{{{{{{\rm{TH}}}}}}2}|\) results that T2 is in saturation regime. e OEND characteristic accessed in current mode calculated with numerical simulations (line) and measured (symbols). A current ramp (\({I}_{{{{{{\rm{OEND}}}}}}}\)) is applied forward and backward, and the voltage \({V}_{{{{{{\rm{OEND}}}}}}}\) is measured. Forward and backward voltages are overlapped. The points \({{{{{{\bf{U}}}}}}}_{{{{{{\rm{ON}}}}}}}=({V}_{{{{{{\rm{ON}}}}}}},{I}_{{{{{{\rm{ON}}}}}}})\) and \({{{{{{\bf{U}}}}}}}_{{{{{{\rm{OFF}}}}}}}=({V}_{{{{{{\rm{OFF}}}}}}},{I}_{{{{{{\rm{OFF}}}}}}})\) define the negative resistance region (NRD).

As \({V}_{{{{{{\rm{OEND}}}}}}}\) increases, the source-drain voltage T1 (\({V}_{{{{{{\rm{SD}}}}}}1}\)) also increases. Once \({V}_{{{{{{\rm{SD}}}}}}1} > {|}{V}_{{{{{{\rm{TH}}}}}}2}|\), T2 turns ON and operates in the saturation region (Fig. 2d). By design, when T2 is ON, its channel resistance (\({R}_{{{{{{\rm{T}}}}}}2}\)) becomes much smaller than \({R}_{2}\), causing \({I}_{{{{{{\rm{OEND}}}}}}}\) to predominantly flow through the branch R1-T2, leading to a sharp increase in current (region 2 in Fig. 2a). This non-linear current enhancement results in a substantial voltage drop across \({R}_{1}\). Being \({V}_{{{{{{\rm{R}}}}}}1}={V}_{{{{{{\rm{GS}}}}}}1}\) (Fig. 2c), the overdrive voltage of T1 decreases and eventually \({V}_{{{{{{\rm{GS}}}}}}1}\, > \,{V}_{{{{{{\rm{TH}}}}}}1}\) causes T1 to turn OFF (Fig. 2c). The voltage required to turn ON T2 is referred to as \({V}_{{{{{{\rm{ON}}}}}}}\). With further increase in \({V}_{{{{{{\rm{OEND}}}}}}}\) beyond \({V}_{{{{{{\rm{ON}}}}}}}\) (region 3 in Fig. 2a), \({I}_{{{{{{\rm{OEND}}}}}}}\) flows through the branch R1-T2, with a linear increase characterized by a slope \({{{{{\rm{d}}}}}}{I}_{{{{{{\rm{OEND}}}}}}}/{{{{{\rm{d}}}}}}{V}_{{{{{{\rm{OEND}}}}}}}=1/({R}_{1}+{R}_{{{{{{\rm{T}}}}}}2})\).

When \({V}_{{{{{{\rm{OEND}}}}}}}\) is swept back to lower voltages, \({I}_{{{{{{\rm{OEND}}}}}}}\) linearly decreases, resulting in a decrease in the voltage \({V}_{{{{{{\rm{GS}}}}}}1}={R}_{1}{I}_{{{{{{\rm{OEND}}}}}}}\). When \({V}_{{{{{{\rm{GS}}}}}}1}\, < \,{V}_{{{{{{\rm{TH}}}}}}1}\), T1 switches ON in the linear region. Under this condition, \({R}_{{{{{{\rm{T}}}}}}1}\) is small, and \({V}_{{{{{{\rm{SD}}}}}}1}={R}_{{{{{{\rm{T}}}}}}1}{I}_{{{{{{\rm{OEND}}}}}}}\) becomes small as well. Given that \({V}_{{{{{{\rm{SD}}}}}}1}={V}_{{{{{{\rm{SG}}}}}}2}\), T2 turns OFF when \({V}_{{{{{{\rm{SG}}}}}}2} < |{V}_{{{{{{\rm{TH}}}}}}2}|\), causing the OEND current to sharply decrease (region 4 in Fig. 2a). The voltage required to turn OFF T2 is referred to as \({V}_{{{{{{\rm{OFF}}}}}}}\). Thus, the OEND exhibits a hysteretic characteristic due to the non-linear current switching from branch R1-T1-R2 to branch R1-T2 in the forward voltage sweep and from branch R1-T2 to R1-T1-R2 in the backward voltage sweep.

It is important to note that the hysteresis is observed under steady-state conditions (DC operation) as it is inherently related to the OEND circuit configuration, ensuring that \({V}_{{{{{{\rm{ON}}}}}}}\) and \({V}_{{{{{{\rm{OFF}}}}}}}\) occur at different voltages. The analytical expressions of \({V}_{{{{{{\rm{ON}}}}}}}\) and \({V}_{{{{{{\rm{OFF}}}}}}}\) as a function of material, geometrical, and device parameters are derived in Supplementary Note 2 and are as follows:

$${V}_{{{{{{\rm{ON}}}}}}} \cong ({{R}_{1}+R}_{2}){g}_{{{{{{\rm{m}}}}}}1}\frac{{V}_{{{{{{\rm{TH}}}}}}2}(2{V}_{{{{{{\rm{TH}}}}}}1}+{V}_{{{{{{\rm{TH}}}}}}2})}{2\,{g}_{{{{{{\rm{m}}}}}}1}\,{R}_{1}{V}_{{{{{{\rm{TH}}}}}}2}\,-\,1}-{V}_{{{{{{\rm{TH}}}}}}2}$$
(1)
$${V}_{{{{{{\rm{OFF}}}}}}} \cong {V}_{{{{{{\rm{TH}}}}}}1}-{V}_{{{{{{\rm{TH}}}}}}2}+\root {{{{\mathrm{\gamma 2}}}}} \of {{\frac{{V}_{{{{{{\rm{TH}}}}}}1}}{{g}_{{{{{{\rm{m}}}}}}2}\,{R}_{1}}}}$$
(2)

where \({g}_{{{{{{\rm{m}}}}}}1}=\left({W}_{1}/{L}_{1}\right){t}_{1}{\mu }_{1}{C}_{{{{{{\rm{V}}}}}}1}/2\) is the transconductance of the OECT T1 normalized to \({V}_{{{{{{\rm{SD}}}}}}1}\), \({g}_{{{{{{\rm{m}}}}}}2}\) the normalized transconductance of the OECT T2, and \({W}_{1}\), \({L}_{1}\), \({t}_{1}\), \({\mu }_{1}\), \({C}_{{{{{{\rm{V}}}}}}1}\) is the width, length, thickness, mobility, and volumetric capacitance of T1, respectively. As depicted in Supplementary Fig. 3 the amplitude of the OAN voltage oscillations can be calculated as \({A}_{{{{{{\rm{Vspike}}}}}}}\,=\,{V}_{{{{{{\rm{ON}}}}}}}\,-\,{V}_{{{{{{\rm{OFF}}}}}}}\). This information can be obtained from the OEND characteristic assessed in voltage mode \({I}_{{{{{{\rm{OEND}}}}}}}\)(\({V}_{{{{{{\rm{OEND}}}}}}}\)). The operation of the OAN necessitates that \({V}_{{{{{{\rm{ON}}}}}}} > 0\,{{{{{\rm{V}}}}}}\), \({V}_{{{{{{\rm{OFF}}}}}}}\, > \,0\,{{{{{\rm{V}}}}}}\), and \({V}_{{{{{{\rm{ON}}}}}}}\, > \,{V}_{{{{{{\rm{OFF}}}}}}}\). Equations (1) and (2) demonstrate that these conditions are satisfied when \({V}_{{{{{{\rm{TH}}}}}}1}\, > \,0\,{{{{{\rm{V}}}}}}\), \({V}_{{{{{{\rm{TH}}}}}}2}\, < \,0\,{{{{{\rm{V}}}}}}\), and \({V}_{{{{{{\rm{TH}}}}}}1} > |{V}_{{{{{{\rm{TH}}}}}}2}|\). Consequently, T1 has to be a depletion-mode OECT, while T2 has to be an accumulation mode OECT. To further validate the operational parameters of T1 and T2 required for the development of a functional OAN, Supplementary Fig. 4 illustrates that the absence of hysteresis in the \({I}_{{{{{{\rm{OEND}}}}}}}({V}_{{{{{{\rm{OEND}}}}}}})\) characteristic, mirrored by the negative resistance region in the \({V}_{{{{{{\rm{OEND}}}}}}}({I}_{{{{{{\rm{OEND}}}}}}})\) characteristic, occurs in the limiting case \({V}_{{{{{{\rm{TH}}}}}}1}=|{V}_{{{{{{\rm{TH}}}}}}2}|\).

When the OEND is operated in current mode, \({I}_{{{{{{\rm{OEND}}}}}}}\) acts as the independent input variable: \({I}_{{{{{{\rm{OEND}}}}}}}\) undergoes forward and backward sweeps while the voltage \({V}_{{{{{{\rm{OEND}}}}}}}\) is recorded. In Fig. 2e the measured data (symbols) and simulated results (line) for the \({V}_{{{{{{\rm{OEND}}}}}}}({I}_{{{{{{\rm{OEND}}}}}}})\) characteristic are presented. The simulations nicely predict the measurements and, unlike the \({I}_{{{{{{\rm{OEND}}}}}}}({V}_{{{{{{\rm{OEND}}}}}}})\) characteristic, the \({V}_{{{{{{\rm{OEND}}}}}}}({I}_{{{{{{\rm{OEND}}}}}}})\) characteristic is non-hysteretic and exhibits S-shape negative differential resistance. The S-NDR behavior is crucial for achieving the spiking behavior of the OAN and can be comprehended based on previous analyses. More in detail, at low input currents, \({I}_{{{{{{\rm{OEND}}}}}}}\) primarily flows through the left branch R1-T1-R2 causing \({V}_{{{{{{\rm{OEND}}}}}}}\) to increase almost linearly with the current (Region 1 in Fig. 2e). When \({V}_{{{{{{\rm{OEND}}}}}}}={V}_{{{{{{\rm{ON}}}}}}}\), T2 turns ON, \({I}_{{{{{{\rm{OEND}}}}}}}\) can also flow in the right branch R1-T2, leading to a lower voltage drop \({V}_{{{{{{\rm{OEND}}}}}}}\) across the OEND and resulting in a negative differential resistance (Region 2 in Fig. 2e). With further increases in \({I}_{{{{{{\rm{OEND}}}}}}}\), the overdrive voltage on T2 rises (\({V}_{{{{{{\rm{SG}}}}}}2}={V}_{{{{{{\rm{SD}}}}}}1}\)), causing the channel resistance \({R}_{{{{{{\rm{T}}}}}}2}\) to decrease and consequently reducing \({V}_{{{{{{\rm{OEND}}}}}}}\) until \({V}_{{{{{{\rm{R}}}}}}1}={I}_{{{{{{\rm{OEND}}}}}}}{R}_{1}\) becomes sufficiently large to deactivate T1. Subsequently, current flows solely through the R1-T2 branch, causing \({V}_{{{{{{\rm{OEND}}}}}}}\) to monotonically increase with \({I}_{{{{{{\rm{OEND}}}}}}}\) and restoring positive resistance (Region 3 in Fig. 2e). The non-linear partitioning of current between the two branches of the OEND explains the S-shaped non-linear characteristic and the currents at the switching points of the OEND are determined (refer to Supplementary Note 2):

$${I}_{{{{{{\rm{ON}}}}}}}\cong {g}_{{{{{{\rm{m}}}}}}1}\frac{{V}_{{{{{{\rm{TH}}}}}}2}(2{V}_{{{{{{\rm{TH}}}}}}1}+{V}_{{{{{{\rm{TH}}}}}}2})}{2\,{g}_{{{{{{\rm{m}}}}}}1}\,{R}_{1}{V}_{{{{{{\rm{TH}}}}}}2}-1}$$
(3)
$${I}_{{{{{{\rm{OFF}}}}}}}\cong \frac{{V}_{{{{{{\rm{TH}}}}}}1}}{{R}_{1}}\quad\;$$
(4)

The inflection points \({{{{{{\bf{U}}}}}}}_{{{{{{\rm{ON}}}}}}}=({V}_{{{{{{\rm{ON}}}}}}},\,{I}_{{{{{{\rm{ON}}}}}}})\) and \({{{{{{\bf{U}}}}}}}_{{{{{{\rm{OFF}}}}}}}=({V}_{{{{{{\rm{OFF}}}}}}},\,{I}_{{{{{{\rm{OFF}}}}}}})\) (Fig. 2e) given by Eqs. (14) define the beginning and the end of the negative resistance region of the OEND and, as will be shown in the next section, are associated with the excitation and inhibition behavior of the OAN.

Organic artificial neuron operation and spiking

The spiking activity of the OAN arises from the coupling between the OEND and the biasing network, as depicted in Fig. 3a. When \({V}_{{{{{{\rm{OEND}}}}}}} < {V}_{{{{{{\rm{IN}}}}}}}\), the bias current \({I}_{{{{{{\rm{B}}}}}}}=({V}_{{{{{{\rm{IN}}}}}}}-{V}_{{{{{{\rm{OEND}}}}}}})/{R}_{{{{{{\rm{L}}}}}}}\) charges the capacitor CL, causing the voltage across the capacitor (\({V}_{{{{{{\rm{spike}}}}}}}\)) to increase (Fig. 3b, green line). Since the OAN topology gives \({V}_{{{{{{\rm{spike}}}}}}}={V}_{{{{{{\rm{OEND}}}}}}}\), as \({V}_{{{{{{\rm{OEND}}}}}}}\) increases and the OEND operates in the negative resistance region, the current \({I}_{{{{{{\rm{OEND}}}}}}}\) decreases, allowing a larger fraction of current \({I}_{{{{{{\rm{B}}}}}}}\) (\({I}_{{{{{{\rm{B}}}}}}}={I}_{{{{{{\rm{C}}}}}}}-{I}_{{{{{{\rm{OEND}}}}}}}\)) to charge CL. This further increases \({V}_{{{{{{\rm{OEND}}}}}}}\) and when \({V}_{{{{{{\rm{OEND}}}}}}}\ge {V}_{{{{{{\rm{ON}}}}}}}\) the OEND current significantly increases, reaching the condition \({I}_{{{{{{\rm{OEND}}}}}}} > {I}_{{{{{{\rm{B}}}}}}}\) (Fig. 3c). Subsequently, CL is discharged and \({V}_{{{{{{\rm{OEND}}}}}}}\) nonlinearly decreases (Fig. 3b, pink line). When \({V}_{{{{{{\rm{OEND}}}}}}}\le {V}_{{{{{{\rm{OFF}}}}}}}\), \({I}_{{{{{{\rm{OEND}}}}}}}\) significantly decreases, and when \({I}_{{{{{{\rm{OEND}}}}}}} < {I}_{{{{{{\rm{B}}}}}}}\), CL is charged again. Therefore, the charging and discharging of the load capacitor depends on the non-linear characteristic of the OEND, input voltage \({V}_{{{{{{\rm{IN}}}}}}}\) (a DC voltage), and load resistor. The asymmetrical profile of the spikes is a result of the hysteretic profile of the OEND. The connection of the two blocks (i.e., OEND and biasing network) are analyzed in Fig. 3d accounting for the OEND (full blue line) and the biasing network (dashed black line) characteristics. Specifically, when the load-line characteristic crosses the OEND characteristic in the negative resistance region (e.g., point \({{{{{\bf{U}}}}}}\) in Fig. 3d), its response bifurcates, producing voltage and current oscillations.

Fig. 3: OAN operation and spiking.
figure 3

a OAN circuit highlighting the current partition when \({V}_{{{{{{\rm{spike}}}}}}}\) increases, and the capacitor is charged. VIN is the load voltage generator, RL represents the load resistor, CL represents the load capacitance. \({I}_{{{{{{\rm{C}}}}}}}\) represents the current of the capacitor CL, while \({I}_{{{{{{\rm{B}}}}}}}\) is the current flowing through RL. b Voltage oscillations (\({V}_{{{{{{\rm{spike}}}}}}}\)) as a function of time. c OAN circuit highlighting the current partition when \({V}_{{{{{{\rm{spike}}}}}}}\) decreases, and the capacitor is discharged. d OEND (blue full line) and load-line (dashed line) characteristics. When the load-line characteristic crosses the OEND characteristic in the negative differential resistance region (e.g., point \({{{{{\bf{U}}}}}}\)), its response bifurcates, producing voltage and current oscillations. The OEND negative differential resistance is highlighted by dot-dashed purple line. The spiking region is defined by the upper and lower points \({{{{{{\bf{U}}}}}}}_{{{{{{\rm{ON}}}}}}}=({V}_{{{{{{\rm{ON}}}}}}},{I}_{{{{{{\rm{ON}}}}}}})\) and \({{{{{{\bf{U}}}}}}}_{{{{{{\rm{OFF}}}}}}}=({V}_{{{{{{\rm{OFF}}}}}}},{I}_{{{{{{\rm{OFF}}}}}}})\), respectively. \({{{{{{\bf{S}}}}}}}_{1}\) and \({{{{{{\bf{S}}}}}}}_{2}\) are two points where the OEND characteristic shows a positive resistance. e Spiking frequency \({f}_{{{{{{\rm{spike}}}}}}}\) as a function of load capacitor (CL). Symbols are the measurements, full line is calculated with the numerical simulations and dashed line is calculated with the analytical model. f Spiking current as a function of time measured (symbols) and calculated with non-linear transient simulations (full line) in the case \({C}_{{{{{{\rm{L}}}}}}}\) = 1 µF, g \({C}_{{{{{{\rm{L}}}}}}}\) = 5 µF, and h \({C}_{{{{{{\rm{L}}}}}}}\) = 10 µF.

Figure 3d highlights that the spiking region (orange area) is defined by the points \({{{{{{\bf{U}}}}}}}_{{{{{{\rm{ON}}}}}}}=({V}_{{{{{{\rm{ON}}}}}}},{I}_{{{{{{\rm{ON}}}}}}})\) and \({{{{{{\bf{U}}}}}}}_{{{{{{\rm{OFF}}}}}}}=({V}_{{{{{{\rm{OFF}}}}}}},{I}_{{{{{{\rm{OFF}}}}}}})\), respectively, which can be analytically calculated with Eqs. (14). Therefore, OAN spiking activity is inhibited when a stimulus causes the crossing point to shift above \({{{{{{\bf{U}}}}}}}_{{{{{{\rm{OFF}}}}}}}\). Conversely, OAN spiking activity is triggered when a stimulus causes the crossing point to shift from below to above \({{{{{{\bf{U}}}}}}}_{{{{{{\rm{ON}}}}}}}\), and the distance between the initial position of the crossing point and \({{{{{{\bf{U}}}}}}}_{{{{{{\rm{ON}}}}}}}\) defines the OAN excitability (i.e., excitation threshold). Importantly, we note that a unique crossing point in the spiking region is required to avoid the simultaneous presence of non-spiking and spiking states. This condition is fulfilled when the OEND negative resistance \({R}_{{{{{{\rm{D}}}}}}} \, < \, {R}_{{{{{{\rm{L}}}}}}}\) (see the dot-dashed violet line and dashed black line in Fig. 3d). When the load line crosses the OEND characteristic in the positive resistance regions (e.g., point \({{{{{{\bf{S}}}}}}}_{1}\) or \({{{{{{\bf{S}}}}}}}_{2}\), Fig. 3d), stable biasing conditions (i.e., no spiking) are obtained. Indeed, in the bias point \({{{{{{\bf{S}}}}}}}_{1}\) the OEND circuit reduces to a resistor \({R}_{{{{{{\rm{OEND}}}}}}({{{{{{\bf{S}}}}}}}_{1})}={R}_{1}+{R}_{{{{{{\rm{T}}}}}}1}+{R}_{2}\) that, according with the OAN architecture, is connected in parallel to CL and in series to RL. The voltage across CL is equal to \({V}_{{{{{{\rm{OAN}}}}}}}={V}_{{{{{{\rm{IN}}}}}}}{R}_{{{{{{\rm{OEND}}}}}}({{{{{{\bf{S}}}}}}}_{1})}/({R}_{{{{{{\rm{OEND}}}}}}({{{{{{\bf{S}}}}}}}_{1})}+{R}_{{{{{{\rm{L}}}}}}})\) and therefore no oscillations are obtained. Analogously, in the bias point \({{{{{{\bf{S}}}}}}}_{2}\) the OEND behaves as a resistor \({R}_{{{{{{\rm{OEND}}}}}}\left({{{{{{\bf{S}}}}}}}_{2}\right)}={R}_{1}+{R}_{{{{{{\rm{T}}}}}}2}\), and once again a constant output voltage \({V}_{{{{{{\rm{OAN}}}}}}}={V}_{{{{{{\rm{IN}}}}}}}{R}_{{{{{{\rm{OEND}}}}}}\left({{{{{{\bf{S}}}}}}}_{2}\right)}/({R}_{{{{{{\rm{OEND}}}}}}\left({{{{{{\bf{S}}}}}}}_{2}\right)}+{R}_{{{{{{\rm{L}}}}}}})\) is obtained. Interestingly, the bias points \({{{{{{\bf{S}}}}}}}_{1}\) and \({{{{{{\bf{S}}}}}}}_{2}\) provide a bistability condition of the OAN.

Focusing on the OAN spiking condition, CL is charged and discharged by a transient current \({I}_{{{{{{\rm{C}}}}}}}=({V}_{{{{{{\rm{IN}}}}}}}-{V}_{{{{{{\rm{OEND}}}}}}})/{R}_{{{{{{\rm{L}}}}}}}-{I}_{{{{{{\rm{OEND}}}}}}}\), and integrating over the voltage-current loop the capacitor current equation \({I}_{{{{{{\rm{C}}}}}}}={C}_{{{{{{\rm{L}}}}}}}{{{{{\rm{d}}}}}}{V}_{{{{{{\rm{C}}}}}}}/{{{{{\rm{d}}}}}}t\), where \({V}_{{{{{{\rm{C}}}}}}}={V}_{{{{{{\rm{OEND}}}}}}}\) and \(t\) is the time, the spiking frequency of the OAN can be calculated as:

$${f}_{{{{{{\rm{spike}}}}}}}={\frac{1}{{R}_{L}{C}_{{{{{{\rm{L}}}}}}}}\left[{\int }_{{V}_{{{{{{\rm{OFF}}}}}}}}^{{V}_{{{{{{\rm{ON}}}}}}}}\frac{{{{{{\rm{d}}}}}}V}{{V}_{{{{{{\rm{IN}}}}}}}-V-{R}_{{{{{{\rm{L}}}}}}}{I}_{{{{{{\rm{spike}}}}}}1}\left(V\right)}+{\int }_{{V}_{{{{{{\rm{ON}}}}}}}}^{{V}_{{{{{{\rm{OFF}}}}}}}}\frac{{{{{{\rm{d}}}}}}V}{{V}_{{{{{{\rm{IN}}}}}}}-V-{R}_{{{{{{\rm{L}}}}}}}{I}_{{{{{{\rm{spike}}}}}}2}\left(V\right)}\right]}^{-1}$$
(5)

where \({I}_{{{{{{\rm{spike}}}}}}1}\, \approx \,\frac{V}{{R}_{1}+{R}_{2}+{R}_{{{{{{\rm{T}}}}}}1}}\) and \({I}_{{{{{{\rm{spike}}}}}}2}\, \approx \, \frac{V+{V}_{{{{{{\rm{TH}}}}}}2}}{{R}_{1}+{R}_{{{{{{\rm{T}}}}}}2}}\), is the OEND current during the charging and discharging of CL, respectively. For the sake of clarity, a detailed derivation of Eq. (5) and its analytical solution providing an explicit approximate expression of \({f}_{{{{{{\rm{spike}}}}}}}\) as a function of the various parameters is reported in Supplementary Note 3. Figure 3e shows the measured (symbols) \({f}_{{{{{{\rm{spike}}}}}}}\) as a function of \({C}_{{{{{{\rm{L}}}}}}}\). The comparison of the measurements with numerical simulations (Fig. 3e, full line) and analytical model given by Eq. (5) (Fig. 3e, dashed line) show that the OAN spiking frequency is accurately predicted by the non-linear transient simulations of the OAN and nicely estimated by the computationally inexpensive analytical model. Interestingly, Eq. (5) accounts for the influence of RL and CL (as expected), and quantifies the (less apparent) impact of the material and device parameters through the quantities \({V}_{{{{{{\rm{ON}}}}}}}\), \({V}_{{{{{{\rm{OFF}}}}}}}\), \({V}_{{{{{{\rm{IN}}}}}}}\), and \({I}_{{{{{{\rm{OEND}}}}}}}\). As a further confirmation, the measured (symbols) and simulated (lines) spiking current as a function of time is displayed in Fig. 3f–h by varying \({C}_{{{{{{\rm{L}}}}}}}\). OAN simulations accurately predict the amplitude, frequency, and shape of the measured spiking current. We note that the spiking frequency is a crucial aspect of neural activity. For instance, in sensory systems, the firing rate of neurons can represent various attributes of a stimulus, such as its intensity, duration, or location. Moreover, the spiking frequency of motor neurons is essential for controlling movement. The rate and pattern of spikes in these neurons determine the strength and timing of muscle contractions. Every class on neuron has its own characteristic range of firing frequency, while deviations might be an indication of pathological conditions. Therefore, it is essential to control the firing frequency of the OANs by design. The developed numerical framework serves as a valuable tool for the design of OANs and Eq. (5) and its analytical formulation given by Supplementary Eq. (26) offers a valuable model for the rational design of \({f}_{{{{{{\rm{spike}}}}}}}\) as a function of the OAN material, geometrical, physical and device parameters.

To gain more insight on the impact of the various material, geometrical, and device parameters on the OAN performance, we take advantage of numerical simulations. Specifically, we performed parametric analysis of the OAN performance as a function of threshold voltage (\({V}_{{{{{{\rm{TH}}}}}}1},\,{V}_{{{{{{\rm{TH}}}}}}2}\)), transconductance (\({g}_{{{{{{\rm{m}}}}}}1},\,{g}_{{{{{{\rm{m}}}}}}2}\)), volumetric capacitance (\({C}_{{{{{{\rm{V}}}}}}1},\,{C}_{{{{{{\rm{V}}}}}}2}\)), width (\({W}_{1},{W}_{2}\)), length (\({L}_{1},{L}_{2}\)), thickness (\({t}_{1},\,{t}_{2}\)), and resistance (\({R}_{1},{R}_{2}\)). To assess the impact of each parameter on the OAN performance, we modified one parameter at a time while keeping all other parameters constant. The parameters are varied accounting for the largest range that ensures OAN oscillation. The corresponding \({I}_{{{{{{\rm{OEND}}}}}}}\)-\({V}_{{{{{{\rm{OEND}}}}}}}\) characteristics accessed in the current mode are displayed in Supplementary Figs. 6, 7, and 8. Numerical simulations in Fig. 4a–h show that \({f}_{{{{{{\rm{spike}}}}}}}\) can be modulated in the range 10–100 Hz by tuning the parameters \({V}_{{{{{{\rm{TH}}}}}}1}\) (Fig. 4a), \({V}_{{{{{{\rm{TH}}}}}}2}\) (Fig. 4b), \({g}_{{{{{{\rm{m}}}}}}1}\) (Fig. 4c), \({C}_{{{{{{\rm{V}}}}}}1}\) (Fig. 4e), \({C}_{{{{{{\rm{V}}}}}}2}\) (Fig. 4f), and \({R}_{2}\) (Fig. 4h) while it is almost insensitive to \({g}_{{{{{{\rm{m}}}}}}2}\), (Fig. 4d) and \({R}_{1}\). (Fig. 4g). More in detail, \({f}_{{{{{{\rm{spike}}}}}}}\) increases by reducing \({V}_{{{{{{\rm{TH}}}}}}1}\), \({g}_{{{{{{\rm{m}}}}}}1}\), \({C}_{{{{{{\rm{V}}}}}}1}\), and \({R}_{2}\), while an opposite trend is obtained with \({V}_{{{{{{\rm{TH}}}}}}2}\). The relation between \({f}_{{{{{{\rm{spike}}}}}}}\) and the material parameters can be linked to the OEND characteristic. \({f}_{{{{{{\rm{spike}}}}}}}\) can be enhanced by minimizing the current flowing in the branch R1-T1-R2 and maximizing the current flowing in the branch R1-T2. Interestingly, \({C}_{{{{{{\rm{V}}}}}}1}\) (Fig. 4e) and \({C}_{{{{{{\rm{V}}}}}}2}\) (Fig. 4f) are very significant material parameters influencing \({f}_{{{{{{\rm{spike}}}}}}}\), with the spiking frequency being maximum when \({C}_{{{{{{\rm{V}}}}}}1}\) and \({C}_{{{{{{\rm{V}}}}}}2}\) are minimized. To further explore this aspect, Supplementary Fig. 9 analyzes the OAN spiking frequency while also considering the relationship between \({C}_{{{{{{\rm{V}}}}}}1}\), \({C}_{{{{{{\rm{V}}}}}}2}\), and load capacitor \({C}_{{{{{{\rm{L}}}}}}}\). We systematically varied both \({C}_{{{{{{\rm{V}}}}}}1}\) and \({C}_{{{{{{\rm{V}}}}}}2}\), and the minimum \({C}_{{{{{{\rm{L}}}}}}}\) required for OAN spiking is calculated. Supplementary Fig. 9a reveals that the minimum \({C}_{{{{{{\rm{L}}}}}}}\) amounts to 10−7 F and is achieved when \({C}_{{{{{{\rm{V}}}}}}2}\) falls in the range 50–100 F cm−3 and \({C}_{{{{{{\rm{V}}}}}}1}\) is approximately 50 F cm−3. Supplementary Fig. 9b highlights that the maximum spiking frequency, \({f}_{{{{{{\rm{spike}}}}}}}=150 \, {{{\rm{Hz}}}}\), is attained when both \({C}_{{{{{{\rm{V}}}}}}2}\) and \({C}_{{{{{{\rm{L}}}}}}}\) are minimized. We observe that minimizing \({C}_{{{{{{\rm{V}}}}}}2}\) leads to a reduction in the capacitance of OECT T2, which can also be achieved by adjusting the geometrical parameters. Supplementary Fig. 10 demonstrates that \({f}_{{{{{{\rm{spike}}}}}}}\) is enhanced when \({W}_{2}\), \({L}_{2}\), and \({t}_{2}\) are minimized. Faster spiking frequencies give rise to a smaller amplitude of oscillations and vice-versa, as confirmed in Supplementary Figs. 1114. As expected, the OAN power consumption (\({P}_{{{{{{\rm{OAN}}}}}}}\)) increases with the spiking frequency (Supplementary Figs. 15 and 16). By contrast, as displayed in Fig. 4 the energy per spike \({E}_{{{{{{\rm{spike}}}}}}}={P}_{{{{{{\rm{OAN}}}}}}}/{f}_{{{{{{\rm{spike}}}}}}}\) decreases with increasing the spiking frequency, indicating that the OAN efficiency is enhanced at faster oscillations. A minimum energy per spike equal to 1 × 10−6 J is obtained by minimizing only \({C}_{{{{{{\rm{V}}}}}}2}\), while keeping all the other parameters at the nominal value. We note that energy per spike of the OAN is significantly higher than the energy consumed by biological neurons46,47 and point out that future work should focus on this direction. In Supplementary Figs. 17 and 18, the analysis of the impact of material parameters on the OAN performance is further extended to several OMIECs suitable for depletion-mode (transistor T1) and accumulation mode (transistor T2) devices. The analysis confirms that OMIECs with large volumetric capacitance and low mobility (e.g., p(gNDI-g2T) in Supplementary Fig. 17), provide OANs with very limited spiking frequency, large energy per spike, and low power consumption. This is further confirmed in Supplementary Fig. 18 where the material with the lowest mobility, i.e., p(gBDT-g2T), provides OANs with minimum spiking frequency and maximum power consumption.

Fig. 4: OAN spiking frequency and energy.
figure 4

Spiking frequency (\({f}_{{{{{{\rm{spike}}}}}}}\), left axis and circle symbols) and energy per spike (\({E}_{{{{{{\rm{spike}}}}}}}\), right axis and triangle symbols) as a function of the OAN parameters. \({E}_{{{{{{\rm{spike}}}}}}}={P}_{{{{{{\rm{OAN}}}}}}}/{f}_{{{{{{\rm{spike}}}}}}}\), where \({P}_{{{{{{\rm{OAN}}}}}}}\) is the power consumption. \({f}_{{{{{{\rm{spike}}}}}}}\) and \({E}_{{{{{{\rm{spike}}}}}}}\) are computed as a function of (a) threshold voltage of T1 \({V}_{{{{{{\rm{TH}}}}}}1}\), b threshold voltage of T2 \({V}_{{{{{{\rm{TH}}}}}}2}\), c normalized transconductance of T1 \({g}_{{{{{{\rm{m}}}}}}1}\), d normalized transconductance of T2 \({g}_{{{{{{\rm{m}}}}}}2}\), e volumetric capacitance of T1 \({C}_{{{{{{\rm{v}}}}}}1}\), f volumetric capacitance of T2 \({C}_{{{{{{\rm{v}}}}}}2}\), g resistance \({R}_{1}\), and h resistance \({R}_{2}\).

Excitability and noise-induced activity

The ability to generate spikes, named excitability, is a fundamental property of neurons that allows them to transmit and process information within the nervous system. Neurons are classified into different categories based on their level of excitability. The excitability of neurons expresses their tendency to elicit spikes when receiving incoming signals or stimuli. Highly excitable neurons have a low threshold for generating action potentials, and even small changes in their membrane potential can trigger an action potential. Some neurons require stronger or more prolonged stimuli to reach the threshold for firing an action potential. These neurons are less sensitive to small fluctuations in their membrane potential and may have a higher threshold for activation.

Taking advantage of the functionality of the OAN, artificial neurons with varying levels of excitability can be achieved. For this purpose, a sinusoidal excitation signal is injected into the electrolyte medium of T1. As depicted in Fig. 5a, a slight increase of a few millivolts in the potential (\({V}_{{{{{{\rm{exc}}}}}}}\)) of the electrolytic medium, corresponding to the same range as extracellular biopotentials48, modulates the excitability of the OAN. During excitation, a phase-locked bursting activity is observed. Furthermore, Fig. 5a demonstrates a high level of agreement between the measurements (gray lines) and numerical simulations (red lines). The numerical simulations reveal that the injected electrolyte voltage \({V}_{{{{{{\rm{E}}}}}}}(t)\) gives rise to a time-varying modulation of the OEND-load line crossing point (Fig. 5b) equal to \(\left[{V}_{{{{{{\rm{IN}}}}}}}+{V}_{{{{{{\rm{exc}}}}}}}+{V}_{{{{{{\rm{E}}}}}}}\left(t\right)\right]/{R}_{{{{{{\rm{L}}}}}}}\). When the crossing point is positioned outside the spiking region (Fig. 5c, thick dark purple segment \(U(t) \, < \, {U}_{{{{{{\rm{ON}}}}}}}\)), a constant bias voltage across CL is obtained during the whole period of the excitation signal and, consequently, the OAN is silent (Fig. 5a top panel). To probe the OAN excitability, \({V}_{{{{{{\rm{exc}}}}}}}\) is increased by a few millivolts, and during the timeframe where the injected signal moves the crossing point \({{{{{\bf{U}}}}}}\) in the negative resistance region (Fig. 5d, thick bright purple segment \(U(t) > {U}_{{{{{{\rm{ON}}}}}}}\)) firing activity is triggered. Excitation can be further enhanced by increasing \({V}_{{{{{{\rm{exc}}}}}}}\) (Fig. 5e), and the bursting width is phase-locked with the electrolyte signal (Fig. 5a), which emulates a fundamental feature of biological neurons49.

Fig. 5: Excitability and phase-locked bursting.
figure 5

a Measured (dot-dashed lines) and calculated (full line) spiking current as a function of time at various input excitation signals \({V}_{{{{{{\rm{exc}}}}}}}\) = 0–40 mV, and \({V}_{{{{{{\rm{E}}}}}}}(t)\) is a sinusoidal signal with amplitude 20 mVpp and frequency 5 Hz. b Current-voltage characteristics \({I}_{{{{{{\rm{OEND}}}}}}}\left({V}_{{{{{{\rm{OEND}}}}}}}\right)\) of the OEND (blue full line) and load-line (dashed line) highlighting the impact of the excitation signal. The excitation signal is depicted by a sinusoidal full line superimposed on the load line. c Zoom showing the range where the crossing point is modulated by the excitation signal (thick dark purple line). \({V}_{{{{{{\rm{exc}}}}}}}\) = 0 V, the crossing point is in the silent region and OAN is not spiking. d \({V}_{{{{{{\rm{exc}}}}}}}\) = 10 mV, the crossing point is modulated by the excitation signal as a function of time. The range where the crossing point is in the spiking region and silent region is marked with a thick brighter and a darker purple line, respectively. e \({V}_{{{{{{\rm{exc}}}}}}}\) = 40 mV, the crossing point is in the spiking. The thick bright purple line shows the range where the crossing point is in the spiking. f Spiking current as a function of time at various input excitation signals \({V}_{{{{{{\rm{exc}}}}}}}\) = 0–40 μV.

The degree of excitability can be adjusted to the desired level by appropriately configuring the OAN parameters. To explore this important aspect, we conducted numerical simulations, systematically varying the material and device parameters. To evaluate the impact of each parameter, we varied one parameter at a time while keeping all others constant. Supplementary Fig. 19 depicts the excitation threshold voltage, \({V}_{{{{{{\rm{exc}}}}}}}\), as a function of OAN parameters. Excitability increases with \({g}_{{{{{{\rm{m}}}}}}1},\,{g}_{{{{{{\rm{m}}}}}}2},\,{C}_{{{{{{\rm{V}}}}}}1},\,{C}_{{{{{{\rm{V}}}}}}2},\,{R}_{1}\) and \({R}_{2}\) while decreases with \({V}_{{{{{{\rm{TH}}}}}}1}\) and \({V}_{{{{{{\rm{TH}}}}}}2}\). These findings suggest that the level of excitability can be finely tuned, ranging from a few microvolts to hundreds of millivolts. This remarkable degree of tunability is achieved by configuring the OAN parameters and allows, for example, phase-locked bursting with excitation input signals as low as a few microvolts. (Fig. 5f).

Modulation of the OAN firing threshold allows to emulate other relevant features of biological neurons, as for example noise-induced activity. Noise, in the context of neurons and neural systems, plays a role in enhancing sensory processing, promoting robustness, generating variability, and shaping network dynamics. Noise couples with neuronal dynamics affecting the neuronal firing properties. To emulate extracellular noise fluctuations, white noise is applied to the electrolyte. The noise-induced transition from tonic to bursting of the OAN activity is shown in Fig. 6. Injected noise sums up to \({V}_{{{{{{\rm{IN}}}}}}}\) and the OAN crossing point dynamically shifts with the amplitude and time evolution of the noise. In the case of a small amplitude noise (e.g., \({V}_{{{{{{\rm{pp}}}}}}}\) = 25 mV, Fig. 6a), a tonic firing is displayed. Numerical simulations allowed us to calculate the \({I}_{{{{{{\rm{OEND}}}}}}}\)-\({V}_{{{{{{\rm{OEND}}}}}}}\) characteristic as well as the load-line crossing point as a function of noise. As displayed in Fig. 6b (gray area around the dashed load line) the amplitude of noise is not sufficient to shift the crossing point outside the spiking region. The OAN is resilient to the noise in the liquid environment and therefore regular firing is obtained. By increasing the level of noise (e.g., \({V}_{{{{{{\rm{pp}}}}}}}\) = 50 mV, Fig. 6c), there is a time-varying probability that the crossing point is shifted in the silent region (Fig. 6d) and in correspondence of such stochastic events OAN spiking is suppressed. Further increasing the noise amplitude (e.g., \({V}_{{{{{{\rm{pp}}}}}}}\) = 150 mV, Fig. 6e), the stochastic modulation of the OAN firing threshold is enhanced (Fig. 6f) and random packets of spikes are triggered by noise. Interestingly, as illustrated in Supplementary Fig. 20, it is possible to design OANs with spiking activity induced by noise: OAN is initially silent (below the excitation threshold) and noise can trigger the spiking activity at various intensities – ranging from a few spikes to nearly tonic firing – by increasing its excitability. This shows the high degree of reconfigurability inherent in the OAN.

Fig. 6: Noise-induced activity.
figure 6

White noise signal with amplitude \({V}_{{{{{{\rm{pp}}}}}}}\) is injected into the electrolyte. a, c, e Spiking current \({I}_{{{{{{\rm{spike}}}}}}}\) as a function of time at various amplitudes (\({V}_{{{{{{\rm{pp}}}}}}}\)) of noise signal injected into the electrolyte. b, d, f Impact of the noise on the current-voltage characteristics \({I}_{{{{{{\rm{OEND}}}}}}}\left({V}_{{{{{{\rm{OEND}}}}}}}\right)\) of the OEND (blue full line) and load-line (dashed line) calculated with numerical simulations. The noise modulates the position of the load line as a function of time (shade area) and elicits interrupted spikes. a, b \({V}_{{{{{{\rm{pp}}}}}}}\) = 25 mV, c, d \({V}_{{{{{{\rm{pp}}}}}}}\) = 50 mV, and e, f \({V}_{{{{{{\rm{pp}}}}}}}\) = 150 mV. Random packets of spikes are triggered by noise.

Neuromorphic ion sensing

Ions are fundamental biological and physiological regulators for vital processes in every living organism. Ion regulation in living organisms involves intra and extracellular fluctuations from the resting condition, and small deviations from the optimal equilibrium levels can be associated with pathological states50. The intracellular-extracellular physiological concentration in mammalian cells is 3–160 × 103 M for potassium (K+) and 10–150 × 103 M for sodium (Na+) ions51. Under resting conditions, neurons maintain a steep gradient between low intracellular free calcium (Ca2+) concentration (0.1–0.5 × 106 M) and high extracellular Ca2+ levels (1 × 103 M)52. In biological neurons, the changes in ion concentration between the intracellular and extracellular medium modulate their excitability and firing threshold53.

The OAN, operating in a liquid environment, displays ion-concentration-dependent spiking properties. In biointerfacing scenarios involving biological cells and OANs, the shared electrolyte serves as a common extracellular space. Information exchange and signaling, such as ionic species transmission, between the biological and artificial domains occur through this shared extracellular medium. Figure 7a shows the measured (symbols) OAN output current (\({I}_{{{{{{\rm{spike}}}}}}}\)) as a function of time by varying the ion concentration in the external medium, i.e., the artificial extracellular space. The OAN excitability and spiking frequency are controlled by the ion concentration of the electrolytic medium. More in detail, at Na+ concentrations below the physiological limit, oscillations are suppressed. Increasing the ion concentration firing is triggered and when the extracellular concentration is varied in the physiological range, \({f}_{{{{{{\rm{spike}}}}}}}\) is modulated from 20 Hz to 50 Hz. The measurements (symbols, Fig. 7a) are accurately predicted by the simulations (solid lines) in the whole range of ion concentration, showing that the numerical framework adequately accounts for the ionic–electronic interaction in the liquid environment. Then, as detailed in Supplementary Note 1, the impact of the ion concentration is accounted for. The spiking frequency as a function of Na+ calculated with the analytical model (Supplementary Note 3) is displayed in Fig. 7b. The results obtained with the computationally inexpensive analytical model are in very good agreement with both the measurements and numerical simulations. The model reveals that the modulation of the spiking frequency is inherently related to the electrochemical properties of the OECTs connected in the OAN configuration. More in detail, the potential drop at the gate/electrolyte interface and the Donnan’s potential at the electrolyte/polymer interface40 are mirrored in a variation of the threshold voltage of T1 (\({V}_{{{{{{\rm{TH}}}}}}1}\)). According to Eqs. (1) and (2), \({V}_{{{{{{\rm{TH}}}}}}1}\) modulates the \({V}_{{{{{{\rm{ON}}}}}}}\) and \({V}_{{{{{{\rm{OFF}}}}}}}\) voltages, which set the maximum charging and discharging voltage of CL, respectively, and eventually define the integration domain in Eq. (5), viz. \({f}_{{{{{{\rm{spike}}}}}}}\). As displayed in Supplementary Fig. 21, the \({I}_{{{{{{\rm{D}}}}}}}\)\({V}_{{{{{{\rm{G}}}}}}}\) characteristics of the OECT systematically shift to more negative voltages with increasing ion concentrations. More negative \({V}_{{{{{{\rm{TH}}}}}}1}\) reduces both \({V}_{{{{{{\rm{ON}}}}}}}\) and \({V}_{{{{{{\rm{OFF}}}}}}}\), thus resulting in a faster spiking activity. In Fig. 7c, d the analysis is extended to other ions, namely K+ and Ca2+, confirming the excitability, firing threshold, and frequency-modulation features of the OAN. Supplementary Figs. 22 and 23 show the corresponding spiking current as a function of time. A modulation of the spiking frequency equal to 38.95 Hz dec1 and 26.69 Hz dec1 is obtained with monovalent and divalent ions, respectively, which correspond to a frequency modulation of about \(\Delta {f}_{{{{{{\rm{spike}}}}}}}/{f}_{{{{{{\rm{spike}}}}}}}=134\,{{{{{\boldsymbol{\%}}}}}}\) for K+ and Na+, and \(\Delta {f}_{{{{{{\rm{spike}}}}}}}/{f}_{{{{{{\rm{spike}}}}}}}=\,\)67% in the case of Ca2+.

Fig. 7: Neuromorphic ion sensing.
figure 7

a Measured (symbols) and modeled (lines) spiking current \({I}_{{{{{{\rm{spike}}}}}}}\) as a function of time. The OAN is exposed to various Na+ concentrations. When the ion concentration is below the minimum physiological level, no spiking is observed. When the ion concentration is above the excitation threshold, OAN spiking is obtained and its spiking frequency is modulated with \({c}_{{{{{{\rm{N}}}}}}{{{{{{\rm{a}}}}}}}^{+}}\). b Spiking frequency \({f}_{{{{{{\rm{spike}}}}}}}\) as a function of Na+ concentration. Symbols are the numerical simulations and the full line is calculated with the analytical model in Supplementary Eq. (26). c \({f}_{{{{{{\rm{spike}}}}}}}\) as a function of K+ concentration. d \({f}_{{{{{{\rm{spike}}}}}}}\) as a function of Ca2+ concentration. In bd the vertical dashed line shows the minimum physiological concentration of each ion. OAN spiking is triggered by the specific ion concentration. e Spiking frequency as a function of ion concentration without and with noise. f Ion-selective spiking: spiking frequency as a function of selected K+ concentration. The OAN is not spiking when the concentration of interfering ions (Na+) is varied. g OEND characteristics \({V}_{{{{{{\rm{OEND}}}}}}}\left({I}_{{{{{{\rm{OEND}}}}}}}\right)\) modulated by the ion concentration (solid lines). Load line (dashed dark line) crosses the OEND characteristics in various points inside the spiking region (orange area), depending on ion concentration. At large ion concentration, the crossing point \({{{{{\bf{U}}}}}}=({V}_{{{{{{\bf{U}}}}}}},\,{I}_{{{{{{\bf{U}}}}}}})\) is outside the spiking region and the OAN is silent. h Crossing voltage \({V}_{{{{{{\bf{U}}}}}}}\) (cross symbols) and current \({I}_{{{{{{\bf{U}}}}}}}\) (plus symbols) as a function of ionic concentration \(c\). \({V}_{{{{{{\bf{U}}}}}}}\) decreases with \(c\) while \({I}_{{{{{{\bf{U}}}}}}}\) increases. In the OAN excitability region (orange region) both \({V}_{{{{{{\bf{U}}}}}}}\) and \({I}_{{{{{{\bf{U}}}}}}}\) logarithmically increase with ion concentration. Outside the excitability region, the OAN is silent.

A key distinctive capability of neuromorphic circuits is their intrinsic robustness against interfering signals. To this aim, noise is intentionally injected into the electrolyte of the OAN to observe its effect on the spiking behavior at different ion concentrations (Supplementary Fig. 24). While the amplitude of the spiking signal is minimally influenced by noise, the information encoded in the spiking frequency remains unchanged. This is evident in Fig. 7e, where the spiking frequency as a function of ion concentration shows nearly identical curves regardless of the presence of noise. To further illustrate this point, we compare our neuromorphic ion-sensing approach with conventional methods based on OECTs54,55,56. When noise of the same amplitude and spectral characteristics as that used in the OAN case is injected into the electrolyte of the OECT, significant differences emerge as shown in the transfer characteristics (\({I}_{{{{{{\rm{D}}}}}}}\)-\({V}_{{{{{{\rm{G}}}}}}}\)) as a function of ion concentration (Supplementary Fig. 25). Unlike the OAN, the OECT drain current is substantially impacted by noise, leading to corruption of information regarding ion concentration. This comparison underscores the intrinsic robustness of neuromorphic approaches, which encode information in the frequency domain, over conventional methods.

Adding another layer of biophysical realism to the OAN response involves integrating selectivity to specific ions, thus obtaining ion-selective excitation and spiking. OANs capable of exhibiting ion-specific oscillatory activities are realized by integrating ionophore-based selective membranes57 (ISM) at the interface between the channel and electrolyte of T1. The ISM generates a voltage at the membrane/electrolyte interface (\({V}_{{{{{{\rm{ISM}}}}}}}\)) in response to the specific ion concentration within the electrolyte (Supplementary Note 1). The integration of ISM in the OAN architecture results in a variation of the OECT T1 threshold voltage \({V}_{{{{{{\rm{TH}}}}}}1}\) (Supplementary Fig. 26). Being the OAN excitability dependent on \({V}_{{{{{{\rm{TH}}}}}}1}\) (Supplementary Fig. 19) variation of the target ion concentration results in ion-selective excitability and spiking of the OAN. Figure 7f shows that a K+ selective OAN exhibits oscillations when immersed in a KCl electrolyte, with the frequency of oscillations increasing with the concentration of the selected ion type. As a control experiment, Na+ is used as interfering ions and the OAN is silent (Fig. 7f), demonstrating ion-selective spiking activity. Analogously, a Na+ selective OAN exhibits oscillations when immersed in a NaCl electrolyte and is silent in KCl (Supplementary Fig. 27).

The ion modulation of the OAN spiking features can be rationalized by analyzing the OEND characteristics as a function of the ion concentration. For the sake of clarity, from here on we considered a variation of the ion concentration in the electrolyte of the transistor T1, and analogous considerations hold for T2. The numerical simulations displayed in Fig. 7g show the OEND characteristics calculated as a function of the ion concentration. The \({I}_{{{{{{\rm{OEND}}}}}}}\)-\({V}_{{{{{{\rm{OEND}}}}}}}\) characteristics are almost perfectly overlapped when the OEND operates in the positive resistance regions because at small currents \({I}_{{{{{{\rm{OEND}}}}}}}\, \approx \,V/({R}_{1}+{R}_{2})\), at large currents \({I}_{{{{{{\rm{OEND}}}}}}}\, \approx \, V/{R}_{1}\) and both \({R}_{1}\) and \({R}_{2}\) are not affected by the ion concentration. Conversely, the OEND characteristics in the negative differential resistance region are significantly modulated by the ion concentration and the crossing point with the load-line characteristic (dashed line in Fig. 7g) systematically shifts with the ion concentration. Specifically, at small ion concentrations, the operating point of the OAN (i.e., OEND-load-line crossing point) lies in the non-excitable region and the OAN is silent. By increasing the ion concentration, a transition of the crossing point in the negative resistance region is obtained (orange region) and the OAN is firing. The crossing point consistently shifts to smaller voltages and larger currents with the increasing of the ion concentration, which enables faster charging and discharging of the load capacitor CL and results in faster OAN oscillations. This behavior is displayed up to a maximum ion concentration where \({V}_{{{{{{\rm{ON}}}}}}}={V}_{{{{{{\rm{OFF}}}}}}}\) (see Supplementary Note 4). Reaching this limit, the OEND characteristic does not show a negative resistance anymore, and therefore the OAN spiking region is lost. The excitability region as a function of the ion concentration is displayed in Fig. 7h. OAN excitability requires a crossing point \(\left({V}_{{{{{{\bf{U}}}}}}},\,{I}_{{{{{{\bf{U}}}}}}}\right)\) located in the negative resistance region of the OEND. The OAN bias voltage \({V}_{{{{{{\bf{U}}}}}}}\) systematically reduces with the ion concentration, while the bias current \({I}_{{{{{{\bf{U}}}}}}}\) increases. According to the previous analysis, larger bias currents and smaller voltages result in enhanced \({f}_{{{{{{\rm{spike}}}}}}}\) when the OAN is operated in the excitability region. Finally, as a convenient design approach, we found that the minimum ion concentration triggering OAN oscillations can be defined with the load-line parameters while the maximum ion concentration of the excitability region can be defined with the concentration-dependent OEND parameters.

Biointerfacing

Direct interfacing with biology has been demonstrated by incorporating a biological membrane in the OAN34. Biological membranes are cellular barriers that control the transport of ions, small molecules, and nutrients through the separated compartments of tissue, thus regulating essential physiological functions in humans, animals, and plants. The transcellular and paracellular pathways are finely regulated by the cellular barriers and the tight-junction proteins allow intercellular sealing and control the paracellular fluxes. A biohybrid OAN is obtained by integrating a cellular barrier between the channel and the gate of T1. Cell medium is used as an electrolyte and the cellular barrier separates the electrolyte into two compartments. Therefore, there is no physical contact between the membrane and the electrode or the device channel. Coupling between the cellular barrier and the OAN components is obtained by means of the cell medium, which represents a shared extracellular space34,58,59. The circuit diagram of the biohybrid OAN is displayed in Fig. 8a. A detailed view of the cellular barrier interface with the transistor T1 is given in Supplementary Fig. 28. We used a non-polarizable Ag/AgCl gate electrode immersed in the apical compartment, which set the potential of the electrolyte to approximately the applied VG42,60,61. We also note that, in general, it is not straightforward to compare qualitatively forced (external) voltages with biopotentials, because the actual voltage that a sensing unit (such as the OAN) experiences is also related to the impedance of a voltage source. For instance, in the case of a source of biological activity (a cell), the potential would depend on the cell-to-device impedance. As detailed in Fig. 8a, the gate electrode is modeled as a resistor, \({R}_{{{{{{\rm{G}}}}}}}\), in series to the ionic resistance of the electrolyte in the apical compartment, \({R}_{{{{{{\rm{Ea}}}}}}}\), due to the ion transport from the gate to the cellular barrier in the electrolyte medium of the apical compartment. The cellular barrier is described by a capacitor \({C}_{{{{{{\rm{M}}}}}}}\) in parallel to a resistor \({R}_{{{{{{\rm{M}}}}}}}\), and the ionic transport in the electrolyte medium of the basal compartment is modeled by the resistor \({R}_{{{{{{\rm{Eb}}}}}}}.\) Electrochemical impedance spectroscopy (EIS) measurements in Supplementary Fig. 29 show that the impedance of intact cellular barrier is much larger than the gate and ionic resistance of apical and basal electrolytes, \(\left|{Z}_{{{{{{\rm{M}}}}}}}\right| \, \gg \, {R}_{{{{{{\rm{G}}}}}}}+{R}_{{{{{{\rm{Ea}}}}}}}+{R}_{{{{{{\rm{Eb}}}}}}}\). Focusing on the cellular barrier model \({C}_{{{{{{\rm{M}}}}}}}\) accounts for the ion accumulation at the basal and apical barrier interfaces. \({C}_{{{{{{\rm{M}}}}}}}\) is relevant when the TJs are intact since in this condition the ions are not flowing through the barrier. \({R}_{{{{{{\rm{M}}}}}}}\) accounts for the ionic transport through the barrier and it is relevant when TJs are disrupted with toxins or transiently open with neuromodulators58. The barrier parameters are obtained from EIS measurements (Supplementary Fig. 29). It is important to observe that the complexity of the equivalent circuit required may vary depending on the types of cells interfacing with the OAN62. The OAN spiking current as a function of time and barrier status is displayed in Fig. 8b, c. Oscillations with a spiking frequency of 23 Hz are obtained without the cellular barrier (blue line). Upon barrier insertion the oscillations are damped, and OAN firing activity is completely suppressed after about 2 s (green line, Fig. 8b, c). The addition of toxin compounds in the cellular medium gives rise to the opening of the TJs and barrier disruption eventually results in a restored spiking of the OAN (pink line, Fig. 8b, c). Importantly, Supplementary Fig. 30 shows that the toxic compounds do not affect the device performance.

Fig. 8: In-liquid biointerfacing.
figure 8

a OAN circuit including the cellular barrier between the apical electrolyte (point A) and the basal electrolyte (point B). The model of the cellular barrier includes the gate \({R}_{{{{{{\rm{G}}}}}}}\) resistance, the apical electrolyte equivalent resistance \({R}_{{{{{{\rm{Ea}}}}}}}\), and the basal electrolyte equivalent resistance \({R}_{{{{{{\rm{Eb}}}}}}}\), the ion resistance \({R}_{{{{{{\rm{M}}}}}}}\) and capacitance \({C}_{{{{{{\rm{M}}}}}}}\) of the membrane, and the channel resistance \({R}_{{{{{{\rm{C}}}}}}}\). b Measurements of the OAN without cellular barrier (blue line), with intact cellular barrier interfaced with the OAN (green line), and with disrupted cellular barrier (purple line). Insertion and disruption of the cellular barrier is highlighted with the vertical dashed lines. c Transient numerical simulations of the OAN without cellular barrier (blue line), with intact cellular barrier interfaced with the OAN (green line), and with disrupted cellular barrier (purple line). Insertion and disruption of the cellular barrier is highlighted with the vertical dashed lines. d Model of the OAN and cellular barrier accounting for an input signal at the apical electrolyte (point A) and output signal at the basal electrolyte (point B). e Calculated transfer function module between A-B \(\left|{H}_{{{{{{\rm{AB}}}}}}}\right|\) for the frequency range 1 Hz–105 Hz. The signal is attenuated in a frequency range relevant to the oscillation (range 10–150 Hz). f Impact of the membrane resistance on the OAN oscillations. g Impact of the membrane capacitance on the OAN oscillations. h Time constant of the oscillation decay \(({\tau })\) after the cellular barrier insertion as a function of RM.

We analyzed the excitability modulation of the OAN due to the cellular barrier by means of numerical simulations. The \(I\)-\(V\) characteristics of the OEND and load line are displayed in Supplementary Fig. 31. The crossing point between the OEND and load line is located in the spiking region but, in contrast to the previous cases, it is neither affected by the barrier insertion nor by the barrier status. Interestingly, Fig. 8b, c show that the OAN output is damped after the insertion of the barrier (green line). This suggests that the spiking activity is dynamically attenuated by the barrier. To investigate this effect, we performed AC analysis of the OAN. As schematically depicted in Fig. 8d, we injected a small input signal (\({v}_{{{{{{\rm{i}}}}}}}\)) in the apical electrolyte and the output signal (\({v}_{{{{{{\rm{o}}}}}}}\)) is obtained in the basal electrolyte. The frequency (\(f\)) of the input signal is scanned from 1 Hz to 105 Hz and the Laplace transform of the signals provides the transfer function \({H}_{{{{{{\rm{AB}}}}}}}={V}_{{{{{{\rm{o}}}}}}}(s)/{V}_{{{{{{\rm{i}}}}}}}(s)\), where \(s={{{{{\rm{j}}}}}}2{{{{{\rm{\pi }}}}}}f\). Figure 8e shows that the signal is attenuated (\(|{H}_{{{{{{\rm{AB}}}}}}}|\) < −0.2 dB) by the cellular barrier in the range of frequency relevant for the OAN oscillation (10–150 Hz). Therefore, at each cycle, an attenuated spike is fed back to the OAN resulting in an overall dumping of the oscillation. The impact of the membrane parameters \({R}_{{{{{{\rm{M}}}}}}}\) and \({C}_{{{{{{\rm{M}}}}}}}\) is analyzed in Fig. 8f, g, respectively. Figure 8f shows that spike attenuation takes place when \({R}_{{{{{{\rm{M}}}}}}}\) > 103 Ω and it progressively increases with the barrier resistance. We note that barrier resistance depends on the cellular barrier status. Conversely, Fig. 8g shows that \({C}_{{{{{{\rm{M}}}}}}}\) has a minor effect on the magnitude of the spike attenuation and larger \({C}_{{{{{{\rm{M}}}}}}}\) extends the frequency range where barrier attenuation is effective. As a result, the excitability of the OAN is controlled by the barrier status (Supplementary Fig. 32) and information on the TJs modulation, mirrored in the barrier resistance, is obtained from the damping time constant \({\tau }\) of the spiking output (Fig. 8h).

Discussion

By combining experiments, new numerical tools, and analytical expressions, we unrevealed the operation of organic electrochemical artificial neurons. Starting from the investigation of the OEND, which is the core of the OAN, we systematically rationalized the OAN fundamental operations focusing on spiking frequency, voltage and current amplitude of the output oscillations, power consumption, and energy per spike. Numerical simulations accurately predicting the measurements have shown that such OAN features intimately depend on the material and device parameters. Table 1 highlights the impact of the various material, biochemical, geometrical, and device parameters on the OAN performance, providing general guidelines to optimize and engineer the OAN and to shape the spiking profile. An overview of artificial spiking neurons providing biophysical realism capabilities (viz. emulation of real neuron spiking behavior) is displayed in Table 2. The wide biorealistic repertoire of the OAN including excitability, noise-induced firing activity, neuromorphic ion sensing, and in-liquid bio-interfacing has been rationalized. We have shown that excitability threshold and noise-induced firing activity are inherently due to the organic electrochemical non-linear device and the key design parameters have been identified and discussed. We demonstrated that the OAN excitability and spiking frequency controlled by the ion concentration of the electrolytic medium can be precisely designed to reproduce the behavior of real neurons. A biohybrid OAN obtained by integrating a cellular barrier has been investigated by means of transient non-linear simulations. The impact of the cellular barrier on the OAN characteristics has been quantified. The analytical expressions describing key aspects of neuronal behavior agree with both measurements and numerical simulations and link the OAN spiking features and its materials/physical parameters, thus bringing closer the domains of artificial neurons and neuroscience. This work provides streamlined and transferable guidelines for the design, development, engineering, and optimization of organic artificial neurons, pushing forward a community-aware and community-wide development of the next-generation OAN-based neural networks, neuromorphic electronics, and intelligent bioelectronics.

Table 1 OAN performance as a function of the material and device parameters
Table 2 Overview of artificial spiking neurons

Methods

Device/OAN fabrication

Standard microscope glass slides (75 mm × 25 mm) were cleaned in a sonicated bath, first in a soap solution (Micro-90 (Sigma–Aldrich)) and then in a 1:1 (vol/vol) solvent mixture of acetone and isopropanol. Source and drain electrodes were made with photolithographically patterned gold (with positive Microposit S1813 photoresist (DOW)) on the cleaned glass slides. A chromium layer was used to improve the adhesion of gold. Each glass slide contains two OECTs T1 and T2, with connections to their respective contacts. The channel dimensions of T1 and T2 are \({W}_{1}\) × \({L}_{1}\) = 50 µm × 20 µm and \({W}_{2}\) × \({L}_{2}\) = 50 µm × 10 µm, respectively. A separate Ag/AgCl electrodes was used for separately gating the OECTs via aqueous electrolytes. Two layers of parylene C (SCS Coatings) were deposited. Soap (Micro-90 soap solution, 1% vol/vol in deionized water) was used for separation between the parylene C layers to enable the peel-off of the upper parylene C layer. The lower parylene C layer insulates the gold electrodes. Silane A-174 (γ-methacryloxypropyl trimethoxysilane) from Sigma–Aldrich was added to the lower parylene C layer to enhance adhesion. The channel dimensions of T1 and T2 were defined in the second photolithography step through the positive photoresist AZ 9260 MicroChemicals (Cipec Spécialités). Reactive ion etching (O2/CF4 plasma, 160 W for 16 min with O2 flow rate of 50 s.c.c.m. and CHF3 flow rate of 5 s.c.c.m.) was used to define the channels of T1 and T2 throughout the photoresist mask. T1 channel is made with the organic mixed ionic–electronic conductor polymer PEDOT:PSS (Clevios PH 1000) mixed with 5.0 wt% ethylene glycol, 0.1 wt% dodecyl benzene sulfonic acid, and 1.0 wt% (3-glycidyloxypropyl)trimethoxysilane. Spin coating was used to produce a film in two steps at 1500 rpm and 650 rpm for 1 min and annealed at 120 °C for 1 min in between. T2 channel is made with the semiconducting polymer p(g2T-TT), synthesized according to another work63. Here p(g2T-TT) was dissolved in chloroform (3 mg ml–1) inside an N2-filled glovebox and spin-coated in ambient conditions at 1000 rpm for 1 min resulting in a thickness of 40 nm. The devices were baked at 60 °C for 1 min. The sacrificial upper parylene C layer was peeled off to confine the polymer to the inside of the channel regions. The devices were subsequently baked at 140 °C for 1 h. Excess soap was rinsed off with deionized water. A diagram of the device fabrication process is provided in Supplementary Fig. 33.

Numerical simulations

The OECT model reported in Supplementary Note 1 was implemented in a Verilog-A module and used for OAN simulations in the electronic design automation software Advanced Design System (ADS) 2023 Update 1 by Keysight Technologies. The file containing this model was included in a “Verilog” folder inside the ADS 2023 project folder, then a symbol was created and linked to the symbol. The symbol must have the same input/output ports as defined in the Verilog-A module. The user must define the same internal parameters with corresponding names (case sensitive). The symbol created can then be used in any valid circuit topology. The model parameters obtained by reproducing the measured transfer characteristics of T1 and T2 displayed in Supplementary Fig. 1 read: \({W}_{1}\) = 50 µm, \({L}_{1}\) = 20 µm, \({g}_{{{{{{\rm{m}}}}}}1}\) = 1.5 × 103 S V1, \({V}_{{{{{{\rm{TH}}}}}}1}\) = 0.477 V, \({\gamma }_{1}\) = 2, \({W}_{2}\) = 50 µm, \({L}_{2}\) = 10 µm, \({g}_{{{{{{\rm{m}}}}}}2}\) = 3.4 × 103 S V1, \({V}_{{{{{{\rm{TH}}}}}}2}\) = –0.23 V, \({\gamma }_{2}\) = 2.18. Numerical simulations are performed in DC, AC, transient, and frequency mode. Data analysis of the simulation results was performed by using the software MATLAB 2023a from Mathworks.