Unveiling the capabilities of bipolar conical channels in neuromorphic iontronics

Abstract

Conical channels filled with an aqueous electrolyte have been proposed as promising candidates for iontronic neuromorphic circuits. This is facilitated by a novel analytical model for the internal channel dynamics [T. M. Kamsma, W. Q. Boon, T. ter Rele, C. Spitoni and R. van Roij, Phys. Rev. Lett., 2023, 130(26), 268401], the relative ease of fabrication of conical channels, and the wide range of achievable memory retention times by varying the channel lengths. In this work, we demonstrate that the analytical model for conical channels can be generalized to channels with an inhomogeneous surface charge distribution, which we predict to exhibit significantly stronger current rectification and more pronounced memristive properties in the case of bipolar channels, i.e. channels where the tip and base carry a surface charge of opposite sign. Additionally, we show that the use of bipolar conical channels in a previously proposed iontronic circuit features hallmarks of neuronal communication, such as all-or-none action potentials and spike train generation. Bipolar channels allow, however, for circuit parameters in the range of their biological analogues, and exhibit membrane potentials that match well with biological mammalian action potentials, further supporting their potential biocompatibility.

---

I. Introduction

Iontronics is an exciting emerging platform that aims to harness the transport of ions for use in signalling. The ability to control ion transport in confined geometries offers unique opportunities compared to traditional electronic systems, such as the ability to mimic biological processes or interface with cells and tissues.^1–3 A particularly exciting direction is that of neuromorphic (brain-inspired) iontronic circuits,^1–6 which offer the unique feature of closely resembling biological processes, as signalling in the brain also relies on ion transport.^11,12 Promising candidates for the realisation of such circuits are ionic microfluidic memristors (memory-resistors).^{7–10,13–24} The dynamical properties of memristors make them artificial analogues of biological synapses, the connections between neurons.^25–29 Over the past few years, a vast array of different memristors have been extensively investigated as components for neuromorphic circuit architectures.^51–53 Not only are memristors analogous to synapses, they are also analogous^12,25 to the biological ion channels present in neuronal membranes that facilitate signalling,^11,30–43 offering even more perspectives for brain-inspired circuits. Due to the prospect of more energy-efficient computers^44,45 and biocompatibility,^26,46–50 memristors and neuromorphic circuit architectures have drastically increased in popularity over recent years.^51–53 However, the emphasis has mostly been on memristors that require electrons or holes as charge carriers,^45,51–53 limiting their applicability in fully ionic fluidic systems.

In the past few years, however, some interest has been shown in microfluidic neuromorphic components.^{7–10,13–23} Candidates of specific interest for the present work are conical channels containing an aqueous electrolyte and a homogeneous surface charge, which are known to act as iontronic microfluidic memristors.^13–23 Recently, an analytical model was derived that explains in a quantitative manner how transient concentration polarisation in such channels is responsible for volatile conductance memory and it was demonstrated that these channels could have the potential to be used in experimentally accessible neuromorphic iontronic circuits.⁸ The underlying functionality that underpins the memristive behaviour of conical channels is that they exhibit current rectification in the steady state.^54–58 Although conical channels with a homogeneous surface charge distribution are desirable as relatively simple model systems for investigating iontronic systems,^{55,57,59–72} they are not necessarily the best-performing channels for current rectification.⁷³ In fact, conical channels with a surface charge distribution that changes sign as a function of the distance to the tip are known to exhibit a much stronger current rectification than homogeneously charged ones.⁷⁴ These so-called bipolar conical channels are therefore promising systems to advance the field of iontronic (neuromorphic) systems.

In this work we present an analytical model that quantitatively captures both the steady-state and dynamical behaviour of conical channels with an inhomogeneous surface charge distribution, based on the methodology in ref. 8 and 55. Our model contains no free parameters and can quantitatively predict the steady-state and time-dependent ionic charge currents as a result of static and dynamic applied voltages, respectively. An understanding of these current–voltage relations and the dependence on system parameters could allow for a targeted development of new circuits of these channels and more effortless scanning of possible applications thereof.

Recently, some fully microfluidic circuits that display neuronal behaviour have been theoretically and numerically explored. These circuits, through which an imposed current can be driven, consist of artificial ion channels, batteries, and a capacitor.^7,8 In ref. 7 a circuit was modelled containing quasi two-dimensional nanochannels that connect aqueous electrolytes, which generated a train of voltage spikes, a feature of neuronal communication.^31–35 In ref. 8, a more experimentally accessible circuit containing conical channels with homogeneous surface charge was proposed that also obeys the defining all-or-none law for action potentials.^11,30,31 Here, we will demonstrate that bipolar conical channels can be used in the circuit from ref. 8 to achieve the same all-or-none action potentials and spike trains, but with more biologically relevant salt concentrations and battery potentials. Furthermore, the circuit's membrane potentials during spiking match typical mammalian values, making it more biocompatible.

II. Bipolar conical channel

