Robust spike timing in an excitable cell with delayed feedback

The initiation and regeneration of pulsatile activity is a ubiquitous feature observed in excitablesystems with delayed feedback. Here, we demonstrate this phenomenon in a real biological cell. We establish a critical role of the delay resulting from the finite propagation speed of electrical impulses in the emergence of persistent multiple-spike patterns. We predict the coexistence of a number of such patterns in a mathematical model and use a biological cell subject to dynamic clamp to confirm our predictions in a living mammalian system. Given the general nature of our mathematical model and experimental system, we believe that our results capture key hallmarks of physiological excitability that are fundamental to information processing.


Introduction
Persistent and repetitive electrical spiking activity underlies information processing in a wide range of excitable biological systems, including neurons [1,2], arthropod muscle fibre [3], cardiac cells [4][5][6], pancreatic β-cells and other endocrine cells [7,8]. Spikes, which are brief, rapid depolarizations of the voltage across cell membranes, are initiated due to the excitability properties of the cell and in response to incoming stimuli, such as electrical impulses induced by other cells. Excitability properties allow cells to selectively respond to inputs in an 'all-or-none' fashion, acting as a nonlinear filter and enabling precise temporal encoding of stimuli. In certain contexts, spiking patterns induced by transient stimuli are repetitively regenerated following stimulus withdrawal [9], suggesting a mechanism by which memories could be encoded in the timing or frequency of spiking events, as has been postulated for the encoding of sensory inputs and in motor control [10,11].
Electrical impulses in the brain traverse, via axonal and dendritic projections, distances that can extend many times the cell body. The finite conduction velocity along such projections induces delays in signal transduction and information processing, including in the encoding and decoding of memories [12][13][14]. Since delayed spiking activity has been linked to working memory [14], it should come as no surprise that demyelinating diseases, in which the electrical insulation around axons degrades, are associated with memory impairment.
Excitable systems subject to delayed feedback are capable of regenerating their own excitable response, resulting in robust multi-pulse patterns dependent on the feedback strength and the length of the delay interval [15][16][17]. The presence of rich dynamics, including coexistence between different types of dynamic behaviour in nonlinear systems with delays [15][16][17], has sparked significant interest in investigating the interplay between excitability and delay in a variety of systems. A prominent example of such a system involves semiconductor micro-cavity lasers [18][19][20]; in an analogous fashion to neural systems, the number of emission pulses per delay period in the laser system is thought to be reflective of its ability to store information. In particular, it has been suggested that this property could be used in the construction of bioinspired 'neuromorphic' photonic resonators that process information through light pulses alone [21,22]. Supporting this perspective are results that show that pulsetiming patterns due to regeneration of excitation exist in a coherently driven laser [23] and a wavelength tuneable semiconductor laser subject to optoelectronic delayed feedback [24] as well as in a driven laser subject to optoelectronic delayed feedback [25], and a micro-laser with integrated saturable absorber [26]. Given the similarities in excitability between neural and laser systems, it is pertinent to inquire as to which features of the regenerative process are common to both and which differ, noting that such dynamics have previously been observed in delay-coupled recurrent neural loops subject to continuous long-term stimulation [27,28].
A mathematical description of neuronal excitability by Hodgkin & Huxley [29] paved the way for theoretical investigations of the nonlinear dynamics associated with neural response. Since then, the essence of excitability has been refined into simpler models, such as the FitzHugh-Nagumo model [30,31] and Morris-Lecar models [3]. Theoretical analysis of such models offers a powerful approach to gain insight into the generation of rhythmic activity in electrically excitable cells by identifying and quantifying the contribution of specific processes to the overall dynamic response. Since the refined models capture the key features of excitability, observations from such analysis can be mapped to a wide range of different cell types.
In this paper, we combine theory and experiments to systematically investigate the interplay between excitability and delayed feedback and its implications for information processing. We append a delayed feedback term to a wellestablished, widely used mathematical model of cellular excitability and use numerical continuation (bifurcation analysis) to characterize the behaviour of the model as the feedback strength and delay period of the feedback are varied. The form of feedback we incorporate has a form akin to that of an autapse, which is a chemical synapse from a neuron to itself. Autaptic connections have previously been shown to be important in establishing and maintaining rhythmic activity in mathematical models of small [32] and large [33] neural networks. Here, we investigate the interplay between autaptic processing and excitability dynamics in a single cell. To test our theoretical results, we employ dynamic clamp [34] and confirm our theoretical predictions in a living mammalian system represented by an excitable cell as a natural platform for studying excitation and regeneration patterns.

