Advanced iontronic spiking modes with multiscale diffusive dynamics in a fluidic circuit

Fluidic iontronics is emerging as a distinctive platform for implementing neuromorphic circuits, characterized by its reliance on the same aqueous medium and ionic signal carriers as the brain. Drawing upon recent theoretical advancements in both iontronic spiking circuits and in dynamic conductance of conical ion channels, which form fluidic memristors, we expand the repertoire of proposed neuronal spiking dynamics in iontronic circuits. Through a modelled circuit containing channels that carry a bipolar surface charge, we extract phasic bursting, mixed-mode spiking, tonic bursting, and threshold variability, all with spike voltages and frequencies within the typical range for mammalian neurons. These features are possible due to the strong dependence of the typical conductance memory retention time on the channel length, enabling timescales varying from individual spikes to bursts of multiple spikes within a single circuit. These advanced forms of neuronal-like spiking support the exploration of aqueous iontronics as an interesting platform for neuromorphic circuits.


Introduction
In the pursuit of brain-inspired circuits the focus is often on the synaptic properties of neuromorphic devices, where synapses are considered as primary computational units in neuromorphic computing [1].Consequently, due to their analogous behaviour to synapses, memristors have significantly shaped and driven research in this domain, where the timeand history-dependent conductance of memristors offers a versatile platform for emulating features of synaptic plasticity [2][3][4].However, synapses are not the only components in the brain which can be emulated with memristors.The biological ion channels responsible for generating action potentials also exhibit memristive behavior [5].This is underscored by the seminal Hodgkin-Huxley (HH) model [6], which mathematically describes the axonal membrane potential by treating the membrane as an equivalent electric circuit in which the ion channels embedded in the axonal membrane are modelled as circuit components.The mathematical models for these ion channels were later recognised as descriptions of memristors [7].Although both synapses and axonal ion channels are neuronal components that can be described and emulated by memristors, they are explicitly distinct biological structures which carry out different tasks.This biological nuance sometimes leads to confusion and inaccurate descriptions of memristive devices in the brain, such as incorrectly associating the HH model with descriptions of synapses [8].Nevertheless, the intriguing connection between memristors and the HH model has also sparked considerable interest [5,9] and neuronal signalling has inspired various circuits that capture various features of neuronal spiking [10,11].
Biological neurons feature a wealth of different spiking modes, which can be clearly categorised and used to judge the quality of neuron models [12].Typically the most basic features to consider are tonic spiking, a regular train of voltage spikes with constant frequency, and phasic spiking, a single isolated voltage spike.In the case of phasic spiking, the neuron model should also obey the all-or-none law [4,13], i.e. a voltage spike is either fully generated upon a sufficiently strong impulse, or the voltage fails to spike, with no intermediate transition in between.However, many more neuronal firing modes are recognised and this signalling behaviour of neurons has inspired various circuits that can emulate a wide array of different modes of neuronal spiking [10,11].Examples, that will also feature in the present study, include phasic bursting, mixed mode spiking, tonic bursting (otherwise known as chattering [14]), and threshold variability [12].In phasic bursting, a single burst of several spikes emerges upon applying a sustained stimulus, after which the system again settles to a steady state, despite the constant and sustained current stimulus.Mixed mode spiking consists of an initial burst of spikes upon a sustained stimulus, followed by tonic spiking.In tonic bursting, short periods of spiking, i.e. bursts, are interchanged by short periods of no spiking at all.Lastly, threshold variability indicates that the threshold for a neuron to spike can depend on the prior activity of the neuron.
The vast majority of neuromorphic devices (both spiking and synaptic) consist (at least partially) of solid-state components [2,3,10,11], which results in fundamental differences with biological neurons.For instance, while solid-state devices typically rely on a single information carrier, such as electrons or holes, driven only by electric forces, neurons employ the transport of various ions and molecules in parallel, while combining electrical and chemical regulation, both for signalling [16] and for synaptic transmission [17][18][19].Additionally, the fast dynamics of solid state components can be a disadvantage when temporal inputs are natural or biological signals as the typical timescales of those inputs can be significantly slower than those of solid-state devices, therefore requiring complicated virtual clocks for synchronisation [20,21].Recent work tries to address and overcome these limitations through electrochemical coupling of solid-state com- eσ(x) g ss