To study the steady-state properties and the voltage-driven dynamics of conical channels with an inhomogeneous surface charge eσ(x), we consider an azimuthally symmetric long conical channel of length L, base radius R_b ≪ L at x = 0, and tip radius R_t < R_b at x = L, with x being the longitudinal coordinate, as depicted schematically in Fig. 1(a). The channel radius is described by R(x) = R_b − xΔR/L for x ∈ [0, L], the central axis being at radial coordinate r = 0 and ΔR = R_b − R_t. On the channel surface, at r = R(x), we assume an inhomogeneous surface charge distribution eσ(x), with a positive surface charge at the base and middle of the channel, and a negative surface charge at the tip, with a linear increase described bywhere throughout this manuscript we set σ′ = −3σ₀/2 with eσ₀ = 0.1e nm⁻², resulting in a bipolar (BP) channel. Unless stated otherwise, we set the channel dimensions as length L = 10 μm, base radius R_b = 200 nm, and tip radius R_t = 50 nm, resulting in a channel geometry similar to that realised before experimentally.⁷⁵ The channel connects two bulk reservoirs of an incompressible aqueous 1 : 1 electrolyte, with mass density ρ_m = 10³ kg m⁻³, viscosity η = 1.01 mPa s, and electric permittivity ε = 0.71 nF m⁻¹. The electrolyte contains ions with monovalent charges ±e, with e being the proton charge, and diffusion coefficients D_± = D = 2 μm² ms⁻¹, a typical value for dilute KCl⁷⁶ in microfluidic channels.^77,78 At the far side of both reservoirs we impose fixed ion concentrations ρ_± = ρ_b = 2 mM, such that the equilibrium Gouy–Chapman surface potential ψ₀(x) = (2k_BT/e)sinh⁻¹[2πλ_Bλ_Dσ(x)] equals ψ₀(0) ≈ 92 mV at the base and ψ₀(L) ≈ −61 mV at the tip. Here we introduced the Bjerrum length λ_B = e²/(4πεk_BT) and the Debye length . An electric double layer (EDL) forms that screens the surface charge with λ_D ≈ 6.8 nm. The far sides of both reservoirs are kept at a constant and equal pressure P = P₀. On the far side of the reservoir connected to the base we impose an electric potential V(t), while the far side of the other reservoir is grounded. The resulting physical quantities of interest in this system are the electric potential profile Ψ(x,r,t), the ionic concentration profiles ρ_±(x,r,t), an electro-osmotic fluid flow with velocity field u(x,r,t), ionic fluxes j_±(x,r,t) with j₊ − j₋ being the charge flux and j₊ + j₋ the salt flux, and the pressure profile P(x,r,t).

Fig. 1 (a) Schematic representation of the azimuthally symmetric bipolar (BP) conical channel (not to scale), with channel length L, base radius R_b, and tip radius R_t < R_b, connecting two bulk reservoirs of a 1 : 1 aqueous electrolyte, with bulk concentration ρ_b. The channel wall carries a surface charge density eσ(x), with . Here σ′ = −3σ₀/2 with eσ₀ = 0.1e nm⁻² such that the surface charge is positive at the base (eσ(0) = eσ₀) and negative at the tip (eσ(L) = −eσ₀/2). A possibly time-dependent electric potential drop V(t) is applied over the channel, driving an ionic charge current I(t) = g(V(t),t)V(t) with g(V(t),t) being the channel conductance that we calculate in this paper. (b) Steady-state current I as a function of the static potential V as predicted by full FE calculations of the PNPS eqn (2)–(5), for a bipolar (BP) channel (blue) and otherwise identical unipolar (UP) channels with uniform surface charges −σ₀/2 (green) and −σ₀ (red). An applied positive (negative) voltage over the channel results in ion depletion (accumulation) as depicted in the insets of (b), responsible for the steady-state diodic behaviour of the cones.⁵⁵

The aforementioned physical quantities can be described by a coupled set of equations. Firstly, the ionic fluxes j_± and concentration profiles ρ_± satisfy the continuity equation,

and the Nernst–Planck equation,where the three terms account for Fickian diffusion, ohmic conduction, and Stokesian convection, respectively. The fluid flow u(x,r,t) satisfies a force balance described by the Stokes equation for an incompressible fluid,where −eρ_e∇Ψ is the electric body force, which depends on the ionic space charge density ρ_e = ρ₊ − ρ₋. This space charge density affects the electric potential Ψ(x,r,t), which satisfies the Poisson equation,where Ψ(−∞,r,t) = V(t) and Ψ(∞,r,t) = 0. Eqn (2)–(5) form the Poisson–Nernst–Planck–Stokes (PNPS) system of equations. To make the system closed we impose the boundary conditions on the channel wall given by the no-slip condition u = 0, the blocking condition n·j_± = 0, and Gauss' law n·∇Ψ(x,R(x),t) = eσ(x)/ε, with n being the normal vector of the wall.

In Fig. 1(b) we show steady-state current–voltage (I–V) curves as determined using finite-element (FE) calculations, not only for the BP channel under consideration (blue) but also for homogeneous and unipolar (UP) surface charge densities −eσ₀/2 (green) and −eσ₀ (red). An applied potential V over the channel leads to a depletion or accumulation of ions, where for our parameters V < 0 results in salt accumulation while V > 0 depletes the channel of salt, as shown in the insets of Fig. 1(b), thereby changing the channel conductance. This concentration polarisation is responsible for the ion current rectification (ICR) found in conical channels.⁵⁵ It is clear that the BP channel exhibits a significantly stronger current rectification, with the ratio ICR = |I(−0.8 V)/I(0.8 V)| of the current I(V) at voltages V = ±0.8 V being as large as ICR ≈ 21 for the BP channel (blue), but as small as ICR ≈ 3 and 2.4 for the UP channels with surface charges −σ₀/2 (green) and −σ₀ (red), respectively.

III. Analytical approximation for bipolar channels