Theoretical predictions
We adopt the canonical Morris-Lecar (ML) model [3] to describe somatic membrane potential dynamics and extend it with a delayed feedback term to describe a transient input stimulus to the model cell. The ML model belongs to the general class of Hodgkin-Huxley-type models, which form the basis for almost all biophysical characterizations of electrical excitability in cell membranes [29]. The feedback term incorporates a fixed delay that could be seen as a representation of conduction delay of the voltage signal, for example. We assume negligible attenuation of the signal and that the delay in the system is due to finite conduction velocity and synaptic transmission. Based on these assumptions, we describe the delayed feedback via a current-based model of post-synaptic activity with I s (t) = κs with the dynamics for s given by [35] Here, V s is the threshold for synaptic activation, V h is the sensitivity of the synapse around the threshold and V pre is the voltage of the pre-synaptic cell, which is the same as the post-synaptic cell in our one cell system. For simplicity, we neglect synaptic dynamics and assume that s is at quasisteady state but note that the synaptic delay may still be incorporated into τ. These observations allow us to write the current induced by the delayed feedback as where the parameter κ scales the feedback strength. This current is added to the standard ionic and applied currents to produce the following complete model: Parameter definitions and values are listed in table 1.

Bifurcation analysis
To identify regions of spiking activity in the ML model with delayed feedback, we appeal to bifurcation theory, which allows for the rapid classification of dynamic behaviour as model parameters are varied. The bifurcation diagram for the system under variation of the feedback strength, κ, and delay, τ, is shown in figure 1a, where we observe coexistence of stable spiking periodic orbits for sufficiently high κ and τ. Along the blue line, we fix κ = 60 and vary τ, upon which we observe a sequence of bifurcations (folds of periodic orbits) with increasing τ, each of which gives rise to a new pair of spiking solutions, one of which is stable. Figure 1b demonstrates that the branches of stable spiking solutions have approximately equal amplitude (except close to the bifurcation at which they originate), concordant with the intuition that action potentials in an individual cell are stereotyped events so that information is conveyed through spike timing rather than amplitude. By contrast, figure 1c highlights that the period of the spiking solutions varies significantly between different branches, while depending approximately linearly on the delay τ along each branch. In particular, each stable branch can be mapped to a solution with n [ N spikes per royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20210029 delay interval, with a common interspike interval (ISI) between each spike. The number of spikes that can be fitted within the delay interval is limited by the minimal period of spiking activity and the neuron's refractory (or recovery) period (since new action potentials cannot be generated during this time). This imposes limits on the information that can be encoded if one regards spikes as bits of information as suggested in [21,22]. Examples of stable one-spike ( per delay interval) and two-spike ( per delay interval) solutions for κ = 60 and τ = 400 are shown as voltage profiles over two delay intervals in figure 1d,e, respectively. The solution to which the system ultimately converges is dependent on the initial condition (given by the history over the delay interval). For a cell initially at rest, this is determined by the timing and form of the current stimulus to which it is exposed. The insets in figure 1d,e depict the Floquet multipliers, which determine the (linear) stability of periodic solutions. We find that the two-spike solution, while indeed stable, possesses a nontrivial multiplier close to the unit circle (the boundary for stability); this scenario is repeated for other n-spike solutions. Hence, equidistant spiking solutions are weakly stable and approached slowly by trajectories over long periods of time.
To examine the effect of multistability on memory encoding, we perform numerical simulations. Here, we use a paired-pulse protocol in which two square-wave current pulses with a finite inter-pulse interval are applied to the cell, as depicted by the purple trace in figure 1f,g. Following these initial pulses, no other external inputs are applied to the system. In figure 1f , the first current pulse elicits a full amplitude spike, but the second does not, because it falls into the refractory period of the action potential. During the second delay interval after the initial excitation, the response of the system to the second current pulse has decayed entirely and the trajectory settles rapidly into a one-spike solution. In figure 1g, the interval between the two current input pulses is well beyond the refractory period, and they both elicit large-amplitude action potentials that persist in subsequent delay intervals with heterogeneous ISI patterns (i.e. non-equispaced in time) approximately conserved across many subsequent delay intervals. Nevertheless, the effect of the weak convergence to the stable solution (figure 1e) is evident: ultimately, the neuron settles to the stable solution with two equally spaced spikes in a similar manner to that observed in laser systems [20].

Experimental verification
We next seek to test our theoretical predictions in an excitable mammalian cell (a GH4 cell; see Material and methods) by using the dynamic clamp technique to control the delay time and feedback strength. A dynamic clamp is an electrophysiological technique used to manipulate the excitability properties of individual cells [34,36]. It is typically used to modify the activity of ion channels (the primary drivers of spiking behaviour) according to pre-defined mathematical models, which allows one to quantify their contribution to overall excitability [37]. Here, we use it to apply delayed feedback to an isolated cell. In this way, the spike generating mechanisms of the cell are preserved, while parameters of the feedback are directly modifiable during the experiment. A schematic of the dynamic clamp protocol is shown in figure 2a. Briefly, the voltage signal recorded via the amplifier is used as an input to our mathematical model of the delayed feedback (2.2). The resulting feedback current is sent back to the amplifier to be injected into the cell to provide real-time closed-loop feedback.
Synthetic data from the augmented ML model are shown in figure 2b, in which we apply a solitary current pulse to the cell of sufficient amplitude and duration to trigger a spike. In the unshaded region, the feedback term is absent (i.e. the synaptic gain is κ = 0). Following the first spike, which occurs in direct response to the current pulse, the cell returns to rest. In the shaded region, the gain is now set to κ = 60. Now, after the application of the initial current pulse and subsequent spike, recurrent excitation builds up in the cell and self-sustained spiking activity is generated with a period approximately equal to the delay interval.
We experimentally emulate the paired-pulse protocol performed in the mathematical model in figure 1f,g and  figure 2b in a real cell. The results of a typical experiment are shown in figure 3, with full traces displayed in electronic supplementary material, figures S1-S3. Figure 3(a,c,e) shows the first 8 s of the overall recording, including the initial current pulse. To facilitate viewing, the green and orange arrows above the spikes indicate by which initial current pulse they were originally generated; figure 3(b,d,f ) shows the same data in a pseudo-space plot where the respective time series is folded over the delay interval to show the evolution of pulses in the delay interval from cycle-tocycle. In figure 3a,b, the synaptic gain is set to κ = 40, which is not large enough for the initial current pulses to generate self-sustained activity, and the cell returns to rest. Upon setting κ = 60, the spikes produced by the initial current pulse (indicated with an orange arrow) die out, while those produced by the second current pulse (indicated with a green arrow) grow in amplitude until self-sustained spiking activity with one spike per delay interval is obtained, as shown in figure 3c,d, echoing the findings in figure 1f. Keeping κ = 60 but increasing the inter-pulse interval results, as reported in figure 3e,f, in the growth of action potentials generated by both initial current pulses until recurrent activity with two spikes per delay interval is established, thus recapitulating the result from figure 1g. This experiment highlights the capability of the cell to encode distinct and persistent multi-spike rhythms (activity patterns) and confirms the predictions from the analysis of the mathematical model.

Discussion
In this study, we demonstrated that a single excitable cell with delayed self-feedback is a rhythm-generating system capable of sustaining regular spiking oscillations of different periods in response to an initial excitation pattern. The number of distinct patterns that can be realized by the neuron is heavily dependent on the strength and delay interval of the self-feedback. In particular, if the feedback strength is too low, the excitation in the cell cannot be regenerated and no persistent activity is generated. If the feedback strength is sufficiently high, the number of supported patterns increases with the delay interval. This latter observation results from the interplay between excitability properties and the delay: the number of spikes that can be accommodated over a fixed delay is limited by the cell's spiking and refractory periods. Mathematical analysis of such dynamic behaviours makes these observations rigorous, whereupon we find that the generation of each spiking pattern can be attributed to a specific bifurcation.
Which of the potentially many spiking patterns is ultimately observed depends on the initial stimulus. The mathematical model predicts that changing the inter-pulse interval in a paired-pulse stimulus experiment is sufficient to switch the resulting activity pattern from a one-spike per delay interval to a two-spike per delay interval solution. The corroboration of these results in a real cell confirms these predictions and validates the modelling approach.
In the experimental verification of our model predictions, we showed that the delayed feedback of the form chosen is sufficient to produce stable rhythms with one-spike and twospike per delay period. Our bifurcation analysis demonstrated that the associated periodic solutions possess Floquet multipliers close to (but inside) the unit circle, suggesting slow rates of attraction to equidistant pulsing rhythms. Additionally, biological systems display considerable complexity and heterogeneity. In particular, neurons are well known to exhibit different classes of excitability [38] and may generate complex spiking rhythms even in the absence of external input or feedback [39,40]. It is therefore a pertinent question for careful future attention to inquire how robust such multiplespike solutions are in real cells. In this direction, we make two remarks. Firstly, our model of choice is general in that we did not tune our analysis to any particular cell type. Moreover, the ML model can, in different parameter regimes, exhibit properties of the different excitability classes. As such, our framework provides flexibility to explore such avenues. Secondly, while the repetitive activity of the kind we discuss here may be fragile in certain experimental systems, recent work in laser physics found sufficiently large parameter ranges over which different multiple pulse rhythms due to excitability and delayed feedback can be found stably over very large timescales [20,41]. This suggests that these dynamics may be quite robust in practice also when the excitable element is a cell of suitable type. A more detailed investigation of how multiple-spike rhythms of the clamped cell experiment depend on the properties of the feedback loop, as well as those of the cell itself, is a challenging task and the subject of ongoing research.
Given the nature of the Hodgkin-Huxley-type model and of the experimental setup, we believe that our results capture key hallmarks of physiological excitability fundamental to neural processing in a wide array of contexts. Notably, the striking similarity between the theoretical predictions and experimental data was not produced by fitting a specific mathematical model to the chosen cell type, but was reflective of general properties of excitability that are present in most neurons, as well as other cell types [4][5][6][7][8]. The parallels of our results with those observed in laser systems [18,20] further highlight the ubiquitous nature of the regeneration of activity in excitable systems with delay. Our experimental setup showcases that a living mammalian system, represented by an excitable cell, is a feasible platform for studying this phenomenon in a neural context. The framework we have established is flexible and will allow further exploration of the interplay between delay and important biophysical aspects such as the aforementioned role of different excitability sub-classes, and of complex intrinsic spiking patterns.
The ability to control multi-spike behaviour in biological cells more widely opens a new avenue towards investigating how memories can be robustly encoded. This may well have significant implications for better understanding of not only normal brain function [42], but also debilitating conditions such as Alzheimer's disease [43,44]. Knowledge of the interplay between excitability and delay is also crucial for designing effective stimulation protocols for brain-machine interfaces. We expect that these issues, among others, will dictate how delays can be incorporated in useful devices based on excitable systems, whether in the neurological, sensory, motor or other domains. In the unshaded region, the feedback is inactive (κ = 0). Here, a current pulse induces a single spike, after which the cell returns to rest. In the shaded region, the feedback is active (κ = 60) and repetitive spiking activity with a period approximately equal to the delay interval is observed even though only one input current pulse is issued. The delay, τ, is marked in the upper trace for reference.

Cell culture
The rat norvegicus pituitary tumour GH4C1 (GH4) cell line was cultured in Hams F-10 nutrient mixture (81.5%), supplemented with horse serum (15%), FBS (2.5%) and Glutamax (1%). Cells were thawed following standard procedures and cultured for 2 weeks in a 5% CO 2 incubator at 37°C, passaging every 3-4 days in T75 vented flasks. Prior to experimentation, cells were transferred to 15 mm glass cover slips for patching, seeded at a sufficiently low density to easily identify isolated cells. Immediately prior to experimentation, cells were washed three times in PBS to remove any residual serum.

Patching solutions
The extracellular solution was comprised of 2.8 mM glucose, 132 mM NaCl, 5 mM KCl, 2.6 mM CaCl 2 , 1.2 mM MgCl 2 , 10 mM HEPES, 100 μM paxilline ( pH 7.4 adjusted with NaOH, osmolarity 295 mOsm l −1 adjusted with sucrose). Paxilline was added to the solution to block BK channels to reduce the probability of bursting activity typically seen in GH4 cells [47]. The intracellular solution was comprised of 100 mM K-gluconate, 20 mM KCl, 4 mM NaCl, 1 mM CaCl 2 , 4 mM MgCl 2 , 10 mM EGTA, 10 mM HEPES ( pH 7.2 adjusted with KOH, osmolarity 285 mOsm l −1 adjusted with sucrose). Immediately before patching, amphotericin B (Sigma) was added to the internal solution at a final concentration of 0.6 mg ml −1 . When not in use, the internal solution was kept on ice and away from light.

Dynamic clamp
The dynamic clamp was achieved through the Teensy 3.6-based breadboard circuit as presented in [34] with updates as presented at the homepage for that project (www.dynamicclamp.com/ updates). Computer code to simulate the delayed synapse model was developed in Arduino and can be downloaded from the GitHub repository (https://github.com/Slowinski-Piotr/MorrisLecarDDE). To account for the more depolarized membrane potential of the GH4 cells relative to the ML model, the value of V s in the experiments was set to −20 mV. The Teensy device was controlled by bespoke software developed in Matlab, which is freely available at www.github.com/Kyle-Wedgwood/DynamicClampController. In each experiment, an initial current pulse is followed after some time by a second current pulse. In (a) and (b), κ = 40, while κ = 60 in all other panels. Across the entire experiment, the delay is fixed at τ = 800. The inter-pulse interval was set to 600 ms in (a) and (e) and 400 ms in (c) (and the same for the corresponding plots to their right). The orange and green arrows above the spikes indicate which initial current pulse induced the spiking behaviour. Observe that in (c), the spikes induced by the orange pulse, alone, die out, while the ones induced by the green pulse grow and ultimately give rise to a self-sustained oscillation with one spike per delay interval. This is in contrast with (e) in which the spikes induced by both the orange and green pulse are sustained, resulting in a two-spike per delay interval pattern. (b,d,f ) The same data as (a,c,e) but with the time series plotted against successive cycles of the delay interval.

Protocol
On the day of recording, thick-walled patch pipettes were pulled to resistance of approximately 3.0 MV and subsequently firepolished. Cover slips were loaded into the recording chamber and allowed to settle for 5 min under constant perfusion at 1.5 ml min −1 in external solution. During this period, isolated cells were marked as candidates for patching. Following the settling period, pipettes were front-filled for 10 s with internal solution without amphotericin B and subsequently back-filled with amphotericin B containing internal solution. Perforated wholecell patch-clamp was then achieved within 60 s of filling the pipette, with seals of >5 GV achieved. Perforation was continued until the series resistance fell below 20 MV, which took around 15 min. Following successful perforated patch, the amplifier was switched to current-clamp. Bridge balance was applied to compensate for residual series resistance. Initial gap-free recordings were obtained to verify that cells were capable of spiking activity. Following this, a holding current of −30 pA was applied to silence spontaneous spiking activity. Paired pulses of 10 ms duration and 50 pA amplitude with varying interpulse intervals ranging from 200 to 600 ms were applied to the cell. Following the initial pulses, recordings were continued for a further 30 s. Recordings were taken at a temperature of 34-36°C with a sampling frequency of 20 kHz and filtered with a Bessel filter at 2.8 kHz. Throughout the recordings, the delay period was set to 800 ms and the feedback gain was set to either 40 or 60.
Data accessibility. Data and computer code related to the mathematical model and dynamic clamp experiments can be downloaded from the GitHub repository https://github.com/SlowinskiPiotr/ MorrisLecarDDE.