E ss
FIG. 1.(a) Schematic representation of the proposed fluidic iontronic circuit featuring four channels of three different types.Two short channels of equal length L ± = 1 µm with fast dynamics on a typical timescale τ ± ≈ 0.042 ms and conductances g ± (t), a longer channel of length L s = 15 µm with slower dynamics on a typical timescale τ s ≈ 9.4 ms and conductance g s (t), and an even longer channel of length L ss = 90 µm with conductance g ss (t) and the slowest dynamics over a typical timescale τ ss ≈ 338 ms.These channels are connected in series with batteries with potential E ± = ±114 mV and E s = E ss = −180 mV, respectively, and in parallel to a capacitor of capacitance C = 0.05 pF.
A time-dependent stimulus current I(t) can be imposed through the circuit and a potential V m (t) forms over the circuit that is equivalent to the neuronal membrane potential [6].Schematic adapted from Ref. [15].(b) Schematic of an individual bipolar channel of length L i , with base radius R i,b and tip radius R i,t , connecting two aqueous 1:1 electrolyte reservoirs of concentration ρ b = 2 mM.The wall of all four channels carries an inhomogeneous surface charge that linearly decreases from 0.1 enm −2 at the base to −0.05 enm −2 at the tip.
ponents to ionic systems, both in the context of synaptic devices [22,23] and for spiking circuits [24][25][26].However, a newly emerging direction proposes to omit solid-state components altogether, and hence the need for any chemical or ionic coupling, by implementing neuromorphic features in an aqueous electrolyte medium [27][28][29][30][31][32][33][34].These (fluidic) iontronic devices have recently garnered significant interest, offering the promise of multiple information carriers, chemical regulation, and bio-integrability [35], although sacrificing on the high speeds obtainable by solid state devices.Unlike traditional solid-state neuromorphic circuits, fluidic iontronic circuits leverage the dynamic interplay of ions within an aqueous electrolyte, mirroring the conductive and fluidic characteristics inherent in biological neuronal environments.This departure from solid-state components introduces a novel dimension to neuromorphic computing, offering the potential for closer emulation of the brain's aqueous dynamics [36,37].Recent advances include chemical regulation [30,31] and initial demonstrations of iontronic neuromorphic computing [38].However, the development of neuromorphic iontronic devices is still in its infancy, requiring further theoretical explorations and experimental investigations to establish their capabilities in emulating complex neuronal functionalities [28,34,35].
In the recent rise of interest in iontronic neuromorphics, spiking circuits also received some attention in the form of theoretical studies, where HH-inspired iontronic circuits are modelled and shown to exhibit features of neuronal spiking [15,29].These proposals feature a circuit composed of an aqueous electrolyte medium, akin to the neuronal medium that the HH model describes, and rely on fluidic iontronic memristors to induce neuronal spiking.Initially, tonic spiking was shown to emerge from a circuit containing angstromscale slits [29], shortly after which an alternative iontronic circuit exploiting conical ion channels was proposed that ex-hibits both the characteristic all-or-none phasic spiking and tonic spiking [15].Thus, the two modes that are typically considered first [11,12] have been theoretically predicted to also emerge from fluidic iontronic circuits.However, no proposals yet exist to also include other spiking modes.
In this work we expand upon the previously reported features of neuronal spiking in fluidic iontronics [15,29].By building upon a previously reported iontronic circuit [15] and a physical description of the dynamical conductance of conical channels with a bipolar surface charge [39], i.e. positive at the base and negative at the tip, we can unlock various new forms of spiking dynamics.Due to the strong dependence of the typical conductance memory retention time on the channel length, we can implement timescales varying from individual spikes to bursts of multiple spikes within a single circuit, thereby enabling new spiking modes.Specifically these spiking modes are the aforementioned phasic bursting, mixed mode spiking, tonic bursting, and threshold variability [12].