The full PNPS eqn (2)–(5) cannot be solved analytically for the system of interest here in the steady-state. However, with a few reasonable assumptions we can simplify them to obtain some closed-form analytical descriptions.⁵⁵ Under the assumption that the Debye length is small compared to the channel radius, i.e. λ_D ≪ R(x), we can make the approximation that for all r at least a few λ_D away from the surface, the salt concentration ρ₊(x,r) + ρ₋(x, r) = ρ_s(x,r) ≈ [small rho, Greek, macron]_s(x) and the electric potential Ψ(x,r) ≈ [capital Psi, Greek, macron](x) are radially independent. With this assumption, as in ref. 55, the slab-averaged electric field −∂_x[capital Psi, Greek, macron](x) and the total salt flux j_s(x,r) = j₊(x,r) + j₋(x,r) can be radially integrated to obtain expressions for the cross-sectional averaged electric field,and the total salt flux through the channel,with Q(V) = −πR_tR_bεψ_effV/(ηL) being the electro-osmotic fluid volume flow, which is similar to the expression for the fluid flow of a UP channel⁵⁵ except for the surface potential term ψ_eff. In BP channels it is not immediately clear how the inhomogeneous surface charge dictates the electro-osmotic flow. For our standard parameter set we use ψ_eff ≈ −25 mV as a fit parameter, which will be discussed further in Section IV. Eqn (7) represents the diffusive, conductive and convective components of the salt flux, respectively. In the steady state the condition ∂_xJ(x) = 0 must hold, yielding for a given σ(x), R(x), and Q(V) a differential equation for the unknown radially averaged salt concentration profile function [small rho, Greek, macron]_s(x). In ref. 55 this differential equation is solved for a conical channel with a homogeneous surface charge and boundary conditions [small rho, Greek, macron]_s(0) = [small rho, Greek, macron]_s(L) = 2ρ_b. Here we consider the case where the surface charge distribution is given by eqn (1). By solving ∂_xJ(x) = 0 for a given static potential V we obtain the following expression for the radially averaged salt concentration:with Pe = Q(V)L/(πDR²_t) being the Péclet number at the narrow end. Note that for our case with a solely voltage-driven flow without any pressure-driven contribution, Q(V) = −πR_tR_bεψ_effV/(ηL) is proportional to V, and hence the ratio Pe/V that appears in eqn (8) does not depend on the static potential V. In Section IV we will see that for our case of ψ_eff < 0 a negative applied voltage (V < 0) will cause an enhancement of the ion concentration in the channel (and hence an increased conductivity), whereas a positive potential (V > 0) gives rise to an ionic depletion and a reduced conductivity, where the effect of ion accumulation and depletion becomes stronger upon increasing |V|. For V > 0 we will see that the profile as predicted by eqn (8) can even become negative, which is obviously an unphysical result that stems from a break-down of the λ_D ≪ R(x) assumption that underlies eqn (8). However, we will discuss in Section IV how we can still ensure good agreement in the current–voltage relation over a wide voltage range.

Interestingly, eqn (8) suggests that it is possible to explain and predict current rectification in cylindrical channels^79–81 as long as σ′ ≠ 0, since in this case a non-trivial source term remains in eqn (8) even for ΔR = 0. Hence our current work suggests to provide a unified theory for non-linear transport through cylindrical and conical channels carrying homogeneous or inhomogeneous surface charges. Additionally we note that eqn (8) seems to suggest that bipolar and conical rectification mechanisms can oppose each other, even to the extent that no current rectification is expected if σ₀ΔR = −σ′R_b (which for our linear surface charge density profile implies σ_t/σ_b = R_t/R_b with σ_t and σ_b being the surface charge at the tip and the base, respectively). Probing this unification will be left for future work, while we will focus here on a more constrained parameter set to investigate the iontronic neuromorphic circuit in Section V.

The static electric conductance of the conical channel can be found by treating the concentration profile as a series of resistors of thickness dx and cross-sectional area πR²(x). Since the electric field scales with the inverse of R²(x) according to eqn (6), the contribution to the resistance of each slab equals (g₀[small rho, Greek, macron]_s(x)/(2ρ_b))⁻¹ dx with the homogeneous channel conductance g₀ = (πR_tR_b/L)(2ρ_be²D/k_BT),^55,82 yielding for the static channel conductivity

In order to account for the possibility of unphysical negative concentration profiles at high positive voltages, we replace [small rho, Greek, macron]_s(x,V) with max[0.2ρ_b,[small rho, Greek, macron]_s(x,V)] in the actual (numerical) evaluations of eqn (9), such that we effectively take the surface conductivity into account by not allowing the concentration profile to drop below 10% of the bulk salt concentration 2ρ_b. This ad hoc cut-off can certainly be improved upon, although the details of the cut-off have limited effects for the system parameters that we use and discuss below. The steady-state current is then given by

I(V) = g_∞(V)V. (10)

As we will show in Section IV, eqn (10) predicts a diodic behaviour of the conical channel through ion depletion (and hence a low conductivity) for V > 0 and ion accumulation (and hence a high conductivity) for V < 0.

It was found in ref. 8 that the process of ion accumulation and depletion is not instantaneous and occurs over a diffusion-like timescale. To derive an expression for the timescale of this process and thus the typical memory retention time τ of a BP conical channel from the PNPS eqn (2)–(5), we apply the same methodology. We consider two quantities, the total number of ions and the net salt flux J_x(0) − J_x(L) into the channel. The change of N given by eqn (8) upon a small voltage perturbation V′ around V = 0 yields

where α < 0 for the standard parameter set of our BP channel, in agreement with the enhanced (reduced) conductance of a negative (positive) potential V′.

At V = 0 the concentration profile is at equilibrium, so for a small voltage perturbation V′ we can assume [small rho, Greek, macron]_s(x) = 2ρ_b. With this assumption the first and third terms in eqn (7) vanish. The net salt flux into the channel, J_x(0) − J_x(L), is then determined by the remaining conductive terms

where γ < 0 for our parameter choices. The typical time taken for ion depletion or accumulation, and thus the typical memory retention timescale, is then approximated by τ = α/γ. This yields, perhaps surprisingly, the purely diffusive timescaleidentical to the expression for UP channels,⁸ which is remarkable as the conductive terms in eqn (7), through which eqn (13) is obtained, differ from those in ref. 8. For our standard parameter set we find τ ≈ 4.17 ms. By assuming that ∂_tg(V(t),t) ∝ g_∞(V(t)) − g(V(t),t), an assumption proven to be effective before,^8,9,83 we can describe the time-dependent conductance g(V(t),t) at a given applied voltage V(t) asand the current I(t) as

