Unveiling the capabilities of bipolar conical channels in neuromorphic iontronics

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 [Kamsma et al., arXiv:2301.06158, 2023], 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 its potential for bio-compatibility.

In this work we present an analytical model that quantitatively captures both the steady-state and the dynamical behaviour of conical channels with an inhomogeneous surface charge distribution, based on the methodology in Refs.[8,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][32][33][34][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, how-ever with more biologically relevant salt concentrations and battery potentials.Furthermore, the circuit's membrane potentials during spiking match typical mammalian values, making it more bio-compatible.

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 the longitudinal coordinate, as depicted schematically in Fig. 1(a).The channel radius is described by 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, a negative surface charge at the tip, with a linear increase described by where throughout this manuscript we set σ = −3σ 0 /2 with eσ 0 = 0.1 enm −2 , 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 as realised before experimentally [75].The channel connects two bulk reservoirs of an incompressible aqueous 1:1 electrolyte, with mass density ρ m = 10 3 kg • m −3 , viscosity η = 1.01 mPa • s, and electric permittivity ε = 0.71 nF • m −1 .The electrolyte contains ions with monovalent charges ±e with e the proton charge and diffusion coefficients D ± = D = 2 µm 2 ms −1 , a typical value for dilute KCl [76] 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 ψ 0 (x) = (2k B T /e) sinh −1 [2πλ B λ D σ (x)] equals ψ 0 (0) ≈ 92 mV at the base and ψ 0 (L) ≈ −61 mV at the tip.Here we introduced the Bjerrum length λ B = e 2 /(4πεk B T ) and the Debye length λ D = 1/ 8πλ B ρ b .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 0 .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 − the charge flux and j + + j − the salt flux, and the pressure profile P(x, r,t).
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 V(t) (b) 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 σ (x) = σ 0 + σ x L .Here σ = −3σ 0 /2 with eσ 0 = 0.1 enm −2 such that the surface charge is positive at the base (eσ (0) = eσ 0 ) and negative at the tip (eσ (L) = −eσ 0 /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) 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 equations (2)-( 5), for a bipolar (BP) channel (blue) and otherwise identical unipolar (UP) channels with uniform surfaces charge −σ 0 /2 (green) and −σ 0 (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 [55].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), that satisfies the Poisson equation where Ψ(−∞, r,t) = V (t) and Ψ(∞, r,t) = 0. Eqs. ( 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 the normal vector of the wall.
In Fig. 1(b) we show steady-state current-voltage (I-V) curves as determined by finite-element (FE) calculations, not only for the BP channel under consideration (blue) but also for homogeneous and unipolar (UP) surface charge densities −eσ 0 /2 (green) and −eσ 0 (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 found in conical channels [55].It is clear that the BP channel exhibits a significantly stronger current rectification, with the ratio ICR = |I(−0.8V)/I(0.8V)| of the current I(V ) at voltages V = ±0.8V being as large as ICR ≈ 21 for the BP channel (blue), while it is as small as ICR ≈ 3 and 2.4 for the UP channels with surface charges −σ 0 /2 (green) and −σ 0 (red), respectively.

III. ANALYTIC APPROXIMATION FOR BIPOLAR CHANNELS
The full PNPS equations ( 2)-( 5) cannot be solved analytically for the system of interest here in steady-state.However, with a few reasonable assumptions we can simplify them to obtain some closed-form analytic descriptions [55].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) ≈ ρ s (x) and the electric potential Ψ(x, r) ≈ Ψ(x) are radially independent.With this assumption, as in Ref. [55], the slab-averaged electric field −∂ x Ψ(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 with Q(V ) = −πR t R b εψ eff V /(ηL) the electro-osmotic fluid volume flow, which is similar to the expression for the fluid flow of a UP channel [55] 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 Sec.IV.Eq. ( 7) represents the diffusive, conductive and convective components of the salt flux, respectively.In steady-state the condition ∂ x J(x) = 0 must hold, yielding for given σ (x), R(x), and Q(V ) a differential equation for the unknown radially averaged salt concentration profile function ρ s (x).In Ref. [55] this differential equation is solved for a conical channel with a homogeneous surface charge and boundary conditions ρ s (0) = ρ s (L) = 2ρ b .Here we consider the case where the surface charge distribution is given by Eq. (1).By solving ∂ x J(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 2 t ) the Péclet number at the narrow end.Note that for our case with solely a voltagedriven flow without any pressure-driven contribution, Q(V ) = −πR t R b εψ eff V /(ηL) is proportional to V , and hence the ratio Pe/V that appears in Eq. ( 8) does not depend on the static potential V .In Sec.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 Eq. ( 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 Eq. ( 8).However, we will discuss in Sec.IV how we can still ensure good agreement on the current-voltage relation over a wide voltage range.
Interestingly, Eq. ( 8) suggests that it can also explain and predict current rectification in cylindrical channels [79][80][81] as long as σ = 0, since in this case a non-trivial source term remains in Eq. ( 8) even for ∆R = 0. Hence our current work suggests to unify the theories for non-linear transport through cylindrical and conical channels carrying homogeneous or inhomogeneous surface charges.Additionally we note that Eq. ( 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 σ 0 ∆R = −σ R b (which for our linear surface charge density profile implies σ t /σ b = R t /R b with σ t and σ b 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 Sec.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 2 (x).Since the electric field scales with the inverse of R 2 (x) according to Eq. ( 6), the contribution to the resistance of each slab equals [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 ρ s (x,V ) by max [0.2ρ b , ρ s (x,V )] in the actual (numerical) evaluations of Eq. ( 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 steadystate current is then given by As we will show in Sec.IV, Eq. ( 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 diffusionlike 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 equations ( 2)-( 5), we apply the same methodology.We consider two quantities, the total number of ions N = π L 0 R 2 (x)ρ s (x)dx and the net salt flux J x (0) − J x (L) into the channel.The change of N given by Eq. ( 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 ρs (x) = 2ρ b .With this assumption the first and third terms in Eq. ( 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 it takes for ion depletion or accumulation, and thus the typical memory retention timescale is then approximated by τ = α/γ.This yields, perhaps surprisingly, the purely diffusive timescale identical to the expression for UP channels [8], which is remarkable as the conductive terms in Eq. ( 7), through which Eq. ( 13) is obtained, differ from those in Ref. [8].For our standard parameter set we find τ ≈ 4.17 ms.By assuming that ∂ t g(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) as and the current I(t) as Despite the fact that Eq. ( 9) needs to be evaluated numerically, we will refer to Eqs. ( 8), ( 10) and ( 15) as an analytic approximation (AA) for the voltage-dependent salt concentration profiles, steady-state current, and time-dependent current, respectively.In Sec.IV we will verify these three equations against full FE calculations of the PNPS equations ( 2)-( 5).

IV. FINITE-ELEMENT VERIFICATION
In Sec.III we derived an AA for the voltage-dependent salt concentration profiles, steady-state current, and timedependent current.Here we will verify these results against full FE calculations of the underlying PNPS equations ( 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 Eq. ( 8) (solid lines) with 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 a neighborhood 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 ρ s by max [0.2ρ b , ρ s (x,V )] in Eq. ( 9).
In Fig. 2(b) we translate the concentration profiles at static potentials V to the steady-state current-voltage relation I(V ) through Eqs. ( 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 of accuracy of the AA for this voltage regime.2)-( 5) and our analytic approximation (AA) of Eqs. ( 8), ( 10) and ( 15), all for our standard parameter set of a bipolar conical channel (see text).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 ρ s (x,V ) by max [0.2ρ b , ρ s (x,V )] in Eq. ( 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 Sec.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 ±1 V and frequency f = 45 Hz, is in line with the prediction that τ f ≈ 0.19 yields the most pronounced memory effect [8].We compare the time-dependent current determined through Eq. ( 15) (red) against FE calculations (blue).In both instances a similar pinched hysteresis loop is found, the hallmark of a memristor [84].We note that this hysteresis loop shows a much more pronounced opening compared to a loop of a similar UP channel [8], 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 Sec.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 ψ(x) and the fluid flow Q(V ) are accurately described by our AA for BP channels for various static V .In the AA −∂ x ψ(x) is given by Eq. ( 6), which shows good agreement for UP conical channels in the present parameter regime [55].In Fig. 3(a) we compare the electric field for various static potentials V as predicted by Eq. ( 6) (solid lines) with FE calculations (circles).For negative and moderately positive potentials we find good agreement, as expected on 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 divergentfree) 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 ψ(x) for V = 300 mV in Fig. 3(a) is not observed in UP channels in a similar parameter regime [55].
The underlying Eq. ( 8) of the reported result is dependent on the fluid volume flow Q(V ) = −πR t R b εψ eff V /(η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 ) ∝ ψ 0 with ψ 0 the (homogeneous) surface potential [55], but in BP channels such a relation is not obvious as the surface potential ψ 0 (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, also the strength of ψ eff = −25 mV seems to be similar to the average surface potential of the tip L 2/3L ψ 0 (x)dx/(L/3) ≈ −34 mV.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 Sec.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 Figs.2(b) and 2(c), which are most relevant in the context of iontronic circuitry.Furthermore, the iontronic circuit presented in Sec.V operates within a voltage regime where the electric fields and fluid flows predicted by the analytic approximation are reasonably consistent with 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 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][32][33][34][35]. Inspired by Hodgkin-Huxley (HH) circuits [36][37][38][39][40][41][42][43], developed by treating the neuronal membrane as a circuit [36], 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 Eq. ( 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 [11].Moreover, the bulk concentration of ρ b = 2 mM is close to typical mammalian extracellular K + concentrations [11].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 mammalian neuronal membrane of area ∼ 2 − 5 µm 2 [85][86][87][88][89][90][91], which is of similar dimensions as the surface area of the channels.
The electric potential V m (t) over the circuit shown in Fig. 4(a) is equivalent to the membrane potential over a neuronal membrane [36] 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 Eq. ( 14) with the corresponding τ i .The voltage arguments with the different signs of the potentials corresponding to the different orientations of the channels as depicted in Fig. 4(a).Eqs. ( 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 .
4. (a) Schematic representation of the circuit proposed in Ref. [8], however 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 by Eq. ( 16), displaying an all-or-none action potential, as can be seen by the jump in spike amplitude around I AP = 17.5 pA as shown in the inset.(c) The membrane potential V m (t) as a result of an imposed subcritical (red) and supercritical (blue) sustained currents I(t) of strengths 18 pA and 18.1 pA respectively, where a spike train emerges for I(t) > I train = 18 pA.The magnitude of the membrane potentials before and during the APs are similar to those observed in mammalian APs [31].
In Figs.4(b) and 4(c) we show that we reproduce the same neuronal behaviour as found in Ref. [8], in the form of all-ornone 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 [31].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 analytic approximation of salt concentration profiles and time-dependent currents are found to be in good agreement with finite-element calculations of the PNPS equations ( 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-Debyelength assumption λ D R(x) that underlies Eq. ( 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 steadystate 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 Refs.[8,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 is previously experimentally demonstrated [79][80][81].Furthermore, since our model allows for any 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 analytic 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 which 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.

FIG. 2 .
FIG.2.Comparisons of finite-element calculations (FE) of the full PNPS equations (2)-(5) and our analytic approximation (AA) of Eqs.(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 by Eq. (8) (solid lines) and by 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 Eq. (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.
FIG.2.Comparisons of finite-element calculations (FE) of the full PNPS equations (2)-(5) and our analytic approximation (AA) of Eqs.(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 by Eq. (8) (solid lines) and by 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 Eq. (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.

FIG. 3 .
FIG. 3. (a) Steady-state electric field −∂ x Ψ(x) inside the channel as predicted by Eq. (6) (solid lines) and as measured on the central axis of the channel through the FE calculations (circles) of the PNPS equations (2)-(5) for various static applied potentials V .(b) Steadystate fluid volume flow Q(V ) as a function of the static potential V as predicted by Q(V ) = −πR t R b εψ eff V /(η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).