Iontronic circuit and bipolar channels
Conical fluidic ion channels act as iontronic volatile memristors [40] and are being investigated as possible candidates for synaptic devices [41] and spiking circuits [15,39].Using theoretical models that quantitatively explain the memristive behaviour of conical channels, it was shown that HH-inspired fluidic circuits containing three conical channels and a capacitor exhibit tonic and phasic spiking [15,39].This modelled circuit was originally composed of conical ion channels with a homogeneous unipolar (UP) surface charge [15] and was later modified by replacing the UP channels with conical channels carrying a bipolar (BP) inhomogeneous surface charge [39], positive at the base and negative at the tip.BP fluidic chan-nels and (Janus) membranes have long drawn great interest as current rectifiers [42][43][44][45][46][47] for applications in e.g.sensing [47] and osmotic energy conversion [47][48][49].BP channels also show potential for iontronic memristors as a modification from a UP to a BP surface charge led, for an individual conical channel, to a much more pronounced current-voltage hysteresis loop upon applying an AC voltage, i.e. a stronger conductance memory effect.
Here we consider a circuit containing several of these conical BP channels, with different lengths L i .An important feature of these BP channel memristors is that their typical conductance memory timescale is dictated by the channel lengths L i according to with D = 2 µm 2 ms −1 the diffusion coefficient of the ions [39], which we assume to be identical for all ionic species for convenience.As we will discuss in Sec. 3 .6,the combination of channels of various lengths in a single circuit gives rise to dynamics on the timescale of individual spikes and of bursts of spikes.
To unlock additional features of neuronal firing, beyond tonic and phasic spiking, we introduce the circuit schematically depicted in Fig. 1(a), containing a capacitor with capacitance C = 0.05 pF, a typical capacitance for a mammalian neuronal membrane with an area of order ∼ 0.1 µm 2 [50,51], i.e. of the same order as the cross-sectional area of a channel.This capacitor is connected in parallel with four BP conical channels with conductances g + (t), g − (t), g s (t), and g ss (t), and four batteries each in series with the conical channels.The channels are taken to be of varying lengths L ± = 1 µm, L s = 15 µm, and L ss = 90 µm.Through Eq. (1) this translates to timescales τ ± ≈ 0.042 ms for the two fast channels, τ s ≈ 9.4 ms for the slow channel, and τ ss ≈ 338 ms for the super slow channel.The batteries have potentials E ± = ±114 mV for the two fast channels, and E s = E ss = −180 mV for the slow and super slow channels.These batteries, which mimic the Nernst potential caused by ionic concentration differences inside and outside the neuron in the HH model [6], are considered to be actual batteries in the microfluidic circuit of interest here, but their potentials are comparable to their biological Nernst potential counterparts [52].
In Fig. 1(b) we show a schematic depiction of a BP channel of length L i , implemented in the circuit in Fig. 1(a), with base-and tip radii R i,b and R i,t = R i,b − ∆R i , respectively, and thus with radius in the channel.The channel connects two 1:1 aqueous electrolyte reservoirs with the viscosity η = 1.01 mPa • s and the electric permittivity ε = 0.71 nF • m −1 of water.The cationic and anionic bulk concentrations are given by ρ b = 2 mM, comparable to the extracellular potassium concentration in biological neurons [52], which gives rise to a Debye length λ D ≈ 6.8 nm.The channels carry a surface charge that linearly decreases from eσ 0 = 0.1 enm −2 at the broad base to −0.05 enm −2 at the narrow tip, thereby changing by σ ′ = −3σ 0 /2 over the channel length and forming a bipolar surface charge profile.These charge densities correspond to Gouy-Chapman zeta potentials that vary between 92 mV and −61 mV.For the short fast channels and the slow channel we fix R i,b = 200 nm and R i,t = 50 nm, while the super slow channel is narrower with R ss,b = 120 nm and R ss,t = 30 nm.Thus, in all cases the channel radii are substantially larger than the Debye length, such that overlap of electric double layers is not prominent.
To fully resolve the dynamics of the circuit depicted in Fig. 1(a), we have to know how the conductances g i (t) of the BP channels evolve.For this we use an analytical model that quantitatively describes the steady-state and dynamical conductance properties of BP channels [39].BP channels exhibit voltage-dependent salt concentration polarisation in steadystate, with the radially averaged salt concentration ρ i,s (x,V i ) described by ) the Péclet number at the narrow end and The system is considered to be at a temperature of 293.15K and the effective surface potential ψ eff = −25 mV is taken to be the same as in Ref. [39] as we consider the same surface charge distributions here.The accumulation or depletion of salt affects the conductance of the channel according to with ) the homogeneous channel conductance.In the numerical evaluation of Eq. ( 3) we replace ρ i,s (x,V i ) by Max 0.2ρ b , ρ i,s (x,V i ) to avoid nonphysical negative concentrations that can emerge due to the strong voltage-dependent salt depletion of BP channels [39].This approach does induce a sharper drop in conductance, compared to full finite-element simulations, when concentrations start to approach the imposed minimum of 0.2ρ b , discussed in more detail in the Supplemental Material.This artefact complicates the circuit equations we introduce below.To help smooth over this sharper drop we employ a third-order interpolation to evaluate Eq. ( 3) between voltages spaced at intervals of 0.025 V, ranging from -0.3125 to 0.3125 V.A more sophisticated theoretical model of individual channels in the future should obviate the need for such an ad hoc approach, but for now this effective method suffices.Since it takes a typical time τ i as per Eq. ( 1) for salt to accumulate or deplete, the channel exhibits a (volatile) memory conductance with typical memory retention time τ i .The resulting dynamic conductance g i (t) was found to be well described by Various modes of voltage spiking (blue curves) extracted by modeling one and the same iontronic circuit driven by different timedependent currents (red curves), with (a) tonic spiking [15,29] for I = 19.05pA and (b) phasic spiking [15] for I = 18.4 pA reported before in iontronic circuits.The newly introduced modes of iontronic spiking dynamics include (c) phasic bursting for I = 19.01pA, i.e. a burst of spikes followed by a return to a steady state upon a sustained stimulus, (d) mixed-mode spiking for I = 19.02pA, i.e. an initial high-frequency burst of spikes followed by a transition into lower frequency tonic spiking, (e) tonic bursting for I = 19.04pA, i.e. a short burst of spiking alternating with periods of quiescence, and (f) threshold variability, with variations in the firing threshold influenced by prior activity.The negative and positive current stimuli are of the same magnitude for I = 18.3 pA but with different time intervals between the negative and the subsequent positive pulse.The firing threshold is temporarily lowered by the negative pulse and therefore the positive pulse only surpasses the (variable) firing threshold when the time between the current pulses is sufficiently short.
where V i (t) is the potential difference between base and tip of the channel, g i,∞ (V i ) is the voltage-dependent steady-state conductance of the channel as per Eq. ( 3), and τ i is the typical conductance memory retention timescale of the channel given by Eq. ( 1) [39].With differential equations for each of the dynamic conductances g i (t), we only need one additional equation to close the set that describes the time-evolution of the "membrane" potential V m (t), here the potential over the capacitor.This additional equation is provided by Kirchhoff's law where i ∈ {+, −, s, ss} and the conductances g i (t) each evolve according to Eq. (4) with their corresponding g i,∞ (V i (t)) and τ i .The voltage arguments V i (t) over the channels are given by with the different signs of the potentials corresponding to the different orientations of the channels as depicted in Fig. 1(a).Using the initial conditions V (0) = −70 mV and g i (0) = g i,0 , with g i,0 as defined below Eq. ( 3), we numerically solve the closed set of Eqs. ( 1), ( 4) and ( 5) for various current stimuli I(t).The system is given at least 10 s to settle into a steady state before applying a current I(t), we offset the time in the results to omit this in the plots.Note that τ ± ≪ τ s , τ ss , additionally τ ± is much faster than the typical response time of V m (t) too as we will see later, so the two short channels actually act as quasi-instantaneous current rectifiers due to their comparatively fast dynamics, rather than memristors, as we will more extensively discuss in Sec. 3 .6.