I(t) = g(V(t),t)V(t) (15)

Despite the fact that eqn (9) needs to be evaluated numerically, we will refer to eqn (8), (10) and (15) as an analytical approximation (AA) for the voltage-dependent salt concentration profiles, steady-state current, and time-dependent current, respectively. In Section IV we will verify these three equations against full FE calculations of the PNPS eqn (2)–(5).

IV. Finite-element verification

In Section III we derived an AA for the voltage-dependent salt concentration profiles, steady-state current, and time-dependent current. Here we will verify these results against full FE calculations of the underlying PNPS eqn (2)–(5). Throughout this section we will use our standard parameter set and vary the applied voltage. Firstly, in Fig. 2(a) we compare for a variety of positive and negative static voltages V the radially averaged concentration profiles as predicted by eqn (8) (solid lines) with those from the FE calculations (circles). For V < 0 we observe ion accumulation and excellent agreement with almost indistinguishable results for AA and FE. For V > 0 the agreement is still reasonable and qualitative; however, a quantitative discrepancy is now clearly visible, especially at larger positive voltages. Whereas the FE concentration profile at the highest voltage (V = 200 mV, purple circles) shows a depletion of salt down to about 30% of the bulk concentration at x/L ≃ 0.75, the FE-generated concentration at this point remains strictly positive, of course. By contrast, the corresponding AA profile (purple line) falls below 10% of the bulk concentration (indicated by the horizontal line) and in fact even becomes negative in the region of x/L ≃ 0.75. As we stated before, the extremely low local salt concentration at high V causes a break-down of the AA-assumption of a small Debye length λ_D (compared to R_t), a problem that we cure in an ad hoc fashion by replacing [small rho, Greek, macron]_s with max[0.2ρ_b,[small rho, Greek, macron]_s(x,V)] in eqn (9).

Fig. 2 Comparisons of finite-element calculations (FE) of the full PNPS eqn (2)–(5) and our analytical approximation (AA) of eqn (8), (10) and (15), all for our standard parameter set of a bipolar conical channel (see text). (a) The radially averaged salt concentration profiles as determined using eqn (8) (solid lines) and the FE calculation (circles) for various static potentials V ∈ [−200, 200] mV as indicated by the colours. (b) Steady-state current–voltage relation as predicted by our AA of eqn (9) and (10) (red) and by the FE calculations (blue), featuring strong (diodic) current rectification. (c) Current–voltage diagram for an applied periodic triangle potential V(t) with amplitudes ± 1 V and frequency f = 45 Hz, revealing a clear pinched hysteresis loop.

In Fig. 2(b) we translate the concentration profiles at static potentials V to the steady-state current–voltage relation I(V) through eqn (9) and (10) (red) and compare them with I(V) as obtained from FE calculations (blue). There is a good agreement; most notably a very similar strongly diodic effect is found through both methods, with quantitatively similar currents. The agreement also seems to hold for strong positive potentials, despite the aforementioned decrease in accuracy of the AA for this voltage regime. We propose that the I–V relation still matches well since the prediction that the channel is locally nearly completely depleted of salt for high static potentials does in fact match with FE calculations. Therefore, replacing [small rho, Greek, macron]_s(x,V) with max[0.2ρ_b,[small rho, Greek, macron]_s(x,V)] in eqn (9) effectively captures, for this parameter set at least, this depletion and ensures an I–V relation agreement over a wider voltage range than perhaps could have been expected. We do note that the circuit we propose in Section V relies on potentials in the range ± 0.2 V, therefore operating on voltages within the AA range of validity.

Lastly, in Fig. 2(c) we plot the current–voltage relation I(t)–V(t) for the case of an applied periodic triangle potential V(t) with amplitudes of ± 1 V and frequency f = 45 Hz, in line with the prediction that τf ≈ 0.19 yields the most pronounced memory effect.⁸ We compare the time-dependent current determined through eqn (15) (red) against FE calculations (blue). In both instances a similar pinched hysteresis loop is found, the hallmark of a memristor.⁸⁴ We note that this hysteresis loop shows a much more pronounced opening compared to the loop of a similar UP channel,⁸ showing that the stronger current rectification of BP channels translates to a stronger memristive effect.

Before we consider iontronic circuits of BP conical channels in Section V, which essentially only involve the AA approximation of the current–voltage relation, let us consider to what extent the radially averaged electric field −∂_x[small psi, Greek, macron](x) and the fluid flow Q(V) are accurately described by our AA for BP channels for various static V. In the AA −∂_x[small psi, Greek, macron](x) is given by eqn (6), which shows good agreement for UP conical channels in the present parameter regime.⁵⁵ In Fig. 3(a) we compare the electric field for various static potentials V as predicted by eqn (6) (solid lines) with FE calculations (circles). For negative and moderately positive potentials we find good agreement, as expected from the agreement we found in Fig. 2(a); however, a clear disagreement is observed for larger positive static voltages V ≳ 0.2 V. As before, we expect this to be due to the strong ion depletion at high positive potentials, typically in the vicinity of x ≈ 2L/3. The resulting overlapping EDLs in combination with the longitudinal dependence of σ(x) create a local buildup of a longitudinally varying ionic charge density, creating a peak in the (no longer divergence-free) electric field around the location of the strongest depletion at x ≈ 2L/3. This explanation relies on the longitudinal electric field within the EDL that is inherently present in BP channels due to the surface charge inhomogeneity; this longitudinal field is not present in UP channels with similar parameters as the surface charge is homogeneous. Moreover, the salt depletion is much weaker in UP channels and thus the underlying assumption that λ_D ≪ R(x) remains valid for a wider voltage regime.^8,55 This is probably why the peak in −∂_x[small psi, Greek, macron](x) for V = 300 mV in Fig. 3(a) is not observed in UP channels in a similar parameter regime.⁵⁵

