Closed-loop stabilization of the Jamming Avoidance Response reveals its locally unstable and globally nonlinear dynamics

Weakly electric fish emit an electric organ discharge (EOD) that is used for electrolocation (Caputi and Budelli, 2006; von der Emde, 2006) and communication (Hopkins, 1974; Fortune, 2006; Hupé and Lewis, 2008; Triefenbach and Zakon, 2008). In wave-type electric fish, the EOD is quasi-sinusoidal and has a relatively stable baseline frequency when undisturbed (Bullock, 1969; Bullock et al., 1975; Moortgat et al., 1998; Zakon and Dunlap, 1999; Zakon et al., 2002). When two or more fish are in close proximity, their EODs interact to produce emergent amplitude modulations (AM) at the difference frequency (df) between the fish. Eigenmannia virescens shift their baseline EOD frequency in response to low (<15 Hz) dfs (Watanabe and Takeda, 1963; Bullock et al., 1972), which have been shown to impair aspects of electrolocation (Heiligenberg, 1973; Bastian, 1987). The direction of the frequency shift is determined by the sign of the df and results in an increase in the magnitude of the df. This behavior is known as the Jamming Avoidance Response (JAR).

The JAR has been analyzed at all levels of organization, from whole-organism behavior down to specific ion channels (Fortune and Rose, 2003; Heiligenberg, 1991; Bastian and Heiligenberg, 1980; Heiligenberg and Bastian, 1980). Despite the fact that the JAR is among the best understood sensorimotor circuits, the sensorimotor responses have not been modeled as a dynamical system. One challenge to modeling the temporal dynamics arises from the intrinsically ‘unstable’ nature of the JAR. This is because the fish shifts its EOD frequency (f₁) in the direction away from the frequency of the conspecific (f₂), resulting in an increase of |df|, where df=f₂–f₁. We overcame this challenge by closing a feedback loop around the natural behavior: the frequency of a conspecific-like signal was calculated and adjusted in real-time to stabilize the response and drive it to any desired frequency in a neighborhood of the fish's original baseline frequency.

Perturbation experiments on the stabilized closed-loop system were used to characterize the dynamics of the JAR. These perturbations included sinusoids, sums of sinusoids, chirps and band-limited noise. Responses to these stimuli were used to estimate a non-parametric frequency response function (FRF). The FRF was then used to infer the frequency response of the open-loop behavior, i.e. the JAR itself. A first-order delayed parametric model was fitted to the behavior near its equilibrium.

This local model does not, however, capture the nonlinear features of the behavior: the biological relevance of the JAR lies in its escape from the unstable equilibrium. To address this, we extended the linear model using additional experiments in which the dfs were ‘clamped’ to furnish a complete nonlinear model. The nonlinear model was parsed into terms that capture competing avoidance and return responses and was validated by comparison with responses to open-loop stimuli. The model was also used to predict a saddle–node bifurcation in the vector field of the system. The bifurcation causes one stable equilibrium to vanish, ultimately manifesting as the ‘snap-through’ of the fish's frequency (f₁) from one stable equilibrium (low |df|) to another (high |df| of opposite sign).

Adult Eigenmannia virescens (Valenciennes 1836) (10–15 cm in length, EODs between 346 and 452 Hz) were obtained from commercial vendors. The fish were housed in group aquarium tanks that had a water temperature of ~27°C and a conductivity in the range 150–500 μS cm⁻¹ (Hitschfeld et al., 2009). All experimental procedures were approved by the Johns Hopkins Animal Care and Use Committee and followed guidelines established by the National Research Council and the Society for Neuroscience.