Advanced iontronic spiking modes
Upon numerically evaluating the membrane potential V m (t) that emerges from the proposed fluidic iontronic circuit introduced in Sec. 2 for various stimuli, we reveal the remarkable diversity of typical neuronal firing modes [12] shown in Fig. 2,  (a which we will discuss individually below.We stress that all spiking modes discussed below originate from one and the same iontronic circuit, with the stimulus current I(t) the only difference between the spiking modes.Additionally, we note that all spikes exhibit voltage amplitudes and spiking frequencies that are typical for mammalian neurons [13].

Tonic Spiking and Phasic Spiking
The earlier reported foundational tonic [15,29] and phasic spiking [15] also emerge from the circuit we consider here.Tonic spiking, characterized by a regular train of voltage spikes as shown in Fig. 2(a), and phasic spiking, featuring a single isolated voltage spike as shown in Fig. 2(b), appear for the present system parameters under sustained current stimuli of 18.40 pA and 19.05 pA, respectively.The phasic spiking current stimulus of 18.40 pA is just above the threshold for any spiking to occur, unless we consider the variability of the threshold as discussed in Sec. 3 .5.The sustained current of Fig. 2(b) does give rise, after the single voltage pulse, to a steady voltage that differs from the initial voltage.An all-ornone spike can also appear upon a pulse stimulus, after which the voltage settles back to its initial steady-state [39].
The dynamics here are governed by the typical RC-like time of the circuit that determines the time it takes for the (de)polarisation of V m (t), while the timescale τ s dictates the typical width of a spike; the short channels respond on such fast timescales that their dynamics can be assumed to be instantaneous [15,53].Although the timescale τ ss does not play a role in these spiking modes as these also appear without the super slow channel [39], a small influence of the super slow channel is still visible in the case of tonic spiking.The spiking frequency initially is slightly higher immediately after the stimulus is applied and then gradually settles into a lower fre-quency over a time ∼ τ ss .This actually corresponds to the spiking mode of spike frequency adaptation [12], but since this effect is so minor in our results, we choose not to explicitly distinguish it as an additional emerging spiking mode.

Phasic Bursting
Imposing a sustained current stimulus to the circuit of 19.01 pA elicits phasic bursting, a spiking mode where a burst of spikes occurs, followed by a return to a (new) steady state, despite the sustained stimulus.This mode is made possible by the super slow channel.The initial burst has a duration of the typical timescale ∼ τ ss of the super slow channel, after which this channel has had sufficient time to increase its conductance to return the system to a steady state.

Mixed Mode
Under a sustained stimulus of 19.02 pA we find mixed mode spiking, i.e. the iontronic circuit transitions from an initial high-frequency burst of spikes with a duration of ∼ τ ss , into a lower frequency tonic spiking, with the individual spikes now separated by ∼ τ ss , as shown in Fig. 2(d).In this case the initial burst is a transient, of typical time ∼ τ ss , as the system settles into the periodic solution of the tonic spiking.

Tonic Bursting
Tonic bursting entails short bursts of spiking interspersed with periods of quiescence.When imposing a sustained stimulus of 19.04 pA we find that the circuit exhibits a periodic behaviour of high frequency burst as shown in Fig. 2(e).The durations of the bursts and periods of quiescence are dictated by the slow dynamics of longest channel, as the super slow channel periodically increases and decreases in conductance, visible by the fact that each burst or quiescence period has a duration of order ∼ τ ss ≈ 338 ms.

Threshold Variability
Our findings also unveil threshold variability, wherein the firing threshold of the neuron is influenced by prior activity.As shown in Fig. 2(f), when imposing a negative and positive stimulus pulse of magnitude ±18.30pA (just below the threshold mentioned in Sec. 3 .1) of duration 0.02 s, separated by 0.18 s (between the end of the first pulse and the beginning of the second) no spike occurs.However, when we impose precisely the same pulses but now separated by 0.01 s we find that a full spike occurs.Thus in the first set of pulses, the threshold for spiking was not reached, but in the second instance it was reached with exactly the same pulses, showing that the prior activity of the circuit can influence the threshold for spiking.This is a result of the slow channel with timescale τ s decreasing in conductance as a result of the negative pulse, while the super slow channel actually plays no role in this spiking mode as it is also observed without the super slow channel.If the interval is much larger than τ s ≈ 9.4 ms, as it is for the first set of pulses, then the slow channel reverts to its steady-state before the second pulse.However, if the interval between the stimuli is of the order of τ s = 9.4 ms (or smaller), as is the case in the second set of pulses where the interval is 10 ms, then the slow channel still has a lowered conductance when the second pulse arrives, making the system more susceptible to stimuli and thereby lowering the firing threshold.

Roles of the different channels
To elucidate the circuit design as shown in Fig. 1(b) we heuristically describe the roles the various channels play.Firstly, since τ ± is much shorter than the typical response time of V m (t) and since τ ± ≪ τ s , τ ss , the two short channels actually act as quasi-instantaneous current rectifiers, rather than memristors, though we solve the dynamic equation for all channels for completeness.Tonic and phasic spiking, which occur without the super slow channel, were already remarked to also emerge by using the instantaneous conductance g ∞,± (V (t)) [15] in Eq. ( 5), reducing such a three channel circuit to a twodimensional dynamic system with dynamic variables V m (t) and g s (t) [53].By extension the results we present here are represented by the three dynamic variables V m (t), g s (t), and g ss (t).In Fig. 3 we show these dynamic variables during the tonic bursting shown in Fig. 2(e).
All channels drive V m (t) toward their respective battery potentials.When V m (t) is near E ± the respective corresponding short fast channel has a high conductance and maintains V m (t) in that state, forming temporary stable states between which V m (t) switches during spiking.The slow channel, which drives V m (t) towards E s = −180 mV, is in a low conductance state when V m (t) is negative, allowing V m (t) to depolarize and switch to the positive voltage state upon a stimulus.Following the increase in V m (t), the slow channel increases in conductance over a timescale τ s as we show in Fig. 3(a), resulting in a consequent downward shift of V m (t).This behavior analogously resembles the delayed activation of K + channels in the Hodgkin-Huxley model [6].The super slow channel, operating over a timescale τ ss , plays a role akin to the slow channel.However, as we show in Fig. 3(b), it takes several spikes for the super slow channel to sufficiently increase in conductance and drive V m (t) towards its negative battery potential E ss = −180 mV, consequently suppressing spiking.The resulting period of quiescence lasts a time ∼ τ ss , forming a bursting process.The salt accumulation and depletion underpinning the conductance change of the super slow channel bear similarity to the slow intracellular Ca 2+ accumulation and depletion implicated in regulating bursting in biological neurons [54][55][56].Combining the three relevant dynamic variables V m (t), g s (t), and g ss (t) in one phase portrait yields Fig. 3(c), revealing a characteristic bursting trajectory with three loops (blue) and a path connecting the third and first loop (red), corresponding to the three spikes and to the periods of quiescence, respectively.

Discussion and conclusion
Previously reported fluidic iontronic circuits have demonstrated tonic spiking [15,29] and phasic spiking [15].In this study, we extend the repertoire of emergent spiking modes by introducing a new HH-like fluidic iontronic circuit, consisting of a capacitor and four iontronic memristors, that exhibits phasic bursting, mixed-mode spiking, tonic bursting, and threshold variability [12], as well as the earlier reported tonic and phasic spiking [15,29].The spikes in our proposed modes exhibit voltages and frequencies that align with those observed in mammalian neurons [13].Moreover, the capacitance, battery potentials and salt concentration in the circuit are comparable to their biological counterparts [52].Our theoretical framework builds upon a previously proposed iontronic circuit that exhibits tonic and phasic spiking [15] and a physical model for conical ion channels with a bipolar (rather than unipolar) surface charge [39].These channels are memristive [39] and their typical conductance memory retention time is dependent on the channel length.By varying the lengths of the four channels we can incorporate timescales on the order of a single spike and of entire bursts in a single circuit, allowing for the spiking and bursting processes that emerge from one and the same circuit.
While our theoretical framework in principle is fully physical, a limitation is the parameter sensitivity of the system, at least for the system parameters we considered.The stimuli strengths that induce different spiking modes are only separated by ∼ 0.01−0.1 pA on the scale of about 20 pA.Notably, if wider current stimuli intervals are found for spiking, then it is possible that class 2 spiking [12] can also be distinguished as a separate feature as the transitions in frequency seem to be discontinuous, but class 1 or 2 spiking is typically evaluated over varying stimuli intervals which in our case are too narrow to meaningfully investigate this.Additionally, although spiking was found to emerge for a wide range of different parameter configurations, the spiking is sensitive to small individual changes of the short fast channels or their respective batteries, where (at least one mode of) spiking only appeared in the tight interval E ± ∈ ± [113, 114.8] mV.The short fast channels play no dynamic roles in the circuit, but rather act as instantaneous current rectifiers that create the stable voltage states between which V m oscillates during spiking due to the dynamic switching of the (super) slow channel.Hence, no memristive properties are required for the short fast channels, which only offer a current rectification of around ≈ 21 [39], and other (perhaps better performing) diodic devices from the wide range of iontronic current rectifiers [57][58][59][60][61] could be considered.Although some devices can be described by similar theoretical models as we use here [38,61], our analysis is limited to devices for which we have an analytical quantitative model, but iontronic devices that lack such models are still feasible for experimental fabrication.For the (super) slow channels we do require the specific (length-dependent) volatile dynamics we find for the iontronic memristors of concern here, but in this case the system is far more stable against parameter shifts, thereby supporting the use of the (super) slow channels as described here.At least one mode of spiking emerges for changing a single parameter at a time in the range E s ∈ [−200, −90] mV, E ss ∈ [−450, −140] mV, while the capacitance can even span orders of magnitude C ∈ [10 −6 , 10 −1 ] pF.Therefore, the three components that govern the circuit dynamics, i.e. the (super) slow channels and the capacitor, are relatively robust once the short fast channels are in order.
The above suggestion for future improvements using other fluidic devices is supported by the fact that the results presented here are already an expansion on results we derived earlier for simpler unipolar conical channels carrying a homogeneous surface charge.Tonic bursting also emerges from a similar circuit with unipolar channels, but with circuit parameters (i.e. higher battery potentials, lower salt concentration, lower capacitance) and spiking voltages that are further removed from their biological analogs.The results and specific parameters for the unipolar channel circuit are laid out in the Supplemental Material.The emergence of tonic bursting in a different circuit with different fluidic memristors shows that the bursting spiking modes we present are not inherently dependent on the bipolar conical channels we consider here.Therefore, possible further improvements can be achieved by considering fluidic iontronic devices with an even wider range of attainable conductances.However, this is an issue of individual device physics and here we mostly focused on the overall circuit architecture and the spiking modes it enables.
In summary, we have considerably expanded the range of spiking modes proposed to emerge from iontronic fluidic circuits, entailing phasic bursting, mixed-mode spiking, tonic bursting, and threshold variability.The alignment of the spikes in our results with typical mammalian neuronal voltages and frequencies, combined with various circuit parameters that are comparable to their biological counterparts, further supports the potential that fluidic iontronics carry for neuromorphic spiking circuits.Moreover, since these biologically realistic spikes emerge from a circuit that is based upon the same aqueous electrolyte medium as in neurons, a unique perspective is the future possible integration with biological systems.However, the present system is rather sensitive to stimulus strengths and other circuit parameters, especially in the short fast channels, a limitation that may be mitigated by implementing fluidic devices with a broader range of available conductances.Nevertheless, we showed that the multiscale diffusive timescales of fluidic iontronic memristors of different lengths facilitate a relatively simple circuit that exhibits various advanced modes of neuronal spiking.Consequently, this work contributes to the ongoing exploration of fluidic iontronics as a promising platform for neuromorphic circuits, providing theoretical insights and proposed applications, thereby paving the way for future advancements in this burgeoning field.
Supplemental Material for: Advanced iontronic spiking modes with multiscale diffusive dynamics in a fluidic circuit In Fig. S1 we compare the steady-state conductance predicted by our analytical approximation (AA, red) of Eq. ( 3) from the main text with finite-element (FE, blue) calculations of the full Poisson-Nernst-Planck-Stokes equations.The precise details of these FE calculations are detailed in Ref. [1].When compared to finite-element calculations we see that the analytical approximation of Eq. ( 3), directly derived from the underlying Poisson-Nernst-Planck-Stokes equations, yields reasonable agreement, however it underestimates the increase in conductance for negative voltages while it overestimates the decrease in conductance for positive voltages.Eq. ( 3) also predicts a relatively sharp conductance drop when concentrations start to approach the imposed minimum salt concentration of 0.2ρ b at around 0.15 V in Fig. S1, which complicates the numerical solutions of the Kirchhoff equations of the main text.With a third-order interpolation of Eq. ( 3) this sharp drop can be somewhat smoothed over by choosing intervals ≳ 0.025 V, in our case in the regime from -0.3125 to 0.3125 V.In the future, a more sophisticated theoretical model of individual channels should address the physics that underlies this detail, which involves a surface contribution to the channel conductance.

Unipolar surface charge conical channel
The tonic bursting described in the main text was also found to emerge from an iontronic circuit containing unipolar (UP) conical channels with a homogeneous surface charge, rather than bipolar (BP) surface charge channels as considered in the main text.These UP channels were found to be memristors that facilitated tonic and phasic spiking upon coupling them in an iontronic circuit containing three UP channels and a capacitor [2].Here we use essentially the same parameters and theoretical model used for tonic and phasic spiking in Ref. [2] as a basis to also find tonic bursting.That is, all UP channels have base and tip radii R b = 200 nm and R t = R b − ∆R = 50 nm, respectively, such that the channel radius is described by R i (x) = R b − x∆R/L i for positions x ∈ [0, L i ] in the channel.Here L i represents the channel length, which varies for different channels in the circuit.The channel connects two reservoirs containing an aqueous 1:1 electrolyte with ionic bulk concentration ρ b = 0.1 mM (so substantially lower than considered in the main text) and with the viscosity η = 1.01 mPa • s and the electric permittivity ε = 0.71 nF • m −1 of water.We assume a uniform surface charge density eσ = −0.0015enm −2 on the channel walls, resulting in a surface potential ψ 0 ≈ −10 mV and an electric double layer that screens the surface charge with Debye length λ D ≈ 30 nm.At the tip the Debye length is not much smaller than the channel radius, which does not fully satisfy the assumption of a small Debye length compared to the channel radius made in the theoretical framework we use to model the channel conductance [3].However, it was shown that the dynamic conductance properties are still reasonably well described by the theoretical model we lay out below [2].The ions in this instance are assumed to all have diffusion coefficients D = 1.75 µm 2 ms −1 .The electro-osmotic volumetric fluid flow rate Q i (V i ) that is driven by an applied voltage V i over the channel is accurately represented by its linearresponse approximation . Our predictions of the dynamic conductance of UP channels is based on an analytical model that describes their steady-state [4] and dynamic [2] conductance properties.UP channels exhibit a voltage-dependent concentration polarisation that can reasonably accurately be described analytically by a slab-averaged salt concentration profile where ∆g ≡ −2w(∆R/R b )Du with w = eDη/(k B T εψ 0 ) the ratio of ionic to electro-osmotic mobility [? ] and the tip Dukhin number Du = σ /(2ρ b R t ).The temperature is set at T = 293.15K throughout.For the negative surface charge we consider here on the channel walls, Eq. (S1) describes an enhancement of the steady-state salt concentration in the channel at V i < 0, and a decrease at V i > 0. Following Ref. [2], we approximate the static conductance g i,∞ (V i ) to follow from the salt concentration profile according to This is a simplification compared to the more accurate dependence on L i / L i 0 (ρ i,s (x,V i )) −1 dx used in Eq. (3) in the main text, although Eq.(S2) was also found to still work reasonably well [2,4].Here g i,0 = (πR t R b /L)(2ρ b e 2 D/k B T ) is the conductance of the UP channel in equilibrium.Eq. (S2) combined with Eq. (S1) predicts that the conical UP channel is a current rectifier since, for surface charge σ < 0, a negative voltage V i < 0 increases g i,∞ (V ) while a positive applied voltage V i > 0 decreases g i,∞ (V ), with respect to g i,0 .The difference is caused by salt accumulation and depletion in the channel for V i < 0 and V i > 0, respectively.
The typical time it takes for salt to accumulate or deplete, and hence the typical conductance memory retention time of the channel, is independent of the surface charge [2] and identical to the time scale of the BP channels of the main text given by The dynamic conductance g i (t) of channels of type i is, similar to BP channels, well-described by the differential equation where V i (t) is the time-dependent voltage drop between base and tip of the channel.

Tonic bursting with unipolar channels
To investigate the emergence of bursting behaviour we consider the same circuit as in the main text, but with UP channels, different circuit parameters and three super slow channels (each connected in series to individual batteries) connected in parallel, which is mathematically equivalent to tripling the conductance of a single super slow channel.The two fastest and shortest channels have an individual battery in series with potential E ± = ±0.975V, and are of length L ± = 8 9 µm ≈ 0.89 µm, hence with timescale τ ± ≈ 0.038 ms.The slow channel is connected in series to a battery with potential E s = −0.5 V, with the channel length L s = 27.75 µm and thus a timescale τ s ≈ 37 ms.Finally, the three super slow (and even longer) channels that are placed parallel of each other each have an individual battery in series with potential E ss = −0.5 V and a channel length of L ss = 277.5 µm, resulting in a timescale of τ ss ≈ 3.7 s.From now on, these three channels are denoted as a single super slow channel with base conductance 3g 0,ss , which is mathematically equivalent since they are described by identical equations.The only additive value of three channels rather than one is their increased total conductance that is necessary for a chattering spiking pattern to emerge in this instance.The capacitor has capacitance C = 4 fF and a stimulus current I(t) can be imposed through the circuit.Via Kirchhoff's law we find the following equation to describe the time-evolution of the voltage V m (t) over the capacitor,

C
dV m (t) dt where i ∈ {+, −, s, ss} and with the voltage arguments V i (t) over the channels given by V − (t) = V m (t) − E − , V + (t) = −V m (t)+E + , V s (t) = −V m (t)+E s , and V ss (t) = −V m (t)+E ss .By numerically solving for V m (t) for sustained current stimuli that undergo a step from 0 pA to 1.495 pA and from 0 pA to 1.485 pA, we find the voltage responses shown in Fig. S2(a) and Fig. S2(b), respectively.In Fig. S2(a) we see after a transient has passed upon applying the stimulus that the circuit settles into tonic spiking.In Fig. S2(b) we show that tonic bursting appears after a similar transient as in Fig. S2(a).Therefore we can also produce the additional spiking mode of tonic bursting with UP channels.However, the typical potentials during spiking are further removed from typical voltages in neurons [6] compared to the results in the main text.Additionally, the higher battery potentials and lower salt concentration and capacitance in this instance are not as similar to their biological counterparts [7][8][9] as the values used in the BP channel circuit discussed in the main text.Lastly, while the approach of connecting three super slow channels in parallel is in principle physical, it further complicates the circuit design.Nevertheless, the existence of tonic bursting in a similar circuit with different fluidic memristors does show that bursting spiking modes can be achieved using multiple types of iontronic devices, suggesting that further improvements can be achieved by using even more desirable devices than the BP channels used in the main text.

FIG. 3 .
FIG.3.Evolution of the three (relevant) dynamic variables V m (t), g s (t), and g ss (t) during the limit cycle of the tonic bursting case shown in Fig.2(e).(a) Dynamic variables V m (t) and g s (t) during a single spike.At low g s (t) the potential V m (t) can switch to a positive voltage.Due to the delayed increase of g s (t) after a timescale τ s in response to the increase of V m (t), the conductance g s (t) drives V m (t) down again, forming the spike.(b) Dynamic variables V m (t) and g ss (t) during a burst of spikes.After three spikes g ss (t) increased sufficiently, in response to the increase of V m (t) during the spikes, to halt the spiking.Without spiking, g ss (t) decreases over a timescale τ ss after which the spiking starts again, creating tonic bursting.(c) The three dynamic variables V m (t), g s (t), and g ss (t) in a phase portrait during tonic bursting, showing a characteristic bursting phase diagram with three loops (blue) where g ss (t) changes little, corresponding to the three spikes, and a trajectory connecting the third and first loop (red) where g ss (t) returns from high to low, corresponding to the periods of quiescence.

Figure S1 .
Figure S1.The voltage dependence of the steady-state conductance g ∞ (V ) as predicted by our analytic approximation of Eq. (3) (red) and by finite-element calculations of the full Poisson-Nernst-Planck-Stokes equations (blue), as detailed in Ref. [1], for a bipolar channel of length L = 10 µm, base radius R b = 200 nm, tip radius R t = 50 nm, and salt concentration ρ b = 2 mM.
and is conveniently characterised by the Péclet number at the narrow end Pe

Figure S2 .
Figure S2.Two modes of voltage spiking extracted by modeling one and the same iontronic circuit with UP channels driven by different time-dependent currents, with (a) tonic spiking, i.e. regular spiking, for a stimulus of 1.495 pA, and (b) tonic bursting, i.e. a short burst of spiking alternating with periods of quiescence, for a stimulus of 1.485 pA.System parameters are given in the text.