Fig. 3 (a) Steady-state electric field −∂_x[capital Psi, Greek, macron](x) inside the channel as predicted by eqn (6) (solid lines) and as measured on the central axis of the channel through the FE calculations (circles) of the PNPS eqn (2)–(5) for various static applied potentials V. (b) Steady-state fluid volume flow Q(V) as a function of the static potential V as predicted by Q(V) = −πR_tR_bεψ_effV/(ηL) with ψ_eff = −25 mV a fit parameter for the linear regime (V ≲ 0.4 V) of Q(V) (red) and as determined by FE calculations (blue).

The underlying eqn (8) of the reported result is dependent on the fluid volume flow Q(V) = −πR_tR_bεψ_effV/(ηL), which we show in Fig. 3(b) (red) compared to FE calculations (blue). The relation of fluid flow Q(V) to surface potential ψ_eff is not immediately clear. In UP channels, Q(V) ∝ ψ₀ with ψ₀ being the (homogeneous) surface potential,⁵⁵ but in BP channels such a relation is not obvious as the surface potential ψ₀(x) is inhomogeneous. In Fig. 3(b) we show that using ψ_eff = −25 mV as a fit parameter based on the linear regime of Q(V) (i.e. for V ≲ 0.4 V) yields good agreement (red) with FE calculations (blue) for roughly the same voltage regime as where we find good agreement for the electric field. Fascinatingly, from Fig. 3(b) we conclude for stronger positive potentials V ≳ 0.4 V that the BP channel acts as a fluidic diode. Remarkably, the tip polarity (here negative) determines the direction of the electro-osmotic flow, positive for positive V and negative for negative V, despite the majority of the channel carrying a positive surface charge. Additionally, the strength of ψ_eff = −25 mV also seems to be similar to the average surface potential of the tip . Whether the tip polarity is a general predictor for the strength and direction of the electro-osmotic flow and whether the fluidic diode behaviour emerges for other parameter configurations requires a more extensive investigation of the parameter space. We leave this topic for future studies and instead here focus on our standard parameter set in order to continue with investigating the iontronic neuromorphic circuit in Section V.

We conclude this section by stating that although we find deviations for the salt concentration profiles, electric field profiles, and fluid volume flow for relatively large positive potentials, these deviations seem to have a limited impact on the overall I(V) relations as demonstrated in Fig. 2(b) and (c), which are most relevant in the context of iontronic circuitry. Furthermore, the iontronic circuit presented in Section V operates within a voltage regime where the electric fields and fluid flows predicted by the analytical approximation are reasonably consistent with those from FE calculations.

V. Neuromorphic microfluidic circuit

We proceed to investigate the use of BP channels in iontronic circuits; specifically, we are interested in neuromorphic circuits. In biological systems the process of neuronal signaling is enabled by the transport of various ionic species through the neuronal cell membrane. Upon application of a stimulus of sufficient strength and duration, a process is set in motion which results in a voltage spike over the membrane due to modulated ionic charge transport through biological ion channels. Such voltage spikes are known as action potentials (APs) and follow the characteristic all-or-none law, meaning that the membrane does not spike at all for stimuli below a critical threshold.^11,30,31 Neurons are also able to generate a series of APs, known as a spike train, which plays a vital role in neuronal communication.^31–35 Inspired by Hodgkin–Huxley (HH) circuits,^36–43 developed by treating the neuronal membrane as a circuit,³⁶ some iontronic HH circuits were proposed that reproduce neuronal spiking features,^7,8 where the circuit in ref. 8 applies UP conical channels. Since the BP channels of interest in this manuscript show more pronounced memristive properties compared to UP channels, we expect to be able to improve upon the circuit described in ref. 8 by considering parameters that are experimentally more accessible and closer to their biological analogues.

In an attempt to reproduce the all-or-none APs and the spike train found in biological neurons and in the iontronic circuit in ref. 8, we consider the circuit architecture presented in ref. 8, shown in Fig. 4(a), where we replace the UP channels with BP channels with conductances g₊, g₋ and g_s and consider a new set of circuit parameters. To separate out the response times of these channels we set the channel lengths to be L_± = 1 μm and L_s = 15 μm. Through eqn (13) this translates to τ_± ≈ 0.042 ms for the two fast channels, while the slow channel has a timescale τ_s ≈ 9.4 ms ≫ τ_±. The batteries, with which the conical channels are connected in series, have potentials E_± = ±114 mV for the two fast channels and E_s = −180 mV for the slow channel. These batteries are the circuit analogues of the Nernst potentials due to concentration gradients over neuronal membranes, where we note that these battery potentials are within the range of typical mammalian Nernst potentials.¹¹ Moreover, the bulk concentration of ρ_b = 2 mM is close to typical mammalian extracellular K⁺ concentrations.¹¹ A capacitor is connected in parallel to the channels, with a capacitance C = 0.05 pF that again is close to typical biological values, as this corresponds to the capacitance of a mammalian neuronal membrane of area ∼2–5 μm²,^85–91 which is of similar dimensions to the surface area of the channels.

Fig. 4 (a) Schematic representation of the circuit proposed in ref. 8, but now with three bipolar rather than three unipolar channels, connected in series to individual batteries and in parallel to a capacitor. The electric potential difference V_m(t) over the capacitor can be driven by an imposed stimulus current I(t). (b) The membrane potential V_m(t) resulting from an imposed subcritical (red) and supercritical (blue) current pulse I(t) of duration 20 ms and strengths 17.5 pA and 17.6 pA, respectively, as determined using eqn (16), displaying an all-or-none action potential, as can be seen by the jump in spike amplitude around I_AP = 17.5 pA shown in the inset. (c) The membrane potential V_m(t) as a result of an imposed subcritical (red) and supercritical (blue) sustained current I(t) of strengths 18 pA and 18.1 pA, respectively, where a spike train emerges for I(t) > I_train = 18 pA. The magnitudes of the membrane potentials before and during the APs are similar to those observed in mammalian APs.³¹