The experimental tank was maintained at a temperature of 25±3°C and a conductivity of 150±25 μS cm⁻¹. Each fish (N=7; n=5 reference-tracking trials, n=2 clamp trials) was tested individually. Each fish was acclimated to the testing tank for a minimum of 24 h prior to the start of the experiment. After the initial acclimation period the fish was restricted in a chirp chamber for 2–3 h to allow the EOD frequency to stabilize. The chirp chamber served to restrict the movement of the fish and prevent changes in orientation during the experiment, resulting in more consistent measurements of the EOD frequency. A pair of measurement electrodes (red) were placed longitudinally (near the fish's head and tail) to record the EOD and a second pair of stimulus electrodes (black) were placed transverse to the fish to provide a frequency-controlled sinusoidal stimulus (Fig. 1A). The distance between the electrodes in each pair was 25 cm.

The EOD of the fish, recorded via the head-to-tail electrodes (Fig. 1A; red circles), was filtered and amplified (0.1 Hz to 1 kHz bandpass, gain 100; A-M Systems Model 1700, Sequim, WA, USA) and input to a frequency-to-voltage converter (F2V; FV-1400, Ono-Sokki, Yokohama, Japan). The F2V calculates the frequency of the signal using measured time differences between zero crossings, digitally averages these values over four cycles, and outputs a voltage that varies linearly with frequency over a set range (e.g. the range ±5 V would correspond to 475–525 Hz). This creates a time lag of no more than four cycles of the EOD. The F2V output was further filtered (Chebyshev low-pass, 30 Hz cutoff). All told, these delays correspond to a worst-case phase lag of 5.2 deg relative to the highest reference stimulus frequency (1 Hz), at the lowest EOD baseline (346 Hz) tested. Both the amplified signal and the filtered F2V output were fed into a Power1401 Mk.II signal acquisition device (CED, Cambridge, UK) which ran a custom-written sequencer script that read the input signal, performed the feedback calculation and generated a sinusoid with the desired output frequency. Using this setup enabled regular temporal sampling intervals and provided more deterministic computation time than including a standard computer in the feedback. The signal acquisition device received parameters and reference signal for each trial from the Spike2 software (CED), which ran simultaneously on a computer. During the trials, this software received and recorded data from the input, output and intermediate channels. The amplitude of the stimulating sinusoid was 100 μV cm⁻¹ (unless otherwise noted), which produced 14–43% contamination across individuals (EOD amplitudes 232–715 μV cm⁻¹) measured at a 1 cm dipole placed adjacent to the head of the fish. We performed an identification experiment on the feedback system to assess what its characteristics were, especially the delay introduced by the equipment. The contribution was minimal (~2 ms delay) due to fast instrumentation. As such, the feedback delay was not incorporated into subsequent calculations.

The fish EOD frequency f₁(t) and the stimulus frequency f₂(t) are functions of time. Under constant lighting, temperature and conductivity, and without conspecific stimulation, the EOD frequency remains relatively stationary over long periods of time (Bullock, 1969). The initial time t=0 for each trial was preceded by a period of no stimulation for at least 300 s, and we defined the baseline (initial) frequency f₁(0) as the EOD frequency at that initial time.

Sinusoidal stimuli were of frequencies 0.01, 0.055, 0.215 and 0.995 Hz and of durations 1000, 1000, 500 and 200 s, respectively. The stimulus durations were chosen to have a sufficient number of beat cycles for spectral analysis.

These stimuli were the sum of 10 logarithmically spaced sinusoids with randomized phase, in the range 0.01–1 Hz. The sum-of-sine stimuli included the four single-sine frequency components. The stimulus duration was 1000 s.

The chirp stimulus was a sinusoid of increasing frequency, from 0.01 to 1 Hz over 1000 s. The increase of frequency was exponential, ensuring sufficient stimulus power across all frequencies. Several weakly electric fish also exhibit EOD frequency changes termed ‘chirps’, not to be confused with the chirp stimulus used in this paper.

The long-chirp stimulus was similar to the chirp stimulus except that the frequency increased from 0.001 to 1 Hz. Consequently, the stimulus duration was increased to 10,000 s.

Band-limited noise stimuli consisted of non-overlapping, 2 Hz wide frequency bins from 0 to 20 Hz. For each bin, a stimulus was generated by summing together sinusoids associated with all uniformly spaced frequencies as dictated by the sampling rate. This ensured that the stimulus has uniform power over all the frequencies analyzed within each bin. Each component had randomized phase and was scaled equally so that the sum would have a maximum magnitude of 1 Hz. The trial duration was 300 s, except for the lowest frequency band of 0–2 Hz, which was 1000 s long.

For closed-loop reference-tracking trials, we provided a reference signal r(t) and applied the controller as described in Eqn 2. At the start of each trial, the fish's baseline was measured. Subsequently, there was a ‘balancing’ period of 100 s wherein the controller aimed to keep the fish at the initial EOD frequency [r(t)=0]. To avoid startling the fish, the amplitude of the stimulus signal was ramped up linearly from 0 to 100 μV cm⁻¹ (unless noted) over the first 50 s of the balancing period. After the balancing period, a stimulus (single sine, noise, etc.) was introduced as the reference. Trials were pseudo-randomized. All stimuli were pre-generated at 1 kHz and input to the sequencer at the start of the trial. A sample interval of a closed-loop trial with a single-sine stimulus is shown in Fig. 2.

Each fish (N=5) completed closed-loop reference-tracking trials (n=49) in a randomized order within a single testing session. The trials consisted of the following: (i) single sine stimuli (n=8), the four frequencies of which were replicated for magnitudes of 1 and 2 Hz; (ii) sum-of-sines stimuli (n=4), with two different component magnitudes (0.2 and 0.3 Hz), with two different sets of randomized component phases each; thus, the maximum stimulus magnitude was 2 or 3 Hz; (iii) chirps (n=2), with magnitudes of 1 and 2 Hz; (iv) white noise stimuli (n=30), with three identical replicates of each of the 10 frequency bands; (v) chirp (n=2) stimuli with signal amplitude 50 and 200 μV cm⁻¹ (typical value=100 μV cm⁻¹); these were to examine sensitivity of the identified system to signal amplitude; and (vi) chirp (n=3) stimuli with controller gains 1.5, 2.5 and 3 (typical value=2); these stimuli were, similarly, to examine sensitivity to the feedback gain. Sample responses to four stimulus types are shown in Fig. 3.

In addition, a subset of three individual fish were presented with long chirp stimuli at two magnitudes (1 Hz, 2 Hz). After each trial, the stimulus was turned off and the EOD frequency was allowed to stabilize over a period of 300 s before the next trial began.

Instead of driving the fish frequency towards a goal, the clamp trials applied a stimulus such that d(t) was maintained at a desired value. Two types of clamp trials were run: (i) static clamps, where the clamp was a constant value d=d_s for 300 s; and (ii) dynamic clamps, where the clamp was kept at a particular value d=d_s for 100 s and then oscillated around the value according to a reference trajectory r(t), such that d(t)=d_s+r(t).

Static- and dynamic-clamp experiments were performed on N=2 individuals that did not complete the reference-tracking trials. These included n=39 static clamps with d_s from −50 to 50 at higher resolution closer to 0. The dynamic clamp trials had single sine, sum-of-sines and chirp stimuli as the reference trajectories.

The objective of these trials was to observe the response of the fish to a stimulus whose frequency trajectory u(t) was predetermined, and not tied to y(t). Two types of open-loop trials were performed. (i) Step inputs: initially, the fish was driven towards 0 Hz in closed loop for the balancing period of 100 s. The amplitude was linearly ramped up to 100 μV cm⁻¹ for 50 s as described previously. At 100 s, the stimulus switched to open loop, and u(t) was commanded to a fixed value for a further 100 s. We performed trials at steps with magnitudes from −5 to +5 Hz. (ii) Ramp inputs: the magnitude of u was ramped down from +30 to −30 Hz or up from −30 to +30 Hz in 300 s at a constant rate of change of 0.2 Hz s⁻¹. The stimulus amplitude was ramped to 100 μV cm⁻¹ in the first 10 s. The initial difference frequency was large, |df|≈30 Hz, well outside the typical range of the JAR.

All data analysis was carried out using scripts custom-written in MATLAB (The MathWorks Inc., Natick, MA, USA). For each trial, we recorded the reference, EOD, F2V output, and output waveform, sampled at 10 kHz. The voltage signal from the F2V was scaled and offset to convert it into frequency. Extremely rapid, transient changes in frequency (commonly caused by fish movement) were eliminated. The known baseline frequency of the fish was subtracted from all signals, so that a frequency of 0 Hz represents the fish's pre-stimulus (baseline) signal. The processed F2V signal was then subsampled to 100 Hz, and used as the output signal y(t). The input for analysis depended on the type of trial: r(t) was used for reference-tracking trials, d(t) was used for the dynamic clamp trials, and u(t) was used for the open-loop trials. These inputs were pre-generated trajectories, as mentioned previously.

The frequency domain representations of the input and output signals were calculated using a fast Fourier transform (FFT) and peaks corresponding to the known number of frequency components in the input were determined (1 for single sines, 10 for sums of sines). The frequency response at a particular frequency was calculated as the ratio of the Fourier transform of the output to the input of the signal at that particular frequency. Thus, we measured eight data points from eight single-sine trials and 40 points from four sums-of-sines trials. Each data point was represented as a phasor, namely a number in the complex plane. The gain (distance of the phasor to the origin) and phase (angle of the phasor from the positive real axis) for all data were computed.

The input and output signals were filtered, subsampled and transformed via FFT as described above. The data points in this case were the input–output ratios of all the frequency components in the chirp frequency range. Thus, with a single trial, we obtained many data points, but each individual data point is somewhat more susceptible to noise. The data points were binned and averaged at 10 bins per decade of frequency, giving us 20 data points per chirp trial and 30 data points per long chirp trial.

White noise stimuli were used to evaluate the range over which the behavior is linear. This was done by evaluating the difference of the square root of the response–response coherence [√(RR)] and the stimulus–response coherence (SR) (Roddey et al., 2000). SR is the coherence between each stimulus and its corresponding response. √(RR) is the square root of the coherence between two responses to the repeated presentation of the same stimulus. In our experiment, as there were three replicates of each frequency band, we obtained three data points per frequency per individual for SR as well as √(RR).

These trials were open loop and were used to compare the fish response with the model response in the time domain. Hence, no frequency domain analysis was done on these trials.

The FRF for single sine, sum of sines, chirp and long chirp closed-loop tracking trials were plotted in a Bode diagram as shown in Fig. 4A. The data were converted to open loop according to Eqn 4 as shown in Fig. 4B. The agreement between single sine, sum of sines and chirp data indicates that the behavior can be modeled approximately as a linear system in the neighborhood of the equilibrium, i.e. near the baseline frequency f₁(0). There are small differences in the closed-loop data among different trial types, particularly between long chirps and the others. Sensitivities of opening the loop tend to amplify these, creating significant differences in the open-loop data. The modest difference in closed-loop responses could come from a variety of sources, including non-linearities and time dependencies in the behavior.

As an important control, we investigated the sensitivity of the open-loop dynamics to experimental parameters. Amplitudes were compared at both half and twice the value of 100 μV cm⁻¹ used in all of the other experiments described in this paper. Gains were tested at values of 1.5, 2.5, 3 and 4 compared with the value of 2 used in all other experiments. There was little-to-no effect of changing either the amplitude of the signal (Fig. 6A) or the feedback gain (Fig. 6B) on the open-loop frequency response. The amplitude insensitivity results in an experimental advantage, as the size of the fish or the specific placement of the fish between the stimulus electrodes would have little effect on the dynamics.

The model fit using the closed-loop experiment is, effectively, a local linearization of the behavior about the unstable equilibrium y=0 (corresponding to the fish's pre-stimulus baseline frequency). Thus, the transfer function described in Eqn 7 is unable to reproduce the full extent of a naturalistic response, which, in addition to avoiding a jamming frequency, involves achieving a steady-state frequency at a higher final df magnitude. Using the linear model as a starting point, and known features of the JAR neural circuit, we fitted a nonlinear model as described below.

The (dimensionless) slopes at the origin for the two curves are remarkably similar (1.68 and 1.73). Both slopes agree well with the predicted value of 1.58 (Fig. 7, orange line), particularly given that the linear system identification data used to predict the escape function slope of 1.58 was categorically different from the clamp experiments, fitting and analytical differentiation used to determine each individual's slope. Note that three trials (Fig. 7, green points), for which the responses were qualitatively different from the other responses, were removed from fitting; including those trials biased the curve toward these outliers, and away from the typical behavior.

As shown in Fig. 7A,B and Eqns 12,13, the escape curves reflect variability between individuals despite the fact that the linearization at baseline is consistent.

As a first validation of the nonlinear model, we examined the dynamic clamp trials, whose FRF should, in theory, match the prediction of the linearization of Eqn 9 at (d_s,y_s). We examined this prediction for two fish (see Fig. 7).

For our prediction, we approximated the nonlinear function e(d) as a straight line through the origin, namely y_s=τkd_s. With this approximation, the FRFs should match the prediction from Eqn 10. The FRFs in Fig. 7Aii,iii and Bii,iii show the frequency responses to dynamic clamps at two points along e(d), marked in Fig. 7Ai and Bi. The response of the linearized model is shown in all four plots and is in good agreement with the frequency responses from the clamp trials. The data shown in Fig. 7 were selected to most clearly illustrate the prediction, but among the experimental paradigms we used, clamp trials were the least robust and repeatable.

For the same two individuals for which the escape functions were fitted, we simulated the model in Eqn 9 in response to (open-loop) step inputs. The simulated step responses matched the data from the open-loop step trials remarkably well, particularly given the qualitatively different nature of open-loop step experiments compared with the experiments used to furnish the nonlinear model. Two examples, showing all the data from one individual for steps u=+4 Hz and u=−3Hz, along with the model response, are shown in Fig. 8.

The model predicts that if the external stimulus frequency were slowly increased or decreased from outside the frequency range typically thought to elicit the JAR, then it should drive the fish's frequency up or down, respectively, until the fish frequency ‘snaps over’ to the other side of the stimulus and then moves away from the stimulus in the opposite direction. Using the open-loop ramp trials (see Materials and methods), we tested this prediction of the model.

We now describe the model's snap-through prediction in more detail. Assuming u(t)=constant, the solution(s) to y(t)=e[u−y(t)] are the equilibria of the system, i.e. the intersections of functions y and e(u−y). The stability of the equilibrium depends on the derivative of the vector field on the right-hand side of Eqn 9 with respect to y: if the derivative is negative, positive or zero, the equilibrium is stable, unstable or a saddle–node, respectively. The subplots along the top in Fig. 9 decompose the vector field at different values of u that produce qualitatively different nonlinear dynamics in the sense that the number and type of equilibria change. Assuming u changes sufficiently slowly, the dynamics should track the nearest stable equilibrium.

The number and type of equilibria depend explicitly on the structure of the escape function e(d) (where, recall, d=u−y). For the specific fit of e(d) in Fig. 7A we tested the predictions of the non-linear model. For this escape function the dynamics are punctuated by two bifurcation points at critical levels of the stimulus input, namely u≈–2.21 and u≈3.71. For constant stimuli below −2.21 there is one equilibrium (Fig. 9A). The single equilibrium (green dashed line) here is stable, as ẏ crosses from positive to negative at the equilibrium. At the critical level u≈−2.21, an additional equilibrium is introduced, namely a saddle–node (Fig. 9B). For −2.21<u<3.71, one stable and one unstable equilibrium branch out of the saddle–node, while the original stable equilibrium persists; thus, in this region there are three equilibria (Fig. 9C). But, at u=+3.71 the unstable branch intersects with the original stable branch creating a saddle–node (Fig. 9D), which vanishes for u>3.71, leaving just one lower stable equilibrium (Fig. 9E). The locus of equilibria at each u is shown in Fig. 9F (green curve); note that the middle branch of this inverted S-shaped curve is the unstable branch.

We simulated the dynamics for increasing and decreasing ramps between −30 and +30 Hz. We started the simulation with an initial condition of y(0)=0, which was near the only equilibrium for the initial values of u=±30 Hz. The red curve in Fig. 9 shows that, indeed, at each time t, y(t) remains near the closest stable equilibrium point until the output y(t) snaps through to the other branch. This occurs when the input u(t) reaches a saddle–node bifurcation. In the case of the increasing ramp from −30, both the simulated (red) and actual (gray, multiple trials) y(t) track the stable equilibrium on the upper branch of the green curve in Fig. 9F, even as the two additional equilibria are introduced at the bifurcation point of u≈−2.21. Eventually, the unstable equilibrium combines with the stable equilibrium that y(t) is tracking, causing both of them to vanish after forming a saddle–node. At that point, the closest (and only) stable equilibrium is the continuation of the lower branch, which causes y(t) to snap through, ultimately converging to the remaining equilibrium. Note that the sign of the df for a given stimulus frequency depends on whether the EOD frequency was initially above or below the stimulus frequency. The fact that the JAR circuitry is driven by the df (including its sign), not merely the stimulus frequency, manifests itself as the hysteresis loop shown in Fig. 9. Additionally, for the decreasing ramps, the snap-through for the actual data are delayed relative to the model. The overshoot of the model is indicative that either e(d) was underestimated in d>0 or there is a velocity dependence on d not captured by our model.

The JAR is a behavior in several species of weakly electric fish that allows an individual to shift its EOD frequency away from an interfering conspecific, whose frequency differs from its own by a low (but non-zero) value. The neural computation of the JAR is conceptually simple yet mechanistically complex: the fish achieve the JAR without an internal reference to their own EOD frequency (Heiligenberg et al., 1978; Bullock et al., 1972) and instead integrate information from receptive fields across the body surface, ultimately evaluating the single parameter df to raise or lower the EOD frequency so as to increase the magnitude of df (Heiligenberg, 1991).

Despite the mechanistic complexity of the JAR, our goal was to capture its conceptual simplicity and express it as a low-order dynamical system. To achieve this goal, we stabilized the JAR using a computer-controlled feedback system. This stabilization allowed us to reduce the complete computational algorithm of the JAR into a simple parsimonious model comprising a delay, a sensory escape function and a motor return (see Fig. 10).

Given the structure of the model described above, comparatively simple experimental measurements can now be used inform the parameters of a model that in turn predicts the responses to novel naturalistic or artificial stimuli. For example, this model captures a snap-through bifurcation that was not previously appreciated. Further, this model can be used to simulate social interactions between multiple individuals, without considering the mechanistic details of the internal dynamics underlying each individual's response.

A stimulus whose frequency perfectly matches that of the fish's own steady-state EOD would not elicit a JAR response; there would be no amplitude modulations or changes in the timings of zero crossings, as is necessary to drive the JAR circuit. This situation is an unstable equilibrium, because any deviations of either the stimulus or response frequency would cause the fish to shift its frequency away from it. The EOD frequency does not escape indefinitely, however, as it settles down into a new equilibrium that is a function of the applied stimulus frequency. The existence of additional stable equilibria reveals the inherent nonlinear nature of the behavior.

Further, it is known that there are two parallel motor pathways, one that shifts the EOD frequency up and the other down (Metzner, 1993). This implies that there could well be a significant non-linearities or even non-smoothness of the dynamics at the switching point between the two circuits, i.e. the baseline. Thus, it is imperative that we understand the local dynamics at baseline, in addition to capturing the nonlinear escape dynamics. Indeed, it is the global asymmetric shape of the nonlinear escape curve e(d) – and not a local discontinuity – that captures the asymmetry that was described previously by Metzner. The local linear model implies that the acceleration and deceleration of the pacemaker signal initiated by these two parallel pathways are similar in magnitude around baseline, causing a ‘smooth’ transition between the circuits. The mechanism of the JAR is based on spatial integration and voting, and not all units agree on the direction of the df (Heiligenberg et al., 1978). Thus, as the stimulus transitions from negative to positive df (or vice versa), the balance in the differential recruitment of the ‘up versus down’ JAR circuits tips, in a manner that is, evidently, graded and linear.

The JAR is an escape response, which makes the examination of the behavior around the baseline challenging. It is the instability of the equilibrium at baseline that underlies this difficulty. However, the stabilization of an unstable system using feedback is a classical application of control theory. In the case of the fish, we stabilized the unstable equilibrium by setting up a closed-loop feedback system that used the error signal between the fish's own frequency and a predefined reference frequency trajectory as the basis for the feedback. This system allows us to ‘dictate’ the frequency of the fish. We showed experimentally (Fig. 2) and confirmed analytically after modeling (Eqn 8) that this feedback did indeed stabilize the open-loop system.

This artificially closed-loop system was systematically perturbed (Fig. 3), and the resulting FRF data were numerically transformed to open loop using the known feedback transformation (Eqn 4, Fig. 4B). This allowed us to identify a model structure and verify it on the experimentally obtained closed-loop FRF data (Fig. 4A). In this manner, by knowing the feedback parameters and the closed-loop response, we were able to fit a frequency domain model to an unstable behavior. This procedure is relevant and applicable to a wide range of unstable biological behaviors, especially robust escape behaviors.

In each of the fish that we tested (five for reference tracking and two for clamp), the linearization at baseline was found to be identical. This was an unexpected result, as jamming avoidance is possible with different slopes for e(d) at the origin, as long as the slopes are negative. Assuming similar motor dynamics, this suggests that small excursions around baseline will be similar across individuals. Further, the linearizations were also found to be independent of stimulus amplitude. Previous work on the JAR (Bullock et al., 1972) has clearly shown that stimulus amplitude affects the ‘strength’ of the JAR in the sense that the magnitude of the frequency excursion is positively correlated with the depth of modulation. This would suggest that, in view of our model, the effect of stimulus amplitude would manifest as a change in the ‘height’ of the escape curve (e.g. see Fig. 7) but that the slope at the origin remains unchanged. The consistency of the slope at the origin also means that the escape curve, which can be identified by clamped-df experiments, serves as a possible ‘signature’ for the JAR behavior. By quantifying the escape curve for two individuals, jamming interactions between them can be modeled and predicted.

Reducing the JAR to this curve also allows us to examine social aspects of the JAR without considering the mechanistic details underlying its dynamics. For instance, previous research, interpreted in the context of our model, suggests that the escape curves may vary across development and be dependent on sex (Kramer, 1987). This means that the distribution of frequencies in a social group might be tailored by the individuals over time by modification of their escape curves. Future work might include identifying a general model structure that works across individuals and across stimulus amplitudes so that a complete interaction between freely swimming individuals in the wild can be modeled.

Finally, the snap-through response shown in Fig. 9 was first reported as part of the work that led to the discovery of the JAR (Watanabe and Takeda, 1963) (Fig. 5), and was described as the “effect of ‘chasing’ the discharge frequency by changing the stimulus frequency towards the discharge frequency”. This sudden reversal of the direction of df can now be understood as a bifurcation in the dynamics of the nonlinear model described in this paper.

The effect of stimulus amplitude on the dynamics has not been fully explored in this work. While we found no significant dependence on constant amplitude (see Fig. 6), the JAR certainly depends on changing amplitude as a changing amplitude is simply an AM. In a natural environment, the amplitude of EOD signals from nearby conspecifics will depend on the time-varying distance and orientation between the animals. If the goal is to model the interaction of groups of fish, it becomes necessary to quantify the effect of stimulus amplitude dynamics on the JAR.

In addition, the escape function may depend on higher order time derivatives of d, taking the form e(d,,,…). An asymmetry in this kind of velocity dependence may explain the difference between the model response and data in the case of the decreasing ramps. Experiments similar to the dynamic clamp trials described in this paper are required to fully understand the contribution of these higher order terms.

In their natural habitat, weakly electric fish are often found in groups of two or more individuals (Stamper et al., 2010; Tan et al., 2005). In these groups, the complex interaction of electric fields can give rise to modulations called envelopes (Stamper et al., 2013), through the relative movement between individuals (Yu et al., 2012), or the higher order interaction of their EODs (Stamper et al., 2012). Evidence of envelope processing has been revealed in the electrosensory system in weakly electric fish (Middleton et al., 2006; Longtin et al., 2008; McGillivray et al., 2012; Savard et al., 2011). In addition, Eigenmannia have a behavioral response to low-frequency ‘social envelopes’ (Stamper et al., 2012). Extending our JAR model not only to incorporate pair-wise differences between individuals but also to capture envelope responses would ultimately provide a powerful tool to predict and interpret complex social behavior in these fish.

This analysis can also be used to further explore the neural circuitry that controls the JAR. These behavioral dynamics are, of course, manifest in CNS circuits. Specifically, our model parses escape as a sensory-driven phenomenon. Based on this, we might expect to discover neural correlates of the escape component in electrosensory areas, such as the electrosensory lateral line lobe, midbrain torus semicircularis, and perhaps the nucleus electrosensorius. Indeed, it is well known that these areas encode and perform the necessary computations for the control of the JAR (Heiligenberg, 1991). In contrast, our model suggests that the return is mediated by motor processes. If so, then we would expect to discover neural correlates of the return in the properties of neurons in the pre-pacemaker nucleus or pacemaker nucleus. This is consistant with previous work in Apteronotus which examined how the baseline EOD frequency is controlled and modulated through activity in the pacemaker (Oestreich and Zakon, 2002). Indeed, our approach could be extended to examine longer term modulation of the motor system by introducing biases or offsets into the feedback controller.

The authors would like to thank Eatai Roth for insightful discussions on system identification, and also Jon Dyhr and Sean Carver for discussions on developing the model fitting scheme.

This section describes the modeling approach for the closed-loop reference tracking trials. Three of the four fish used in the reference frequency tracking trials were tested using long chirp stimuli whose reference frequency started at 0.001 Hz and increased to 1 Hz; the resulting FRF data binned into 30 frequency bands (see Materials and methods). Two chirp magnitudes (1 and 2 Hz) were tested for each individual for a total of six trials that were used for fitting.

The next step was to determine the model structure. A model-fitting criterion such as reduced χ², or information criteria such as Akaike information criterion (AIC) or Bayesian information criterion (BIC) can be generally minimized to determine the model order. However, these criteria assume independent data points. As the chirp is a time-varying frequency stimulus, the response data are clearly covariant, and the relatively small number of trials (six) is insufficient to determine the covariance structure between 30 bins. So, we implemented a selection criterion decision technique that takes into account two factors: model fit and model consistency.

Leave-one-out cross-validation was used as the technique to fit a model while penalizing over-fitting. This technique is asymptotically equivalent to AIC (Stone, 1977). For each model, the Matlab non-linear optimization function fminsearch was run with the error function (Eqn A4) to fit parameters to data from five trials, leaving the sixth trial out. This optimization was initialized 100 times (using Matlab to generate pseudo-random initial parameter values), and the set of parameters with the lowest fit error among them was then compared with the left-out trial, using the same error function. The leave-one-out error determined for one trial was averaged over all six trials being left out, providing a combined leave-one-out error for the model structure (see Fig. A1, y-axis).

Even though a certain model might have low leave-one-out error, the parameters fitted during each of the six leave-one-out minimizations may vary significantly. This variance, which is a measure of the uncertainty of the parameters based on our data, needs to be penalized, as we desire to have a consistent model fit with all six trials. To measure model consistency, we calculate the maximum singular value of the K×6 matrix where each of the six columns contains an estimate of all K parameters, with one column for each fit (one for each of the six trials left out). The poles and zeros of each fit are sorted, and the mean of each parameter across trials is subtracted, before finding the maximum singular value.

The higher order models have lower error in some instances, but are much less consistent. As shown in the inset of Fig. A1, among the two ‘good’ models, we chose the model closest to the origin, (0,1,1), corresponding to a delayed first-order system.

After determining the model structure, the parameters for each model were ultimately estimated without cross-validation (i.e. using all available data). Optimization was performed 100 times for each model structure with random initial parameter guesses. The lowest error fit among these was designated the fit error, the parameters associated with which were used. For the two ‘good’ models (see Table A1), we inspected the histogram of all 100 final parameters and found that most initial parameter guesses converged to the same final parameter values. This suggests that, given the wide set of initial guesses, our final estimate was not trapped in a local minimum. For clarity, we present the best parameters and errors of only the 10 lowest order models shown in Fig. A1 and Table A1. Fig. A2 compares experimental FRFs with the four lowest order models. The (0,0,0) and (0,0,1) models are obviously poor fits to the data, whereas the two-parameter model (0,1,0) fits the data well except at the highest frequencies. The (0,1,1) model compensates for this by introducing a phase lag due to the delay. There is no discernible improvement in the frequency responses of models of higher order than (0,1,1), so they are not shown. The best-fit parameters for the (0,1,1) model are reported in Results (see Eqn 7).

 
• AIC

  Akaike information criterion

 
• AM

  amplitude modulation

 
• BIC

  Bayesian information criterion

 
• D(s)

  Laplace transform of difference signal

 
• d(t)

  difference signal (Hz)

 
• df

  difference frequency (Hz)

 
• e(d)

  escape function

 
• EOD

  electric organ discharge

 
• f₁

  EOD frequency (Hz)

 
• f₂

  stimulus frequency (Hz)

 
• FRF

  frequency response function

 
• G(s)

  behavior transfer function

 
• H(s)

  closed-loop transfer function

 
• J(s)

  open-loop transfer function

 
• JAR

  Jamming Avoidance Response

 
• K

  number of parameters in candidate model

 
• k

  transfer function gain (Hz)

 
• M(s)

  model transfer function

 
• N

  number of fish

 
• n

  number of trials

 
• n_d

  existence of delay (0/1) in candidate model

 
• n_p

  number of poles in candidate model

 
• n_z

  number of zeros in candidate model

 
• p

  pole of transfer function (Hz)

 
• R(s)

  Laplace transform of reference signal

 
• r(t)

  reference signal (Hz)

 
• RR

  response–response coherence

 
• s

  Laplace complex frequency variable

 
• SR

  stimulus–response coherence

 
• T

  delay of transfer function (s)

 
• t

  time (s)

 
• U(s)

  Laplace transform of input signal

 
• u(t)

  input signal (Hz)

 
• Y(s)

  Laplace transform of output signal

 
• y(t)

  output signal (Hz)

 
• α

  feedback gain

 
• ρ(y)

  return function

 
• τ

  time constant (s)

 
• ω

  frequency (rad s⁻¹)