The electric potential V_m(t) over the circuit shown in Fig. 4(a) is equivalent to the membrane potential over a neuronal membrane³⁶ and responds to the imposed stimulus current I(t), which acts as the control parameter and determines whether spiking occurs. The time-evolution of V_m(t) is provided by Kirchhoff's law

where the conductances g_i(V_i(t),t) each evolve according to eqn (14) with the corresponding τ_i. The voltage arguments V_i of g_i,∞(V_i) are given by V₋(t) = V_m(t) − E₋, V₊(t) = −V_m(t) + E₊ and V_s(t) = −V_m(t) + E_s, with the different signs of the potentials corresponding to the different orientations of the channels as depicted in Fig. 4(a). Eqn (13), (14) and (16) form a closed set, which we numerically solve with initial conditions V(0) = −0.1 V and g_i(0) = g_0,i.

In Fig. 4(b) and (c) we show that we reproduce the same neuronal behaviour as found in ref. 8, in the form of all-or-none action potentials (Fig. 4(b)) and a spike train (Fig. 4(c)). Excitingly, the membrane potentials before and during the APs range from ∼ −70 mV to ∼50 mV and are therefore of similar magnitude to those observed in mammalian APs.³¹ This, combined with the biologically more relevant battery potentials E_i and bulk concentration ρ_b compared to the circuit with UP channels from ref. 8, may prove to be crucial for the integration of such an iontronic circuit with biological systems in future applications.

VI. Conclusion and outlook

In summary, we have presented a theoretical approximation of the voltage-dependent steady-state current and the dynamic conductive properties of conical channels that are filled with an aqueous electrolyte and carry an inhomogeneous surface charge. Specifically, we focus on a channel with a positive surface charge at the base and middle, and a negative surface charge at the tip, thus forming a bipolar channel. This channel exhibits significantly improved current rectification when compared to unipolar conical channels with homogeneous surface charges and otherwise identical parameters. For negative and moderately positive static potentials V ≲ 0.2 V, our analytical approximation of salt concentration profiles and time-dependent currents are found to be in good agreement with finite-element calculations of the PNPS eqn (2)–(5), providing a solid foundation for further investigation of the use of these channels in (neuromorphic) iontronic circuits. While the steady-state and time-dependent current–voltage relations also show good agreement for large potentials, we do observe some qualitative deviations for V ≳ 0.2 V in the salt concentration profiles and electric field profiles compared to finite-element calculations. Additionally, for large static potentials V ≳ 0.4 V we observed a non-linearity in the relation of fluid volume flow and applied potential, where the bipolar channel acts as a fluidic diode. We hypothesize that this is due to the strong salt depletion that bipolar channels exhibit at large potentials, which implies that the small-Debye-length assumption λ_D ≪ R(x) that underlies eqn (8) becomes increasingly less accurate. Although the microscopic salt concentration profiles and electric field profiles are not accurately predicted by the present analytical model, the overall steady-state and time-dependent conductance is still in good agreement, indicating that our presented analytical approximation is an effective tool for the exploration of bipolar channels for iontronic circuits.

By extending the analytical methodology of ref. 8 and 55 to bipolar conical channels, we have demonstrated its generalizability and potential for predicting features such as current rectification in a wider range of geometries and surface charge distributions. Our derived equations suggest that the model we present here may be directly applicable to predicting current rectification in bipolar cylindrical channels, rather than solely conical geometries, which was previously experimentally demonstrated.^79–81 Furthermore, since our model allows for any general linear increase in surface charge along the longitudinal axis, our approach may also aid in identifying optimized surface charge values, distributions, and geometries for iontronic systems, beyond the parameter set on which we focus in this work. These findings point towards the generality, utility and potential of this analytical methodology in the field of iontronics.

In addition to the implications for optimizing and understanding individual channel properties, this work has also highlighted the potential of this analytical approximation method in the context of exploring iontronic circuits. By modeling a Hodgkin–Huxley circuit with bipolar channels we are able to present a system that relies on battery potentials and on salt concentrations comparable to their biological analogues, and that produces all-or-none action potentials and spike trains with voltage membranes that closely resemble the values observed in biological systems. This suggests that further research in this direction may prove beneficial in the development of advanced iontronic devices with improved performance.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work is part of the D-ITP consortium, a program of the Netherlands Organisation for Scientific Research (NWO) that is funded by the Dutch Ministry of Education, Culture and Science (OCW). T. M. K. performed the calculations; T. M. K. and W. Q. B. conceptualized the work; T. M. K. and W. Q. B. developed the theory under the supervision of C. S. and R. v. R.

References

1. S. H. Han, M.-A. Oh and T. D. Chung, Chem. Phys. Rev., 2022, 3, 031302.
2. C. Yang, K. Hu, D. Wang, Y. Zubi, S. T. Lee, P. Puthongkham, M. V. Mirkin and B. J. Venton, Anal. Chem., 2019, 91, 4618.
3. K. Hu, D. Wang, M. Zhou, J. H. Bae, Y. Yu, H. Xin and M. V. Mirkin, Anal. Chem., 2019, 91, 12935.
4. A. Noy and S. B. Darling, Science, 2023, 379, 143.
5. C. Li, T. Xiong, P. Yu, J. Fei and L. Mao, ACS Appl. Bio Mater., 2021, 4, 71.
6. B. Xie, T. Xiong, W. Li, T. Gao, J. Zong, Y. Liu and P. Yu, Chem.–Asian J., 2022, 17, e202200682.
7. P. Robin, N. Kavokine and L. Bocquet, Science, 2021, 373, 687.
8. T. M. Kamsma, W. Q. Boon, T. ter Rele, C. Spitoni and R. van Roij, Phys. Rev. Lett., 2023, 130(26), 268401.
9. P. Robin, T. Emmerich, A. Ismail, A. Niguès, Y. You, G.-H. Nam, A. Keerthi, A. Siria, A. Geim and B. Radha, et al. , Science, 2023, 379, 161.
10. T. Xiong, C. Li, X. He, B. Xie, J. Zong, Y. Jiang, W. Ma, F. Wu, J. Fei and P. Yu, et al. , Science, 2023, 379, 156.
11. L. Squire, D. Berg, F. Bloom, S. du Lac, A. Ghosh and N. Spitzer, Fundamental Neuroscience, Academic Press, 3rd edn, 2008, ch. 11.
12. M. P. Sah, H. Kim and L. O. Chua, IEEE Circuits Syst. Mag., 2014, 14, 12.
13. D. Wang, M. Kvetny, J. Liu, W. Brown, Y. Li and G. Wang, J. Am. Chem. Soc., 2012, 134, 3651.
14. Y. Li, D. Wang, M. M. Kvetny, W. Brown, J. Liu and G. Wang, Chem. Sci., 2015, 6, 588.
15. D. Wang, J. Liu, M. Kvetny, Y. Li, W. Brown and G. Wang, Chem. Sci., 2014, 5, 1827.
16. D. Wang and G. Wang, J. Electroanal. Chem., 2016, 779, 39.
17. D. Wang, W. Brown, Y. Li, M. Kvetny, J. Liu and G. Wang, Anal. Chem., 2017, 89, 11811.
18. Q. Sheng, Y. Xie, J. Li, X. Wang and J. Xue, Chem. Commun., 2017, 53, 6125.
19. W. Brown, Y. Li, R. Yang, D. Wang, M. Kvetny, H. Zheng and G. Wang, Chem. Sci., 2020, 11, 5950.
20. W. Brown, M. Kvetny, R. Yang and G. Wang, J. Phys. Chem. C, 2022, 126, 10872.
21. W. Brown, M. Kvetny, R. Yang and G. Wang, J. Phys. Chem. C, 2021, 125, 3269.
22. D. Wang, W. Brown, Y. Li, M. Kvetny, J. Liu and G. Wang, ChemElectroChem, 2018, 5, 3089.
23. P. Ramirez, J. J. Perez-Grau, J. Cervera, S. Nasir, M. Ali, W. Ensinger and S. Mafe, Appl. Phys. Lett., 2021, 118, 181903.
24. G. Sun, Z. Slouka and H.-C. Chang, Small, 2015, 11, 5206.
25. L. Chua, Nanotechnology, 2013, 24, 383001.
26. Y. van De Burgt, A. Melianas, S. T. Keene, G. Malliaras and A. Salleo, Nat. Electron., 2018, 1, 386.
27. S. T. Keene, P. Gkoupidenis, and Y. Van de Burgt, in Organic Flexible Electronics, Elsevier, 2021, pp. 531–574.
28. E. Chicca and G. Indiveri, Appl. Phys. Lett., 2020, 116, 120501.
29. D. V. Christensen, R. Dittmann, B. Linares-Barranco, A. Sebastian, M. Le Gallo, A. Redaelli, S. Slesazeck, T. Mikolajick, S. Spiga and S. Menzel, et al. , Neuromorphic Comput. Eng., 2022, 2, 022501.
30. K. Lucas, J. Physiol., 1909, 38, 113.
31. B. P. Bean, Nat. Rev. Neurosci., 2007, 8, 451.
32. L. Squire, D. Berg, F. Bloom, S. du Lac, A. Ghosh and N. Spitzer, Fundamental Neuroscience, Academic Press, 3rd edn, 2008, ch. 6.
33. G. S. Cymbalyuk, Q. Gaudry, M. A. Masino and R. L. Calabrese, J. Neurosci., 2002, 22, 10580.
34. E. Marder and V. Thirumalai, Neural Networks, 2002, 15, 479.
35. S. M. Sherman, Trends Neurosci., 2001, 24, 122.
36. A. L. Hodgkin and A. F. Huxley, J. Physiol., 1952, 117, 500.
37. W. Rall, Comprehensive Physiology, 2011, ch. 3, suppl. 1, pp. 39–97,  DOI:10.1002/cphy.cp010103.
38. R. FitzHugh, J. Theor. Biol., 1973, 40, 517.
39. W. Rall, Ann. N. Y. Acad. Sci., 1962, 96, 1071.
40. J. A. Halter and J. Clark Jr, J. Theor. Biol., 1991, 148, 345.
41. E. Hay, S. Hill, F. Schürmann, H. Markram and I. Segev, PLoS Comput. Biol., 2011, 7, e1002107.
42. M. L. Hines and N. T. Carnevale, Neural Comput., 1997, 9, 1179.
43. M. H. Kole, S. U. Ilschner, B. M. Kampa, S. R. Williams, P. C. Ruben and G. J. Stuart, Nat. Neurosci., 2008, 11, 178.
44. A. Mehonic and A. J. Kenyon, Nature, 2022, 604, 255.
45. V. K. Sangwan and M. C. Hersam, Nat. Nanotechnol., 2020, 15, 517.
46. S. T. Keene, C. Lubrano, S. Kazemzadeh, A. Melianas, Y. Tuchman, G. Polino, P. Scognamiglio, L. Cinà, A. Salleo and Y. van de Burgt, et al. , Nat. Mater., 2020, 19, 969.
47. P. C. Harikesh, C.-Y. Yang, D. Tu, J. Y. Gerasimov, A. M. Dar, A. Armada-Moreira, M. Massetti, R. Kroon, D. Bliman and R. Olsson, et al. , Nat. Commun., 2022, 13, 1.
48. I. Krauhausen, D. A. Koutsouras, A. Melianas, S. T. Keene, K. Lieberth, H. Ledanseur, R. Sheelamanthula, A. Giovannitti, F. Torricelli and I. Mcculloch, et al. , Sci. Adv., 2021, 7, eabl5068.
49. P. D. Marasco, J. S. Hebert, J. W. Sensinger, D. T. Beckler, Z. C. Thumser, A. W. Shehata, H. E. Williams and K. R. Wilson, Sci. Robot., 2021, 6, eabf3368.
50. L. Yuan, S. Liu, W. Chen, F. Fan and G. Liu, Adv. Electron. Mater., 2021, 7, 2100432.
51. C. D. Schuman, T. E. Potok, R. M. Patton, J. D. Birdwell, M. E. Dean, G. S. Rose, and J. S. Plank, arXiv, 2017, preprint, arXiv:1705.06963.
52. T. Venkatesan and S. Williams, Appl. Phys. Rev., 2022, 9, 010401.
53. J. Zhu, T. Zhang, Y. Yang and R. Huang, Appl. Phys. Rev., 2020, 7, 011312.
54. C. Wei, A. J. Bard and S. W. Feldberg, Anal. Chem., 1997, 69, 4627.
55. W. Q. Boon, T. E. Veenstra, M. Dijkstra and R. van Roij, Phys. Fluids, 2022, 34, 101701.
56. H. S. White and A. Bund, Langmuir, 2008, 24, 2212.
57. L. Jubin, A. Poggioli, A. Siria and L. Bocquet, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 4063.
58. I. Vlassiouk, T. R. Kozel and Z. S. Siwy, J. Am. Chem. Soc.T, 2009, 131, 8211.
59. L.-J. Cheng and L. J. Guo, Nano Lett., 2007, 7, 3165.
60. Z. S. Siwy, Adv. Funct. Mater., 2006, 16, 735.
61. S. N. Bush, T.T. Volta and C. R. Martin, Nanomaterials, 2020, 10, 571.
62. Z. Siwy, Y. Gu, H. Spohr, D. Baur, A. Wolf-Reber, R. Spohr, P. Apel and Y. Korchev, Europhys. Lett., 2002, 60, 349.
63. A. Fuliński, I. Kosińska and Z. Siwy, New J. Phys., 2005, 7, 132.
64. Z. Siwy, I. Kosińska, A. Fuliński and C. Martin, Phys. Rev. Lett., 2005, 94, 048102.
65. D. Duleba, P. Dutta, S. Denuga and R. P. Johnson, ACS Meas. Sci. Au, 2022, 2(3), 271–277.
66. W.-J. Lan, M. A. Edwards, L. Luo, R. T. Perera, X. Wu, C. R. Martin and H. S. White, Acc. Chem. Res., 2016, 49, 2605.
67. I. Vlassiouk, S. Smirnov and Z. Siwy, ACS Nano, 2008, 2, 1589.
68. J. Liu, M. Kvetny, J. Feng, D. Wang, B. Wu, W. Brown and G. Wang, Langmuir, 2012, 28, 1588.
69. C. Kubeil and A. Bund, J. Phys. Chem. C, 2011, 115, 7866.
70. S. Dal Cengio and I. Pagonabarraga, J. Chem. Phys., 2019, 151, 044707.
71. A. R. Poggioli, A. Siria and L. Bocquet, J. Phys. Chem. B, 2019, 123, 1171.
72. Y. Uematsu, Physics of Fluids, 2022.
73. X. Huang, X.-Y. Kong, L. Wen and L. Jiang, Adv. Funct. Mater., 2018, 28, 1801079.
74. I. Vlassiouk and Z. S. Siwy, Nano Lett., 2007, 7, 552.
75. M. L. Kovarik, K. Zhou and S. C. Jacobson, J. Phys. Chem. B, 2009, 113, 15960.
76. D. R. Lide, CRC Handbook of Chemistry and Physics, CRC press, 2004, ch. 5, vol. 85.
77. E. Choi, C. Wang, G. T. Chang and J. Park, Nano Lett., 2016, 16, 2189.
78. M. Shen, H. Yang, V. Sivagnanam and M. Gijs, Anal. Chem., 2010, 82, 9989.
79. R. Karnik, C. Duan, K. Castelino, H. Daiguji and A. Majumdar, Nano Lett., 2007, 7, 547.
80. H. Daiguji, Y. Oka and K. Shirono, Nano Lett., 2005, 5, 2274.
81. Z. Meng, Y. Chen, X. Li, Y. Xu and J. Zhai, ACS Appl. Mater. Interfaces, 2015, 7, 7709.
82. B. L. Werkhoven and R. van Roij, Soft Matter, 2020, 16, 1527.
83. V. S. Markin, A. G. Volkov and L. Chua, Plant Signaling Behav., 2014, 9, e972887.
84. L. Chua, Semicond. Sci. Technol., 2014, 29, 104001.
85. L. J. Gentet, G. J. Stuart and J. D. Clements, Biophys. J., 2000, 79, 314.
86. W. Niles, R. Levis and F. Cohen, Biophys. J., 1988, 53, 327.
87. C. Solsona, B. Innocenti and J. M. Fernández, Biophys. J., 1998, 74, 1061.
88. V. L. Sukhorukov, W. M. Arnold and U. Zimmermann, J. Membr. Biol., 1993, 132, 27.
89. G. Major, A. U. Larkman, P. Jonas, B. Sakmann and J. Jack, J. Neurosci., 1994, 14, 4613.
90. D. Thurbon, H.-R. Lüscher, T. Hofstetter and S. J. Redman, J. Neurophysiol., 1998, 79, 2485.
91. R. A. Chitwood, A. Hubbard and D. B. Jaffe, J. Physiol., 1999, 515, 743.

---