In high-dimensional stimulus spaces, the same response from a sensory neuron can be elicited by a continuum of equally effective optimal stimuli rather than a single optimal stimulus (**Figure 1**). Therefore, in some experiments it may be of interest to find sets of equivalent stimuli known as *iso-response surfaces* which yield the same response. One possible way of formalizing an optimization problem for this class of experiments is to formulate it as finding stimuli

$$\mathbf{x}^\* = \operatorname\*{arg\,min}\_{\mathbf{x} \in X} d(f(\mathbf{x}), \boldsymbol{\varepsilon}), \tag{4}$$

which *d*(· ,·) is some metric measure (e.g., squared error) quantifying the discrepancy between the desired response *c* and the neuronal response *f*(**x**). Multiple optimization runs from different starting locations and for different values of the desired constant response *c* permit the experimenter to determine families of iso-rate surfaces for the neuron under study. The geometrical shapes of the iso-rate surfaces can help to determine how stimulus variables *x*1,...,*xn* are integrated, and thus provide a useful tool for comparison of hypothetical models. For instance, linear integration of stimulus energy would yield iso-response surfaces which are hyperplanes of the form *n i*=1 *xi* = *c*, whereas non-linear integration would yield non-planar iso-response surfaces. **Figure 3** illustrates iso-response surfaces for two different hypothetical sensory processing models.

The iso-response surface method was used by Gollisch et al. (2002) to test several competing hypotheses about how spectral energy is integrated in locust auditory receptors. The iso-response contours to combinations of two or three pure tone stimuli with fixed frequencies and variable amplitudes were of elliptical shape, consistent with an energy-integrator model of spectral integration. Further work extended the iso-response method to incorporate temporal integration, yielding a complete cascade model of auditory transduction (Gollisch and Herz, 2005).

A more recent study applied this technique to study the integration of visual contrast over space in salamander retinal ganglion cells, revealing a threshold-quadratic non-linearity in the receptive field center as well as a subset of ganglion cells most sensitive to spatially homogeneous stimuli (Bölinger and Gollisch, 2012). The iso-response surface method has also been applied fruitfully in mammalian sensory systems as well. A recent study by Horwitz and Hass (2012) utilized this procedure to study integration of color signals from the three retinal cone types in single neurons in the primary visual cortex. It was found that half of the neurons had planar iso-response surfaces, consistent with linear integration of cone signals. However, the other half showed a variety of non-linear iso-response surfaces, including cup-shaped surfaces indicating sensitivity to only narrow regions of color space.

Although the iso-response surface method has been applied successfully in stimulus spaces of low dimensionality (two or three dimensions), tracing out level hyper-surfaces in higher-dimensional spaces may pose a formidable computational

challenge (Han et al., 2003; Willett and Nowak, 2007). In future research, dimensionality reduction procedures might be useful for extending the iso-response surface method to high-dimensional stimulus spaces like pixel space or auditory frequency space (Yu and Young, 2000), as well as for high-dimensional spaces defining complex naturalistic stimuli like 3D shapes or species-specific communication sounds (DiMattina and Wang, 2006; Yamane et al., 2008).

### **MAXIMALLY INFORMATIVE STIMULUS ENSEMBLES**

It has been proposed that one of the major goals of sensory coding is to efficiently represent the natural environment (Barlow, 1961; Simoncelli, 2003). In this spirit, another class of closed-loop stimulus optimization methods has been developed to find the optimal *ensemble* of sensory stimuli for maximizing the mutual information between stimuli and neural responses (Machens, 2002). This method differs from efforts to find the optimal stimulus or efforts to find iso-response surfaces because the goal is not to find an individual stimulus **x**∗ which optimizes the desired criterion (i.e., Eq. 1), but rather to find the optimal distribution *p*∗(**x**) which optimizes the mutual information *I*(y; **x**), where y denotes the observed neural response (typically the firing rate of a single neuron). Mathematically, we write

$$p^\*(\mathbf{x}) = \operatorname\*{arg\,max}\_{p(\mathbf{x}) \in P} I(\mathbf{y}; \mathbf{x}) = \int\_X \int\_Y p(\mathbf{y} \mid \mathbf{x}) p(\mathbf{x}) \ln \frac{p(\mathbf{y} \mid \mathbf{x})}{p(\mathbf{y})} d\mathbf{x} d\mathbf{y} \,, \tag{5}$$

where *P* is the space of probability densities on the stimulus space *X*, and *p*(*y* | **x**) and *p*(*y*) are determined experimentally by observing neural responses to stimuli. In practice, one starts with an assumption of a uniform distribution with finite support and then applies the Blahut–Arimoto algorithm (Blahut, 1972; Arimoto, 1972) to iteratively update the sampling weights (Machens, 2002). This method has been applied experimentally to characterize grasshopper auditory receptor neurons, demonstrating optimality for processing behaviorally relevant species-specific communication sounds (Machens et al., 2005; Benda et al., 2007).

## **ADAPTIVE OPTIMIZATION FOR SENSORY MODEL ESTIMATION**

An ideal gold standard for sensory neuroscience is to obtain a complete and accurate functional stimulus–response model of the neuron under study. In theory, once such a model is attained, one can then numerically or analytically calculate from this model the neuron's optimal stimulus, its iso-response surfaces, and its maximally informative stimulus ensembles. That is, if one identifies the system, one gets"for free" other information one may be interested in. However, despite its conceptual appeal, the problem of system identification is of great practical difficulty. This is because one needs to specify an accurate yet experimentally tractable model whose parameters can be estimated from data obtained during the available observation time. Unfortunately, research in computational neuroscience has shown that tractable (e.g., linear and quadratic) models are not accurate, whereas more biologically accurate models (deep, multiple layer neural networks incorporating dynamics, recurrence, etc.) often pose intractable parameter estimation problems.

"fncir-07-00101" — 2013/6/5 — 10:46 — page 6 — #6

It is well known from the fields of statistics and machine learning that one can more quickly and accurately estimate the parameters of a function using adaptive data collection, where new observations are generated in an iterative, adaptive manner which optimize the expected utility of the responses given the goal of estimating the model parameters (Lindley, 1956; Bernardo, 1979; MacKay, 1992). Mathematically, the optimization problem is to find at each iteration

$$\mathbf{x}\_{n+1}^{\*} = \underset{\mathbf{x} \in X}{\text{arg}\max} \ U\_{n}^{(\mathcal{E})}(\mathbf{x}),\tag{6}$$

where *U*(E) *n* (**x**) is the estimation utility function based on the data of the first *n* stimulus–response pairs. There are many choices for this function, including expected squared error (Müller and Parmigiani, 1995), expected prediction error (Sugiyama, 2006), and mutual information between stimuli and model parameters (Paninski, 2005). The generic name for this strategy is *optimal experimental design* or OED (Federov, 1972; Atkinson and Donev, 1992; Cohn et al., 1996), and it is often studied in a Bayesian framework (MacKay, 1992; Chaloner and Verdinelli, 1995). Recent theoretical and experimental work has shown that such methods can potentially be fruitfully applied in neuroscientific experiments (Paninski, 2005; Paninski et al., 2007; Lewi et al., 2009, 2011; DiMattina and Zhang, 2011). Not only can optimal experimental design make the estimation of high-dimensional models practical (Lewi et al., 2009), but can also make it tractable to estimate highly non-linear models which cannot be readily identified from random "white noise" data of the kind typically used in system identification experiments (DiMattina and Zhang, 2010, 2011). We first discuss the application of such methods in psychology and cognitive science, and then discuss recent theoretical and experimental work on applications of OED methods to sensory neurophysiology experiments

#### **ADAPTIVE STIMULUS OPTIMIZATION IN PSYCHOLOGY AND COGNITIVE SCIENCE**

Psychophysics has long utilized adaptive data collection, with the classic example being the staircase method for threshold estimation (Cornsweet, 1962). More recently, an adaptive Bayesian approach to threshold estimation (QUEST) which chooses new stimuli for each trial at the current Bayesian estimate of the threshold was developed (Watson and Pelli, 1983), and subsequent work extended this approach to permit simultaneous estimation of both the threshold and slope of the psychometric function (Snoeren and Puts, 1997). Another line of work applied an information-theoretic approach to estimating the slope and threshold parameters, where stimuli were chosen at each trial to maximize the expected information gained about the slope and threshold parameters (Kontsevich and Tyler, 1999). More sophisticated methods of this kind have been utilized for psychometric

"fncir-07-00101" — 2013/6/5 — 10:46 — page 7 — #7

functions defined on two-dimensional stimuli (Kujala and Lukka, 2006), with these procedures being applied for estimating contrast sensitivity functions (Lesmes et al., 2010) and color sensitivity of human observers (Kujala and Lukka, 2006). In addition to finding widespread application in sensory psychophysics, adaptive methods have also been used more broadly in the cognitive sciences (Wixted and Ebbesen, 1991; Rubin and Wenzel, 1996; Nosofsky and Zaki, 2002; Opfer and Siegler, 2007; Kujala et al., 2010; McCullough et al., 2010).

#### **GENERALIZED LINEAR MODELS AND BIOPHYSICAL MODELS**

More recently, investigators in computational neuroscience have demonstrated that adaptive information-theoretic sampling where stimuli are chosen to maximize the expected information gain between a stimulus and the model parameters can be a highly effective means of estimating the parameters of sensory processing models (Paninski, 2005; Paninski et al., 2007). A fast information-theoretic algorithm has been developed for the generalized linear model which applies a static non-linearity to the output (Lewi et al., 2009). The generalized linear model has been utilized in numerous studies (Simoncelli et al., 2004) and enjoys a likelihood function with no local maxima (Paninski, 2004). Their algorithm relied on a Gaussian approximation to the posterior density, permitting fast recursive updates, with the calculations for finding the optimal stimulus growing only as the square of the stimulus space dimensionality. Numerical simulations demonstrated that their procedure was asymptotically efficient, with the empirically computed variance of the posterior density converging to the minimum theoretically possible variance.

One issue which potentially affects studies of stimulus optimization is neuronal adaptation due to the stimulus history (Ulanovsky et al., 2003; Asari and Zador, 2009). In sensory neurons, this may be manifested as the system actually changing its underlying parameters which we seek to estimate as the experiment progresses. However, the procedure developed by Lewi et al. (2009)was demonstrated to be robust to parameter drift in numerical simulations, suggesting the ability to compensate for changes to the system brought about by adaptation effects. Furthermore, their model also permits the estimation of a spike-history filter, allowing neuronal response history to influence predictions to new stimuli.

A further study by this group applied this algorithm to fitting generalized linear models to avian auditory neurons probed with conspecific song samples, and it was found that accurate estimation could be obtained using vastly fewer samples when they were chosen adaptively using the algorithm then when they were chosen non-adaptively (Lewi et al., 2011). Although this procedure has yet to be applied in real on-line experiments, it provides experimenters working on a variety of systems with a powerful tool for quickly characterizing neurons whose responses are well described by generalized linear models (Chichilnisky, 2001) or related models (Pillow et al., 2008).

More recently, this group has also applied optimal experimental design to the cellular neuroscience problem of accurately estimating voltages from dendritic trees using measurements suffering from low signal-to-noise ratio (Huggins and Paninski, 2012). Using simulated compartmental models, these authors demonstrated that by adaptively choosing observation locations which minimize the expected squared error of the voltage measurement, a substantial improvement in accuracy was obtained compared to random sampling. This procedure is potentially of great experimental usefulness because techniques like two-photon imaging permit spatially complete observations of dendrites, but with low signal-to-noise ratios (Djurisic et al., 2008; Canepari et al., 2010).

## **MULTIPLE LAYER NEURAL NETWORKS**

Since many sensory neurons are non-linear (Young et al., 2005;Wu et al., 2006), it is of interest to characterize neurons using various non-linear models, including quadratic and bilinear models (Yu and Young, 2000; Berkes andWiskott, 2006;Ahrens et al., 2008a,b), neural networks (Lau et al., 2002; Prenger et al., 2004; Cadieu et al., 2007) and basis function networks (Poggio and Girosi, 1990). A generalized linear model is also a non-linear model because it employs a static non-linearity at the output stage. Although a generalized linear model allows limited non-linearities, it enjoys tractable and consistent estimation procedures without problems of local minima (Paninski, 2004). Identifying more complex nonlinear models like hierarchical neural networks from physiological data tends to be harder due to problems like local minima and plateaus in the error surface (Amari et al., 2006; Cousseau et al., 2008; Wei and Amari, 2008; Wei et al., 2008).

For studies aimed at estimating generalized linear models, the use of a fixed white-noise stimulus set is often quite effective and is theoretically well-justified (Chichilnisky, 2001; Paninski, 2004;Wu et al., 2006). However, recent theoretical work suggests that using fixed stimulus sets like white noise may be deeply problematic for efforts to identify non-linear hierarchical network models due to *continuous parameter confounding* (DiMattina and Zhang, 2010). This problem is illustrated for a very simple non-linear neural network model shown in **Figure 4A**. In this example, the goal is to recover the parameters (*w,v*) of the network by performing maximum likelihood (ML) estimation given noisy stimulus–response observations. When the input stimuli *x* only drive the hidden unit over a region of its gain function which can be well approximated by a power function (**Figure 4B**, top), the estimates obtained by ML for different datasets lie scattered along the continuum *vw*α = *C*, as one would expect for a power law gain function *g*(*u*) = *Au*α (**Figure 4C**, top). (Here the constant *C* = *v*T*w*α T, where *w*T and *v*T are the true values of the input and output weights.) In contrast, when the input stimuli *x* drive the hidden unit over a full range of its gain so that the power law approximation is poor (**Figure 4B**, bottom), the true parameters are accurately recovered for different datasets (**Figure 4C**, bottom).

A hypothetical experiment which suffers from this problem is illustrated in **Figure 5**. We see that when the network in **Figure 5A** is probed with random stimuli (**Figure 5B**, right), the hidden unit is driven over a limited range of its gain function which may be well approximated by an exponential, so that the sigmoidal gain (**Figure 5C**, black curve) may *de facto* be replaced by a new exponential gain function *g*(*u*) = *Ae*α*u* (**Figure 5C**, red curve). With this new gain, it follows that a continuum of different values of the output weight *v* and bias *w*0 lying on the curve *ve*α*w*0 = *C* will yield models whose responses to the training data are

"fncir-07-00101" — 2013/6/5 — 10:46 — page 8 — #8

output weight parameters (*w*,*v*) we wish to estimate from noisy stimulus–response data. Noise is drawn from a Poisson distribution. **(B)** *Top*: The input stimuli *x* ∈ [−0.5, 0.5] only drive the hidden unit over a limited region of its gain function (black curve) which may be well approximated by a stimuli like that in the top of **Figure 4B**, the estimates (black dots) lie scattered along the curve predicted by the power law confounding theory. *Bottom*: When trained with sets of stimuli like those in the bottom panel of **Figure 4B**, the true parameter values (red triangle) are more reliably recovered.

indistinguishable and therefore multiple estimates of these parameters from different random training sets will lie scattered along this curve (**Figure 5D**). (Here the constant *C* = *v*T*e*α*w*0T where *V*T and *w*0T are the true values of the output weight and hidden unit bias.)

Adaptive stimulus optimization methods like informationtheoretic sampling (Paninski, 2005) can in principle overcome this problem of continuous parameter confounding, as we see in **Figure 5D** where the correct network parameters are reliably recovered when optimally designed stimuli (**Figure 5B**, left) are used. This simple example suggests that adaptive stimulus optimization may make it tractable to reliably recover the parameters of complex hierarchical networks needed to model non-linear neurons, whereas it is much harder to recover these networks using standard stimulus sets like white noise.

Many previous studies in the statistics and machine learning literature have demonstrated that faster convergence and smaller generalization error may be obtained when neural networks are trained adaptively using optimally designed stimuli (Lindley, 1956; MacKay, 1992; Cohn et al., 1994). Recently, we have developed a practical method for implementing the information-theoretic stimulus optimization approach derived for generalized linear models (Lewi et al., 2009) for arbitrary nonlinear models like hierarchical neural networks. Although this method employs numerous approximations, it has been shown in simulated experiments to be effective at recovering non-linear neural networks having multiple hidden units, and is fast enough to utilize in real experiments (Tam et al., 2011; Dekel, 2012; Tam, 2012).

## **ADAPTIVE OPTIMIZATION FOR SENSORY MODEL COMPARISON**

Quite often the appropriate model for describing a sensory neuron or perceptual quantity is unknown. Therefore, an important experimental goal may be to discriminate between two or more competing models. Mathematically, the optimization problem is to iteratively find stimuli

$$\mathbf{x}\_{n+1}^{\*} = \underset{\mathbf{x}}{\text{arg}\max} \ U\_{n}^{(\text{C})}(\mathbf{x}),\tag{7}$$

which optimize a model comparison utility function *Un* (C) (**x**), one choice of which may be the expected reduction in model space entropy (Cavagnaro et al., 2010; DiMattina and Zhang, 2011). This equation may be regarded as the optimal comparison counterpart of the equation for optimal estimation (Eq. 6). We now briefly discuss recent studies making use of adaptive stimulus optimization for model selection.

#### **PSYCHOPHYSICAL MODEL COMPARISON**

Although standard model comparison methods like the Bayesian Information Criterion (BIC; Schwarz, 1978) or predictive

"fncir-07-00101" — 2013/6/5 — 10:46 — page 9 — #9

input from a broadly integrating interneuron (I-unit). **(B)** Examples of optimally designed (left) and random (right) stimuli. Note that the optimally designed stimuli exhibit complex correlated structure. **(C)** Random stimuli (green dots) only drive the E-unit over a limited range of its gain function (black curve)

cross-validation may be applied *post hoc* (Vladusich et al., 2006; Wu et al., 2006), numerous studies suggest that performing experiments using stimuli optimized for model comparison may be far more effective (Atkinson and Fedorov, 1975a,b). One method for model comparison developed recently for psychophysical experiments is known as MAximum Differentiation (MAD) competition (Wang and Simoncelli, 2008). Given two perceptual models which relate stimulus parameters to a perceptual quantity, this method generates a pair of stimuli which maximizes/minimizes the response of one model while holding the other model's response fixed. Next, this procedure is repeated with the role of the two models reversed. Testing human subjects on the two pairs of synthesized stimuli can determine which model is "better" in the sense of telling us which model's max/min pairs are simpler to

the exponential theory (black curve). In contrast, estimates attained from optimally designed stimuli accurately recover the true parameters (red

triangle).

"fncir-07-00101" — 2013/6/5 — 10:46 — page 10 — #10

discriminate. This procedure has been fruitfully applied to comparing image quality assessment models which aim to predict human perception of image quality (Wang and Bovik, 2006)

An information-theoretic method for model comparison was recently derived by Cavagnaro et al. (2010). Given a set of models with the *i*-th model having prior probability *P*0(*i*), stimuli are chosen to maximize the mutual information between the stimulus and the model index *i* by minimizing the expected model space entropy in a manner directly analogous to informationtheoretic model estimation (Paninski, 2005), except that in this case the unknown variable is a discrete model index *i* rather than a continuous parameter value **θ**. This method was applied to competing models of memory retention from the cognitive science literature (Rubin et al., 1999) and was shown to permit much more accurate discrimination than standard non-adaptive methods.

#### **NEURAL MODEL COMPARISON**

In general, the correct parameters of competing sensory processing models are not known beforehand. Therefore, it is of interest to consider how to conduct experiments which estimate and discriminate competing models. Typically, investigators in neurophysiology and neuroimaging have applied model-comparison techniques *post hoc* (David and Gallant, 2005; Vladusich et al., 2006; Penny, 2012), particularly in the system identification literature (Prenger et al., 2004; David and Gallant, 2005; Wu et al., 2006; Sharpee et al., 2008; Rabinowitz et al., 2012; Schinkel-Bielefeld et al., 2012). However, a fundamental limitation with *post hoc* analysis is that it is not possible to generate and test critical stimuli which are optimized for model comparison, as this is only possible while the system is under observation. This limitation can only be overcome by fitting multiple models to a sensory neuron during the course of an experiment and then using the fitted models to generate and present critical stimuli which are optimized to best discriminate the models. Although previous work has presented stimuli on-line to test or verify a single model (deCharms et al., 1998; Touryan et al., 2002), very little work in single-unit *in vivo* sensory neurophysiology has presented stimuli optimized for model comparison in real-time (Tam et al., 2011).

A recent study considered a two-stage approach for combining the goals of model estimation and comparison in neurophysiology experiments, illustrated schematically in **Figure 6A** (DiMattina and Zhang, 2011). In the first stage, stimuli are adaptively optimized for parameter estimation, with the optimal stimulus for each model being presented in turn. In the second stage, stimuli are generated adaptively in order to optimally discriminate competing models making use of an information-theoretic criterion (Cavagnaro et al., 2010) or a likelihood-based criterion. In the special case of two models *f*1 (**x**, **θ**1), *f*2 (**x, θ**2) and Gaussian response noise, it can be shown that under a likelihood criterion the best stimulus for model discrimination is the stimulus which maximizes the quantity *f*1 (**x, θ**1) − *f*2 (**x, θ**2) 2 , and furthermore this stimulus will maximally increase the BIC in favor of the best model (DiMattina and Zhang, 2011).

**Figure 6** illustrates a numerical experiment making use of this two-stage procedure for the problem of discriminating an additive and multiplicative model of neural responses (**Figure 6B**), where the additive model is assumed to be the true model. After the estimation phase, the BIC does not have a strong preference for either model, only being correct about half the time (**Figure 6C**). However, after presenting 500 stimuli optimized for discriminating the additive and multiplicative model and applying the BIC to all available data, the correct (additive) model is preferred for 24 of 25 Monte Carlo trials (red curve). As a control, presenting additional stimuli optimized for model estimation only improves final model selection moderately (blue curve), while presenting random stimuli does not at all improve model selection performance (green curve). This procedure has now been applied in neurophysiology experiments to generate critical stimuli to distinguish between two competing models of spectral processing by single neurons in the primate inferior colliculus (Tam et al., 2011; Tam, 2012).

## **DISCUSSION**

With increasing computer power, it is becoming practical for neuroscience experiments to utilize adaptive stimulus optimization where stimuli are generated in real-time during the course of the experiment (Benda et al., 2007; Newman et al., 2013). Various experiments have utilized adaptive stimulus optimization in order to break the "curse of dimensionality" and find the optimal stimulus for a sensory neuron in spaces which are too large for factorial exploration (O'Connor et al., 2005; Yamane et al., 2008). However, simply characterizing the optimal stimulus for a sensory neuron provides at best only a partial description of neural coding (Olshausen and Field, 2005). Therefore, in addition to helping to find the optimal stimulus, adaptive stimulus optimization makes it possible to pursue engineering-inspired approaches to sensory neurophysiology which may yield greater functional insights, for instance finding stimulus ensembles maximizing information between stimuli (Machens, 2002; Machens et al., 2005) and neural responses or fitting and comparing multiple non-linear models to neural responses (Lewi et al., 2009; DiMattina and Zhang,2011). **Table 1** summarizes the various closed-loop stimulus optimization paradigms discussed in this review, and **Figure 7** schematically illustrates the closed-loop experimental approach.

The vast majority of the work to date has applied closedloop methods to studying scalar firing rate responses measured from single neurons. However, as closed-loop approaches are continuing to develop, and as new techniques like optical imaging (Ohki et al., 2005; Bock et al., 2011) are making it increasingly feasible to observe large numbers of neurons simultaneously, it is of great interest for future investigations to apply these methods to neural populations and to measurements beyond scalar firing rate. Here we briefly discuss some possible directions for future research.

While the notion of the optimal stimulus is well-defined for single neurons, it is not well-defined for neural populations. However, an alternative approach to stimulus optimization for a population of neurons is to find the stimulus at which the population is best at discriminating nearby stimuli, as opposed to the stimulus yielding the highest firing rate response. Indeed, it has been suggested by a number of investigators that highslope regions of tuning curves, where nearby stimuli are best

"fncir-07-00101" — 2013/6/5 — 10:46 — page 11 — #11

discriminated, are much more important in sensory coding than tuning curve peaks (Seung and Sompolinsky, 1993; Harper and McAlpine, 2004; Butts and Goldman, 2006; Bonnasse-Gahot and Nadal, 2008). Under reasonable assumptions of independent Poisson responses, the one-dimensional stimulus *x* at which a neural population can best discriminate nearby stimuli *x* + δ*x* is the stimulus which maximizes the Fisher information *I*F (*x*) = -*N i* = 1 *f i* (*x*) 2 / *fi* (*x*), where *fi*(*x*) is the tuning curve of the *i*-th neuron (Dayan et al., 2001). It is relatively straightforward to extend this Fisher information formalism to higher dimensional stimulus spaces (Zhang and Sejnowski, 1999; Johnson et al., 2001; Bethge et al., 2002). Local approximation of the Fisher information matrix has been used in previous work aimed at stimulus optimization in a single neuron (Bandyopadhyay et al., 2007b), and this technique could readily generalize to find the stimulus which is best discriminated from nearby stimuli by a population code.

Extension of the definition of iso-response surfaces (Gollisch et al., 2002) to multiple neurons is relatively straightforward. In particular, if we can view each neuron as implementing a function *f*(**x**) on the stimulus space, then the region of stimulus space which simultaneously satisfies multiple constraints *f*1 (**x**) = *c*1, ··· , *fN* (**x**) = *cN* should simply be the (possibly empty) intersection of the regions of stimulus space satisfying each individual constraint. It would be interesting to extend the maximally informative ensemble approach (Machens, 2002) to multiple neurons as well. One potential difficulty is that the number of possible responses which one needs to measure to compute

**Table 1 | Summary of various closed-loop stimulus optimization approaches utilized in sensory systems neuroscience.**

"fncir-07-00101" — 2013/6/5 — 10:46 — page 12 — #12

the probability distribution *p*(y|**x**) increases exponentially with the number of neurons in the population. Indeed, this exponential increase in the number of symbols with the dimensionality of the response space is a well-known problem with applications of information-theoretic methods in neuroscience (Rieke et al., 1997). It would be desirable to develop more efficient computational techniques for studying neuronal populations in the future (Yarrow et al., 2012).

In addition to considering neural populations, another direction for extending the closed-loop paradigm is to consider neural responses more sophisticated than firing rates, for instance the temporal patterns of neural responses (Optican and Richmond, 1987; Victor and Purpura, 1996), first spike latency (VanRullen et al., 2005; Gollisch and Meister, 2008), or synchronous responses in neural populations (Brette, 2012). Since a temporal pattern is a vector but not a scalar, one needs to extract a scalar quantity from a temporal pattern in order to define the optimal stimulus. For example, synchrony can be defined as a scalar quantity (Steinmetz et al., 2000) and can in principle be optimized over a stimulus space in the same manner as firing rate. The iso-response paradigm (Gollisch et al., 2002) would generalize quite well to both spike pattern and synchrony measures. In this case of spike pattern, the goal would be to find the equivalence class of all stimuli which could elicit a desired pattern of spiking, and theoretical efforts have demonstrated that it is possible to design stimuli to produce a desired spike pattern (Ahmadian et al., 2011). Similarly, for iso-synchrony curves one could find equivalence classes of stimuli yielding the same degree of synchrony in the population by utilizing algorithms similar to those developed for firing rate.

One of the most powerful applications of the closed-loop paradigm is the ability to move sensory neurophysiology toward a model-based paradigm, where experiments are performed with the goal of identifying and comparing multiple competing nonlinear models (Paninski, 2005; Lewi et al., 2009; DiMattina and Zhang, 2011; Tam et al., 2011). One advantage of model identification is that successful identification gives the experimenter a variety of biologically important information about the neuron or neuronal population "for free." That is, once one has determined an accurate model for a sensory neuron, the optimal stimulus for maximizing firing rate, the iso-response surfaces, or the stimulus ensemble maximizing information transmission can be predicted from this model, and these predictions can be tested experimentally. However, the model-based approach is not without its difficulties, as many sensory neurons are poorly described by tractable linear and quadratic models and may be better described by more complex models like basis functions and neural networks. Recent work has demonstrated that in principle, adaptive stimulus optimization methods long utilized in machine learning and psychophysics can be applied in sensory neurophysiologyfor purposes of model estimation and comparison (Paninski, 2005; Lewi et al., 2009; DiMattina and Zhang, 2011). In particular, our recent study has presented a practical two-stage experimental

"fncir-07-00101" — 2013/6/5 — 10:46 — page 13 — #13

method for generating stimuli which are optimal for estimating the parameters of multiple non-linear models and then generating stimuli on-line in order to critically compare the predictions of different models (DiMattina and Zhang, 2011). This method is presently being applied in ongoing auditory neurophysiology studies (Tam et al., 2011; Dekel, 2012; Tam,

#### **REFERENCES**

Bandyopadhyay, S., Young, E. D., and Reiss, L. A. J. (2007b). "Spectral edges as optimal stimuli for the dorsal cochlear nucleus," in *Hearing – From Basic Research to Applications*, eds B. Kollmeier, G. Klump, V. Hohmann, U. Langemann, M. Mauermann, S. Upperkamp, et al. (New York: Springer-Verlag), 39–45.

2012), and may be applicable to a broad variety of sensory systems.

#### **ACKNOWLEDGMENTS**

Supported by grant NSF IIS-0827695. Thanks to Sarah Osmer DiMattina for her assistance with graphic design.

"fncir-07-00101" — 2013/6/5 — 10:46 — page 14 — #14

determining visual receptive fields. *Vis. Res.* 14, 1475–1482.

response strategies, selective attention, and stimulus generalization. *J. Exp. Psychol. Learn. Mem. Cogn.* 28, 924–940.

"fncir-07-00101" — 2013/6/5 — 10:46 — page 15 — #15

Sharpee, T., Rust, N., and Bialek, W. (2004). Analyzing neural responses to natural signals: maximally informative dimensions. *Neural Comput.* 16, 223–250.

Shepherd, G. (2003). *The Synaptic Organization of the Brain*, 5 Edn. New York: Oxford University Press.

Simoncelli, E. (2003). Vision and the statistics of the visual environment. *Curr. Opin. Neurobiol.* 13, 144–149.

cat and monkey visual cortex. *Vis. Res.* 23, 775–85.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 30 August 2012; accepted: 08 May 2013; published online: 06 June 2013.*

*Citation: DiMattina C and Zhang K (2013) Adaptive stimulus optimization for sensory systems neuroscience. Front. Neural Circuits 7:101. doi: 10.3389/ fncir.2013.00101*

*Copyright © 2013 DiMattina and Zhang. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in other forums, provided the original authors and source are credited and subject to any copyright notices concerning any thirdparty graphics etc.*

"fncir-07-00101" — 2013/6/5 — 10:46 — page 16 — #16

## A Hebbian learning rule gives rise to mirror neurons and links them to control theoretic inverse models

## *A. Hanuschkin1,2, S. Ganguli 3 and R. H. R. Hahnloser 1,2\**

*1 Institute of Neuroinformatics, University of Zurich and ETH Zurich, Zurich, Switzerland*

*3 Department of Applied Physics, Stanford University, Stanford, USA*

#### *Edited by:*

*Ahmed El Hady, Max Planck Institute for Dynamics and Self Organization, Germany*

#### *Reviewed by:*

*Yang Dan, University of California, USA Klaus R. Pawelzik, University*

#### *Bremen, Germany*

*\*Correspondence: R. H. R. Hahnloser, Institute of Neuroinformatics, University of Zurich and ETH Zurich, Winterthurerstrasse 190, CH-8057 Zurich, Switzerland e-mail<EMAIL_ADDRESS>
Mirror neurons are neurons whose responses to the observation of a motor act resemble responses measured during production of that act. Computationally, mirror neurons have been viewed as evidence for the existence of internal inverse models. Such models, rooted within control theory, map-desired sensory targets onto the motor commands required to generate those targets. To jointly explore both the formation of mirrored responses and their functional contribution to inverse models, we develop a correlation-based theory of interactions between a sensory and a motor area. We show that a simple eligibility-weighted Hebbian learning rule, operating within a sensorimotor loop during motor explorations and stabilized by heterosynaptic competition, naturally gives rise to mirror neurons as well as control theoretic inverse models encoded in the synaptic weights from sensory to motor neurons. Crucially, we find that the correlational structure or stereotypy of the neural code underlying motor explorations determines the nature of the learned inverse model: random motor codes lead to causal inverses that map sensory activity patterns to their motor causes; such inverses are maximally useful, by allowing the imitation of arbitrary sensory target sequences. By contrast, stereotyped motor codes lead to less useful predictive inverses that map sensory activity to future motor actions. Our theory generalizes previous work on inverse models by showing that such models can be learned in a simple Hebbian framework without the need for error signals or backpropagation, and it makes new conceptual connections between the causal nature of inverse models, the statistical structure of motor variability, and the time-lag between sensory and motor responses of mirror neurons. Applied to bird song learning, our theory can account for puzzling aspects of the song system, including necessity of sensorimotor gating and selectivity of auditory responses to bird's own song (BOS) stimuli.

**Keywords: mirror neurons, inverse problem, linear models, songbird, sensory motor learning**

## **INTRODUCTION**

Complex vertebrate motor behaviors are generated by dedicated cortical circuits. The organization of these circuits and the plasticity rules that lead to their development and that guarantee their maintenance are functionally related to neural activity in single units and across larger populations (Gallese et al., 1996; Rizzolatti et al., 1996; Rizzolatti and Craighero, 2004; Harvey et al., 2012). For example, neural activity often strongly co-varies with motor behavior, allowing for estimation of detailed limb movement parameters from mere single-neuron recordings (Georgopoulos et al., 1986; Schwartz et al., 1988) and facilitating neural prosthesis (Santhanam et al., 2006; Ethier et al., 2012). However, in other cases, the amount of firing variability in single neurons can be dramatically dissociated from behavioral variability. For example, in songbirds, two distinct premotor areas are responsible for the generation of different aspects of the same vocal behavior. On the one hand, the cortical area HVC is involved in generating stereotyped adult song; lesions of HVC lead to degradation of typical adult song toward more unstructured subsong typical of very young birds (Nottebohm et al., 1976; Aronov et al., 2008). On the other hand, its counterpart, the lateral magnocellular nucleus of the anterior nidopallium (LMAN) in very young birds is involved in subsong production and in adults it is involved in the production of very subtle song variability that is barely noticeable to the human ear (Aronov et al., 2008). Lesions of LMAN in juveniles abolish song learning (Bottjer et al., 1984), and lesions in adults reduce the already small variability of adult undirected songs (the songs not direct toward another bird), manifest for example by reduced fluctuations of sound pitch (Kao et al., 2005; Stepanek and Doupe, 2010).

These lesion studies ascribing differential roles of HVC and LMAN to song production, are paralleled by findings from electrophysiology. In HVC of singing birds, single principal neurons fire highly stereotyped spiking patterns associated with a given song syllable, with precision of individual action potentials in the sub millisecond range (Hahnloser et al., 2002; Kozhevnikov and Fee, 2007). By contrast, in LMAN of birds singing undirected songs, neurons fire very variable spike patterns, patterns

*2 Neuroscience Center Zurich (ZNZ), Zurich, Switzerland*

that fluctuate on a trial-to-trial basis between loosely timed high-frequency bursts of action potentials and no spiking at all (Olveczky et al., 2005; Kao et al., 2008). Thus, stereotyped adult song is subserved by precise firing in HVC whereas subtle variability of adult song is subserved by large firing variability in LMAN, **Figure 1**. The diverse neural codes in LMAN and HVC are integrated in a dedicated nucleus that mediates both differential influences from these stereotypy and variability generators. Both HVC and LMAN project to the robust nucleus of the arcopallium (RA), which is the cortical output nucleus that directly innervates syringeal and respiratory motor neurons.

Whether stereotyped or variable, internal motor patterns responsible for generating behavior cannot be fully understood without considering the sensory input reaching the motor system.

stereotyped song motif (in a different bird), three exemplary motif oscillograms are shown on top. In each rendition of the motif the neuron

produces a brief burst of spikes at precisely the same time.

Indeed, the very development of motor systems as well as the formation of motor plans are profoundly shaped by sensory inputs. For example, the development of the mirror neuron system depends on sensorimotor experience (Catmur, 2012) and, the successful development of birdsong depends on intact HVC and LMAN activity during sensory exposure (Basham et al., 1996; Roberts et al., 2012).

We have learned much about the integration of sensory inputs into motor systems from single neuron studies examining responses during motor production and during matched sensory states.   Mirror neurons also exist in HVC of songbirds; these neurons fire at a precise time in the song, both when the bird sings the song and when it hears a similar song produced by another bird (Prather et al., 2008).

Mirror neurons establish a link between the observation of an act in another and self-generation of that same act. Such a remarkable correspondence between sensory and motor roles in single neurons has led to numerous suggestions about the function of mirror neurons in communication, imitation learning, cultural learning, and language development (Rizzolatti and Craighero, 2004; Oztop et al., 2012). Most importantly, mirrored responses have been proposed to be causally related to streams of motor and sensory activity (Oztop et al., 2006, 2012). A recent proposal is to tie properties of the mirror neuron system to correlative learning rules (Cooper et al., 2012). Accordingly, sensory responses in mirror neurons could develop from the contingency of motor-related firing and its sensory consequences feeding back to motor areas. Here we develop this idea and propose a simple mathematical theory of mirror neuron formation from correlational learning rules. To examine the critical role of motor variability, we study, based on earlier work (Hahnloser and Ganguli, 2013), mirror neuron formation for both motor codes with strongly correlated firing patterns among neurons, as in HVC, as well as for motor codes with uncorrelated firing patterns among neurons, as in LMAN.

We are particularly interested in relating mirror neuron properties to their computational role in control theoretic inverse models. Mirror neurons have previously been recognized as direct evidence of inverse models, which are models that transform desired sensory states into motor commands that can achieve those states and may be used for action generation (Oztop et al., 2012). From the control-theoretic perspective, internal inverse models give rise to mirrored responses because of the precise correspondence between a desired sensory target, the motor commands for producing that target, and the resulting sensory feedback. We pursue this idea and elucidate the conditions under which inverse models can arise from correlational learning during sensory feedback-dependent motor explorations.

We assume inverse models form in a context without prior knowledge of structure of either the motor apparatus or the delayed sensory feedback. We design an eligibility-weighted correlational learning rule that allows for the formation of both inverse models and mirror neurons. In the rule we propose, synaptic strengthening depends on contiguous co-activation of pre-and postsynaptic neurons, whereas synaptic weakening depends on heterosynaptic competition between sensory afferents innervating the same motor neuron. We argue that from a synaptic perspective, this rule is considerably simpler and more plausible than previously proposed rules and computational approaches toward systems-level inverse models based on error backpropagation (Jordan and Rumelhart, 1992). Our rule is most closely related to direct inverse model approaches (Miller, 1987; Slotine, 1987), in which, however, the possibility of unknown feedback delays has not been adequately addressed. Most importantly, we find that whether the formed mirror neuron system and inverse model is suitable for action imitation depends on the correlational structure of the neural code associated with motor production. Whereas a variable (explorative) motor code leads to causal inverse models and is suitable for mirror-neuron dependent action imitation, a stereotyped (repetitive) motor code leads to predictive inverse models and is not suitable for action imitation. Thus, our work provides an interesting link between the correlational structure of motor behavior, its underlying neural code, and fine-grained temporal properties of mirror neuron responses and their suitability for flexible action imitation.

Furthermore, these conceptual connections suggest a set of natural experiments designed to probe for the existence, and characterize the causal nature of, inverse models by measuring the fine grained temporal properties of the sensory and motor responses of mirror neurons. As we discuss below, when applied to the bird song system, these experiments make a specific, testable prediction about the existence and temporal properties of mirror neurons in the variable motor circuit LMAN, as well as explain the origin of previously observed temporal properties of mirror neurons in the stereotyped motor circuit HVC.

## **RESULTS**

#### **A LINEAR FRAMEWORK**

We develop our theory in a simple linear framework in which the sensory response *a(t)* in a sensory brain area at time *t* is a vector of firing rates that is linearly related to the motor cause *m(t* − *τ)* at an earlier time *t* − *τ* , where *m(t* − *τ)* is a vector of firing rates in a motor area such as HVC or LMAN. The *time delay of sensory feedback τ* = *τm* + *τa* is the sum of the time *τm* needed to translate motor activity into behavioral (vocal) output and the time *τa* it takes for a vocalization to elicit a sensory response. We assume a linear *motor-sensory mapping* modeled by the matrix **Q**, allowing us to specify the form of delayed sensory feedback as *a(t)* = **Q***m(t* − *τ)*, **Figure 2**.

Note that for simplicity we assume linearity of the motorsensory mapping **Q**. However, the simple linearity assumption inherent in **Q** need not be inconsistent with the existence of nonlinearities between motor neuron activity and behavioral output (for example, song) and also with non-linearities between behavioral output and sensory responses. While it is the case that each of these transformations is highly non-linear, the dimensionality of motor behavior patterns realizable by muscle activity, or

**FIGURE 2 | Delayed feedback and inverse model, illustrated by vocal production in birds.** In our model of delayed sensory feedback the auditory response *a(t)* in a sensory area at time *t* depends linearly on motor activity *m(t* − *τ)* in a motor brain area at an earlier time *t* − *τ* according to *a(t)* = **Q***m(t* − *τ)*, where **Q** is the unknown motor-sensory mapping and *τ* the unknown delay of auditory feedback. An inverse **V** is a mapping from sensory neurons back onto motor neurons that inverts the action of **Q**: **V** = **Q**−1.

recorded by early sensory responses, is much smaller than the dimensionality of sensory or motor activity patterns deep within the cortex, by virtue of the fact that cortical motor and sensory neurons largely outnumber the few muscles and sensory receptors involved in the composite motor to sensory feedback loop. So for example, within the bird song system, it is thus probable that the low dimensional, composite non-linear transformation from cortical motor patterns, to muscle activity in the syrinx, to song, to cochlear response, back to cortical sensory feedback, could be well-approximated by a direct high dimensional linear map from the cortical motor area back to the cortical sensory area. This is in exact analogy to the theory of support vector regression approaches from machine learning, in which low dimensional non-linear maps can be well-approximated by high dimensional linear maps (Smola and Schölkopf, 2004). Thus, for our purposes, all we assume is that there exists at least one high dimensional linear map from cortical motor patterns to cortical sensory feedback patterns that approximates the composite feedback pathway implemented through the non-linear processes of motor generation and perception.

Now, an inverse model in this context is a mapping **V** = **Q**−1 expressed in the synaptic weights **V** from sensory onto motor neurons. Such a mapping allows sensory neurons to postdict the possible motor cause *ma* of a sensory target (vector) *a* (either driven externally, recalled from memory, or resulting from a planning strategy) according to *ma* = **V***a*. Such a postdiction ability of inverse models can be used in feedforward motor control in which the appropriate stream of motor commands *ma(t)* can be computed for a given desired sensory target sequence *a(t)* according to *ma(t)* = **V***a(t)*.

The goals of our theory are to outline a biologically plausible, local mechanism for learning of the synaptic mapping **V** and to characterize the associated emergence of mirror neurons in this process.

#### **ELIGIBILITY-WEIGHTED HEBBIAN LEARNING**

We designed a simple learning rule in which potentiation of sensory-to-motor synaptic connections **V** arises from correlated firing in pairs of sensory and motor neurons. Because sensory feedback is delayed, synapses must be able to detect correlated firing within some non-zero time window, which we achieve by introducing an eligibility trace *e(s)* that establishes a link between activity at time *t* in a motor neuron and activity in a sensory neuron at a later time *t* + *s* (see also **Figures 3A,C**). The eligibility trace modulates the change in synaptic strength associated with correlated pre- and postsynaptic firing—it is a biophysical process that resides on the postsynaptic side of **V** synapses and is triggered by activity (i.e., spikes) in the postsynaptic (motor) neuron. Intuitively, we imagine that the spiking of a motor neuron, elicited for example from an internal source of motor variation that generates exploratory motor behavior, triggers the eligibility trace that in turn makes all synapses from sensory neurons onto that motor neuron eligible for future modification. Thus, if the delayed sensory feedback arrives to the sensory area within the window of eligibility, sensory to motor synapses can potentially learn to postdict the motor cause by correlating the current sensory feedback with past motor activity that might have generated it. We further assume that the eligibility is monotonically decaying in time, implying that sensory inputs preferentially connect onto motor neurons that were recently and reliably activated rather than motor neurons that were activated a long time ago. Necessarily, the decay of the eligibility trace must be slow enough to be able to attribute significant eligibility to sensory inputs with motor-to-sensory delays *τ* , which we subsume in the condition *e(τ)* -0.

The full correlational learning rule describing changes in synaptic strength V*ij* from auditory neuron *j* onto motor neuron *i* reads:

$$\text{V}\,\text{8V}\_{i\dot{j}} = \int\_0^\infty \text{ds} \left[ e(s)m\_i(t-s)a\_{\dot{j}}(t) \right] - \hat{m}\_i(t)a\_{\dot{j}}(t),\qquad(1)$$

where *m*ˆ *i(t)* = *k* V*ikak(t)* is the (silently) postdicted motor activity, corresponding to the summed auditory input to neuron

**FIGURE 3 | Cross-correlation functions for variable and stereotyped motor codes. (A)** In a variable motor code *m(t)* (shaded area). Activity bursts *m*1 (black) and *m*2 (blue) of width *t*0 in two example motor neurons occur at diverse time lags relative to each other across renditions of the song motif. Auditory tuning in the shown sensory neuron is such that it responds *a*1 to bursts *m*1 after a time lag *τ* . Repeated co-activation *m*1 → *a*1 and non-zero eligibility *e(τ)* (red bar) at time lag *τ* leads to increased synaptic weight V11 (red arrow) and to a causal inverse. Lack of correlation between *m*2 and *a*1, as well as heterosynaptic competition, prevents V21 from similarly increasing (blue thin arrow). **(B)** The cross-correlation function *Cij(t)* for variable codes is flat except the

auto-correlation peak at zero time lag (motor activity is uncorrelated among neuron pairs). Note: based on square activity pulses in motor neurons in **(A)** the true cross-correlation shape is triangular (blue dotted line) which we approximate by a square pulse of width *t*0 10 *ms*. The auto-correlation peak height is *C*0. **(C)** In a stereotyped motor code *m(t)* (shaded area), bursts *m*1 (black) and *m*2 (blue) occur at a fixed time lag relative to each other across renditions of the song motif (traveling pulse of activity). Repeated co-activation *m*2 → *a*1 at higher eligibility (red bar) than the eligibility of *m*1 → *a*1 leads to strengthening of synapse V21 (red arrow) and to a predictive inverse. **(D)** The cross-correlation function *Cij(t)* for stereotyped codes peaks also at non-zero time lags.

*i* at time *t*. The subtractive term *m*ˆ *iaj* provides an equal time heterosynaptic depression (Lynch et al., 1977; Chistiakova and Volgushev, 2009) among all sensory afferent synapses onto a motor neuron. The strength of this depression depends on the amount of presynaptic activity but does not depend directly on postsynaptic activation. The utility of such depression is not only to stabilize activity but also to force synaptic connections towards inverse mappings as we will see. Note that we assume **V** synapses are "silently" correlating pre- and postsynaptic activity in Equation 1, i.e., **V** synapses do not contribute either to postsynaptic depolarization or to postsynaptic hyperpolarization. In other words, while the inverse model is being learned, the motor activity *mi(t)* is entirely driven by some other source than the afferent auditory input. Thus, from the perspective of extracellular physiology, it would appear that sensory feedback arriving to the motor area through the inverse model, is gated out of the motor area while that motor area is engaged in internally generated motor explorations.

In the following we examine the outcome of this learning rule in response to various forms of motor codes, with the goal of computing the synaptic weight matrix **V** at a steady-state of the learning rule, d dt**V** = 0, where denotes averaging over time (e.g., over different renditions of the song). To simplify the calculations, we assume motor codes with narrow spiketrain cross-correlation functions, i.e., the width *t*0 of spike-train cross-correlation functions is much smaller than the characteristic decay time of the eligibility trace. Although such functions have not been extensively studied due to the difficulty of simultaneously recording from several neurons during singing, narrow cross-correlation is plausible for RA and HVC neurons because pseudo simultaneous recordings can be constructed from serial recordings thanks to high firing stereotypy in these cells, yielding cross-correlation widths on the order of 10 ms (Leonardo and Fee, 2005). Note that in LMAN, because of high firing variability, similar estimation of cross-correlation width is virtually impossible.

We model motor codes with diverse inherent levels of randomness. We model stereotyped motor codes by assuming that spike-train cross correlations extend over large time lags, in agreement with a traveling pulse of activity (Hahnloser et al., 2002; Harvey et al., 2012). We model variable motor codes by assuming that cross-correlations vanish except in a peak at zero time lag (white noise assumption), **Figure 3B**.

#### **A VARIABLE NEURAL CODE YIELDS CAUSAL INVERSES**

If motor activity is uncorrelated among different neuron pairs, the resulting sensory to motor map **V** = *e(τ)t*0**Q**−1 equals the inverse of **Q** weighted by the eligibility at time lag *τ* (Equation A5, for the derivation see Appendix A3). Hence, **V** is a causal inverse that maps sensory representations onto their motor causes (in **Figure 3A**, auditory neurons map onto those motor neurons whose firing correlates most strongly with their own).

For example, during singing the motor cause *m*1*(t* − *τ)* (say a neuron that generates a 4 kHz tone) will frequently be followed by auditory response *a*1*(t)* (a 4 kHz detector neuron), leading to strengthening of synapse V11. By contrast, due to high variability of the motor code, associations between *m*2*(t* − *τ)* (say a neuron that generates a 3 kHz tone) and *a*1*(t)* are much less frequent (because the bird randomizes the production of 3 and 4 kHz tones). Hence, synapse V21 from the 4 kHz detector onto the 3 kHz generator will lose to synapse V11 due to heterosynaptic competition (**Figure 3A**).

#### **A STEREOTYPED NEURAL CODE YIELDS PREDICTIVE INVERSES**

If the motor code is stereotyped and different motor neuron pairs are correlated at even very large time lags (extending over the full range of the eligibility trace and possibly beyond), then **V** *e(*0*)t*0**H***τ***Q**−1 is approximately a concatenation of the inverse of **Q** and a shifter matrix **H***τ* that maps motor activity at one time onto motor activity at a time lag *τ* later, i.e., **V** is a predictive inverse of **Q** (Equation A10). Under a predictive inverse **V**, a sensory neuron maps onto those motor neurons that were most recently active (and reliably follow in activation other motor neurons that give rise to the sensory neuron's response).

For example, during singing, the motor cause *m*1*(t* − *τ)* of a 4 kHz tone will frequently occur before the cause *m*2*(t)* of a 3 kHz tone (because the bird produces stereotyped downsweep syllables). Hence, the 4 kHz auditory detector response *a*1*(t)* will find much higher eligibility in motor neuron 2, leading to strengthening of V21 at the expense of V11, i.e., the 4 kHz detector neuron connects onto the 3 kHz generator neuron (**Figure 3C**).

#### **LACK OF RESPONSE TO PERTURBED AUDITORY FEEDBACK AND SELECTIVITY FOR THE BOS**

During Hebbian learning of **V** in Equation 1 we required that synapses **V** are not able to drive spike responses in motor neurons during singing (**V** synapses learn silently). The main intuitive reason for the necessity of silent learning is that the learning goal of the inverse model synapses are to silently correlate the motor and sensory streams, without perturbing the motor stream that would result if sensory feedback were to pass through and drive spikes in the motor area. If the inverse model synapses allowed sensory feedback to significantly drive motor spikes, then the incoming sensory signals would serve to drive motor activity resulting in cyclic motor output with cycle time approximately equal to *τ* , i.e., birds would unavoidably produce repetitive motor output (stuttering).

Interestingly, there is much evidence for the gating out of sensory information in song motor nuclei. Principal motor neurons in LMAN and HVC do not respond to playback of white noise stimuli during singing (Leonardo, 2004; Kozhevnikov and Fee, 2007) and during states of high arousal (Cardin and Schmidt, 2003), though there are reports of distorted feedback responses in HVC interneurons in Bengalese finches (Sakata and Brainard, 2008). Lack of feedback sensitivity in principal motor neurons is usually ascribed to a form of gating caused by specific thalamic or neuromodulatory mechanisms (Dave et al., 1998; Schmidt and Konishi, 1998; Shea and Margoliash, 2003; Cardin and Schmidt, 2004; Coleman et al., 2007; Hahnloser et al., 2008), see also the Discussion.

By contrast, LMAN (Doupe and Konishi, 1991; Doupe, 1997; Solis and Doupe, 1999; Roy and Mooney, 2007) and HVC neurons (Katz and Gurney, 1981; Margoliash, 1983, 1986; Williams and Nottebohm, 1985) respond to auditory stimulation while birds are anesthetized or asleep, which we model as gating on of **V** synapses, i.e., we assume that auditory responses in motor neurons are driven via the learned inverse models.

The puzzling observation of the gating out of sensory inputs to motor areas during motor exploration, is naturally accounted for in our theory by necessity of correlating the current presynaptic sensory stream with past postsynaptic motor streams, to learn an unbiased inverse model (of unperturbed motor stream).

Also, interestingly, in both HVC and LMAN sensory responses are strongest for bird's own song (BOS) stimuli compared to other stimuli including the tutor song or the BOS played back in reverse time (McCasland and Konishi, 1981; Margoliash, 1986; Lewicki, 1996; Solis and Doupe, 1999). Such selectivity follows naturally from our model assumptions: for both stereotyped and variable motor codes, the mappings, whether causal or predictive, can only invert sensory responses that lie in the image of **Q** and cannot invert the full space of responses orthogonal to the image of **Q**. Such a restriction arises because only sensations that could arise through combinations of previously experienced sensory feedback during singing can actually be inverted into appropriate motor commands. In other words, the inverse model synapses map prior sensory feedback generated by the bird's own previous song into appropriate motor commands, but necessarily fails to map sensory activity patterns that are very different from the BOS into coherent motor patterns. Thus, assuming HVC and LMAN can be thought of as downstream of the output of an inverse model, our Hebbian learning rule generating inverse models can naturally account for the preference of sensory responses in HVC and LMAN for BOS; sounds very different from BOS are not appropriately inverted, and therefore presumably do not lead to coherent activation of motor patterns via sensory inputs propagating through the inverse model synapses.

#### **INVERSE MODELS AND SENSORIMOTOR MIRRORING**

The Hebbian learning rule in Section Eligibility-Weighted Hebbian Learning determines the wiring of sensory afferents into motor areas based on sensorimotor experience. How could one experimentally test for the existence of such wiring without painstaking, detailed inspection of anatomical connections and characterization of the sensorimotor mapping **Q**? Here we outline the design of experiments to probe for the existence of either causal or predictive inverse models. We propose to record from single neurons both in sensory and motor states and to compare motor activity and sensory-evoked responses using cross-correlation functions: as we will show, the time lag of peak cross correlation provides evidence for either predictive or causal inverses.

In such mirroring experiments that we propose, a single neuron is first recorded during singing and then during playback of the just recorded songs while the bird is asleep in the dark (during which the auditory gate is open and motor neurons become responsive to auditory stimuli, presumably through an inverse model from an upstream sensory area). In our model, sensory responses *ma i (t)* = ˆ*mi(t)* = *k* V*ij aj(t)* during playback are driven via synaptic weights **V** (assumed to be at a steadystate of Equation 1, d dt**V** = 0). Computing the cross-correlation

functions Corr*(s)* of the sensory response *ma i (t)* with motor activity *mi(t)* (as a function of time lag *s)* yields that (see **Figure 4**):

For derivations and model assumptions see Appendices A2–A4. In particular, here we assumed no synaptic delay between auditory and motor neuron, though this assumption can be relaxed. In summary, for both stereotyped and variable motor codes, sensory responses mirror motor activity. The amount of randomness in the motor code dictates the time lag of peak cross-correlation between motor activity and sensory-evoked responses, which we refer to as the *mirroring offset*. The mirroring offset thus serves as an important experimental observable that provides a window into fundamental differences in the types of inverse models that are computed by Hebbian learning, **Figure 4**.

Note that variable motor codes are associated with weaker mirroring than stereotyped codes, i.e., the cross-correlation functions for variable codes exhibit lower peak amplitudes than crosscorrelation functions associated with stereotyped codes: In our model, the ratio of peak cross correlation is given by the eligibility at time lag *τ* divided by the eligibility at time lag zero (Equation A12 derived in the Appendices A3, A4). Thus, the steeper the

*ma*

eligibility trace, the weaker the mirrored response in case of variable motor codes. By contrast, the shape of the eligibility trace is expected to have almost no influence in case of stereotyped codes.

1 in case of a causal inverse (middle column) and to motor neuron response

2 in case of a predictive inverse (right column) after an additional time lag *τs* (spike propagation time from auditory to motor area) that is assumed to be 0

Note that the auditory response *ma i (t)* = *k* V*ik ak(t)* in a motor neuron to song playback is mathematically identical to the (silently) postdicted motor activity *m*ˆ *i(t)* = *k* V*ik ak(t)* defined after Equation 1, and used in learning the inverse model. Nevertheless, we use different symbols for these quantities to disentangle their meaning, i.e., the former being a superthreshold sensory response elicited in a quiet non-singing state of the bird, the latter being a subthreshold subtractive term that stabilizes synaptic learning during singing. The biophysical underpinnings of these two terms might largely be identical, with the silent nature of the posticted activity arising from some form of response gating (see also the Discussion).

#### **GRADIENT DESCENT**

*ma*

*ma*

We note that the learning rule in Equation 1 corresponds to gradient descent on the following error function:

$$\mathcal{E}(t) = \frac{1}{2} \sum\_{i} \int\_{0}^{\infty} \left[ m\_{i}(t - s) - \sum\_{k} \mathcal{V}\_{ik} a\_{k}(t) \right]^{2} \mathcal{e}(s) \, \text{ds} \tag{2}$$

For a derivation, see Appendix A1. Thus, synaptic weights **V** converge such as to yield optimal postdiction *m*ˆ *i(t)* = *k* V*ik ak(t)* of motor activity from sensory feedback. The origin of our eligibility-weighted Hebbian learning rule with heterosynaptic competition, from gradient descent of an energy function, confers a degree of robustness to the learning, as well as suggests generalizations to situations in which the synaptic transformation from sensory to motor areas is non-linear.

2 is selective to the sound feature that during singing was generated by the much earlier burst *m*1 in a different neuron (black burst in **(A)**, right panel), but

not the feature generated by *m*2 (blue burst in **(A)**, right panel).

#### **PROBABILISTIC MODELS**

More realistic neuron models are non-linear and contain spikes that are potentially probabilistic and certainly binary events. Also, more realistically, we may want to explicitly model intrinsic noise in motor and sensory-related responses rather than deal with motor variability only through their effects on cross correlations. As a first step to dealing with such realism, we have derived two probabilistic neuron models in which inverse models and mirroring can be studied in similar manners as in the linear model, outlined in the following.

In one of these models we calculate the influence of probabilistic (binary) responses on the strength of mirroring. We consider a random motor area that at any time can only be in one of two possible states *M* = 1 and *M* = 0 with prior probability *p(M* = 1*)* = 1 2 . Assume analogously that the sensory area is such that a particular sensory feature is either detected (*S* = 1) or not detected (*S* = 0). We then model the relationship between motor activity and sensory consequence in terms of conditional dependencies between these two random variables. We assess the strength of mirroring in this model in terms of the cross-correlation coefficient between the two random variables (as derived in Appendix A5) and find the following result:

Thus, the simple probabilistic model shows that the strength of mirroring may also be strongly reduced by the amount of intrinsic noise present in sensory and motor systems.

## **DISCUSSION**

We have presented a simple model for the development of mirror neuron systems that is mathematically tractable, allowing us to relate mirror neuron properties such as the correlative strengths and the time lag of peak mirrored responses to the stereotypy (the correlation structure) of motor-related firing. Mirroring properties depend on the variability of the neural motor code which may be dissociated from apparent variability of the motor behavior as is the case in LMAN neurons that fire highly variable spike patterns despite high song stereotypy in adults. Our conclusions are valid for arbitrary sensory systems, provided they are able to signal sensory feedback from motor actions with sufficient sensitivity matched to the behavioral richness generated by the motor system (and of course provided that sensory afferents are subject to correlative Hebbian learning). In our derivation we have assumed that cross-correlation functions among motor neuron pairs are narrow, which was a simplifying assumption that allowed us to derive simple analytical forms of the sensory-to-motor mapping **V** and of mirroring properties. Approximate inverses should also result for motor codes with more complex time dependence, because by construction, the learning rule we considered corresponds to a gradient-descent rule that achieves minimal inversion error.

Although inverse models are attractive as models for vocal learning (Guenther et al., 2006; Hahnloser and Ganguli, 2013), they have previously been judged to be inappropriate for vocal learning in songbirds because of mainly two reasons: (1) young birds require many song repetitions with auditory feedback (Doya and Sejnowski, 2000), and (2) the learning schemes proposed either used a biologically implausible algorithm (Jordan and Rumelhart, 1992) or assumed the preexistence of an approximate inverse model (Kawato, 1990). Here we suggested a resolution to both of these issues and shown that in contrary to

previous beliefs, inverse models constitute a potentially plausible framework for vocal learning in birds, too: the many song explorations used by young birds could be required to actually learn the high dimensional inverse model; and, the correlational learning we proposed is quite plausible and simple (but non-trivial nevertheless). This suggests potentially opening up the hypothesis space for learning rules operating within cortico-basal ganglia circuits, in both mammalian and bird song systems, to include models spanning the range from pure reinforcement learning (RL) to pure inverse model learning. Of particular interest would be intermediate learning rules that synergistically incorporate both dopamine-dependent plasticity thought to underlie RL as well, as Hebbian based plasticity shown here to mediate inverse model learning, in order to implement sophisticated model-based RL strategies. For example, a simple proposal would be that dopamine delivered to striatal synapses from the ventral tegmental area (VTA) might not be released purely nonspecifically, but instead might be delivered by an inverse model that can partially map errors in sensory coordinates to errors in motor coordinates, thereby guiding learning in ways more sophisticated than pure RL (O'Reilly and Frank, 2006).

The key to learning causal inverse models is motor variability. In motor areas such as HVC that fire stereotyped patterns, auditory afferents cannot disentangle cause-and-effect, leading to preferential formation of predictive inverses rather than causal ones. Predictive inverses have limited usefulness for action imitation from action observation, because under a predictive inverse, observation of a particular motor gesture will lead to imitation of the subsequent gesture in the imitator's motor repertoire, which may not be part of the actions to be imitated. For example, if a bird repeatedly sings *ABCD* during formation of the inverse and wants to later imitate repetitions of *ABDB*, then its predictive inverse will constrain it to produce repetitions of *BCDA* because perception of *A* maps to production of *B*, perception of *B* maps to production of *C*, etc.

Small temporal delays between motor activity and activity evoked by playback of BOS or BOS-resembling sounds have been reported previously. Prather et al. (2008) showed there is a small mirroring offsets of just a few milliseconds in HVCX neurons of awake swamp sparrows and report similar (not quantified) results in Bengalese finches. Furthermore, Dave and Margoliash (2000) observed a small time lag of auditory-evoked activity also in RA neurons of sleeping zebra finches. Both these experimental findings reflect a predictive inverse. While predictive inverses have limited usefulness for action imitation they might provide stability in sequential vocalization. Indeed Sakata and Brainard report that perturbation of auditory feedback can change song syntax in Bengalese finches (Sakata and Brainard, 2006, 2008; Hanuschkin et al., 2011).  Indeed, preliminary results indicated a large non-zero mirroring offsets in LMAN (Giret et al., 2012).

An important element of our theory is the eligibility trace. To endow Hebbian learning with such a trace is necessary in realistic situations in which effects (sensory feedback) follow their cause (motor command) with some non-zero time lag arising from signal propagation delays, from twitch times of muscles, and from sensory and synaptic receptor latencies. In humans such a lag could span up to several hundreds of milliseconds, whereas in birds it may be as short as several tens of milliseconds. Eligibility traces also appear in RL theories (Seung, 2003; Fiete et al., 2007) and seem to be a general prerequisite for learning in the context of delayed feedback or delayed reward. We can imagine that neurons and synapses may hold decaying eligibility traces in terms of dedicated molecules such as calcium. Action potential generation is associated with rapid calcium entry that decays over the time course from several hundreds of milliseconds to a few seconds (McGeown et al., 1996; Wallace et al., 2008). The monotonic decay of intracellular calcium is well-suited to modeling a monotonically decaying eligibility trace. However, a monotonic decay of eligibility harbors both advantages and disadvantages. The disadvantage, as discussed, is the problem associated with stereotyped motor generators that can only hold predictive inverse models; to make inverse models causal, motor variability is required. Another way to guarantee causal inverse models—even under stereotyped motor explorations—would be to consider eligibility traces that do not monotonically decay but that peak at precisely the time delay inherent in closed sensorimotor feedback loops. The main caveat of such eligibility traces is that it may be questionable whether different muscles recruited for the same behavior must necessarily be associated with the same sensorimotor delay—and it is presently unclear how such variable delays could be matched to variable eligibility traces across synapses in a way that would ensure the learning of a causal inverse model. Moreover, phenomena such as speech co-articulation make it unlikely that there exists a constant sensorimotor delay across a large range of premotor neurons. The advantage, on the other hand, of a decaying eligibility trace is that sensorimotor contingencies and inverses can be learned regardless of sensorimotor latencies, providing robustness of sensorimotor learning.

Convergence of the sensory to motor synaptic weights toward inverses depends on details of the heterosynaptic competitive term. Heterosynaptic competitive terms have a certain appeal because of the useful normalization they provide (Fiete et al., 2010). In the context of this work, such terms imply locally available information at a single synapse about sensory inputs to other synapses. Though this information need not be provided instantaneously, we can only speculate about possible mechanisms for sharing such information among different synapses onto the same postsynaptic neuron. One possibility is that some form of intracellular signaling conveys this information. Another possibility to be explored is whether there exists an entire class of such competitive terms with a similar effect. For example, provided that motor and sensory codes are sufficiently sparse, it is conceivable that very simple subtractive terms might suffice for inverse formation. Whether other (even simpler) competitive terms result in approximate inverses needs to be further explored. We would like to point out preliminary evidence that inverses can be learned with Hebbian rules that include no heterosynaptic competitive terms (Senn and Pawelzik, pers. communication).

Our Hebbian learning theory has been analyzed so far in linear circuits, but we have indicated ways to overcome linearity by pinpointing extensions of our work to include nonlinear mappings and probabilistic neuron models. Further work will be required to test whether our correlative learning approach is suitable also for inverse model learning employing detailed biophysical models of the avian syrinx.

## **ACKNOWLEDGMENTS**

We acknowledge support by the European Research Council (ERC-Advanced Grant 268911) and the Swiss National Science Foundation (Grant 31003A\_127024), and support from the Swartz, Sloan, and Burroughs-Wellcome foundations, and Defense Advanced Research Projects Agency (DARPA). R. H. R. Hahnloser thanks Walter Senn for helpful discussions on inverse models and S. Ganguli thanks Michael Brainard and Kris Bouchard for useful discussions on birdsong learning, and the transition from causal to predictive inverse models.

#### **REFERENCES**

*Science* 290, 812–816. doi: 10.1126/science.290.5492.812

S. Gazzaniga (Cambridge, London: MIT Press), 469–484.

syntax generation in the Bengalese finch. *J. Comput. Neurosci.* 31, 509–532. doi: 10.1007/s10827-011- 0318-z

white-crowned sparrow. *J. Neurosci.* 3, 1039–1057.

Mooney, R. (2012). Motor circuits are required to encode a sensory model for imitative learning. *Nat. Neurosci.* 15, 1454–1459. doi: 10.1038/nn.3206

*J. Neurophysiol.* 104, 2474–2486. doi: 10.1152/jn.00977.2009

Wallace, D. J., Meyer, S., Astori, S., Yang, Y., Bausen, M., Palmer, A. E., et al. (2008). Single-spike detection *in vitro* and *in vivo* with a genetic Ca2+ sensor. *Nat. Methods* 5, 797–804. doi: 10.1038/nmeth.1242

Williams, H., and Nottebohm, F. (1985). Auditory responses in avian vocal motor neurons: a motor theory for song perception in birds. *Science* 229, 279–282.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 01 November 2012; accepted: 15 May 2013; published online: 19 June 2013.*

*Citation: Hanuschkin A, Ganguli S and Hahnloser RHR (2013) A Hebbian learning rule gives rise to mirror neurons and links them to control theoretic inverse models. Front. Neural Circuits 7:106. doi: 10.3389/fncir. 2013.00106*

*Copyright © 2013 Hanuschkin, Ganguli and Hahnloser. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in other forums, provided the original authors and source are credited and subject to any copyright notices concerning any third-party graphics etc.*

## **APPENDICES**

## **A1. GRADIENT DESCENT DERIVATION OF ELIGIBILITY-WEIGHTED HEBBIAN LEARNING**

We can derive the learning rule in Equation 1 by gradient (steepest) descent on an error function E. The differential change δ**V** in synaptic weight is proportional to the gradient and we can write:

$$
\delta \mathbf{V} = -\frac{\mathrm{dE}}{\mathrm{d} \mathbf{V}}.\tag{A1}
$$

The error function E*i(t)* for neuron *i* we define as the square difference between motor activity *mi(t* − *s)* and postdicted motor activity *m*ˆ *i(t)* = *k* V*ik ak(t)*, weighted by the eligibility associated with the time lag *s*.

$$\mathcal{E}\_i(t) = \frac{1}{2} \int\_0^\infty \mathrm{d}s \left[ m\_i(t-s) - \sum\_k \mathcal{V}\_{ik} \, a\_k(t) \right]^2 \, \mathcal{e}(s)$$

The total error is simply the sum of errors over all neurons E = *i Ei*.

By taking the gradient with respect to the *i,j* th weight only we find

$$\frac{d\mathcal{E}\_i}{d\mathcal{V}\_{\vec{\eta}}} = -\int\_0^\infty \text{ds} \left[ m\_i(t-s) - \sum\_k \mathcal{V}\_{ik} a\_k(t) \right] e(s) a\_j(t)$$

$$= -\int\_0^\infty \text{ds} \left[ e(s) m\_i(t-s) \left. a\_j(t) + \sum\_k \mathcal{V}\_{ik} a\_k(t) \right| a\_j(t) e(s) \right].$$

Assume a normalized eligibility trace ∞ 0 *e(s)*ds = 1 :

$$\frac{\mathrm{dE}\_i}{\mathrm{dV}\_{ij}} = -\int\_0^\infty \mathrm{ds}\,\boldsymbol{e}(\boldsymbol{s}) m\_i(t-\boldsymbol{s}) \,\boldsymbol{a}\_j(t) + \sum\_k \mathrm{V}\_{ik} \,\boldsymbol{a}\_k(t) \,\boldsymbol{a}\_j(t).$$

$$\Rightarrow 8\mathrm{V}\_{ij} = \int\_0^\infty \mathrm{ds}\,\boldsymbol{e}(\boldsymbol{s}) m\_i(t-\boldsymbol{s}) \,\boldsymbol{a}\_j(t) - \sum\_k \mathrm{V}\_{ik} \,\boldsymbol{a}\_k(t) \,\boldsymbol{a}\_j(t) \,\text{ (A2)}$$

## *A1.1. Extension to non-linear network*

Note that our linear approach can be extended by introducing a nonlinear function *f* in the auditory to motor mapping in Equation 2:

$$\mathcal{E} = \frac{1}{2} \sum\_{i} \int\_{0}^{\infty} \mathrm{d}s \left[ m\_{i}(t - s) - f \left( \sum\_{k} \mathcal{V}\_{ik} \, a\_{k}(t) \right) \right]^{2} e(s)$$

$$\Rightarrow \delta \mathcal{V}\_{ij} = \int\_{0}^{\infty} \mathrm{d}s \Big[ m\_{i}(t - s) - f\_{i} \Big] f\_{i}^{'} \, a\_{k}(t) e(s),$$

$$= \int\_{0}^{\infty} \mathrm{d}s \, m\_{i}(t - s) f\_{i}^{'}(t) \, a\_{k}(t) e(s) - f\_{i}(t) f\_{i}^{'}(t) \, a\_{k}(t)$$

Where *fi* = *f k* V*ik ak(t)* and *f i* = d*f k* V*ik ak(t) /*dV*ij*.

## *A1.2. Probabilistic derivation of Hebbian learning rule*

We derive a version of the Hebbian learning rule in Equation 1 that is based on the following probabilistic Boltzmann neuron model. For simplicity, we do not include the time dependence in the derivation (*τ* = 0). The auditory feedback response *a* given a motor activation *m* is given by the conditional probability

$$\mathbb{P}\_{\mathbf{Q}}(a|m) = \frac{\ell^{a^T \mathbf{Q}m}}{Z\_{\mathbf{Q}}(m)},$$

parameterized by the matrix **Q**, which is the motor-sensory mapping (as before) and where

$$Z\_{\mathbf{Q}}(m) = \sum\_{a} e^{a^T \mathbf{Q}^m}$$

is the partition function. The posterior probability of *m* is given by

$$\mathbb{P}\_{\mathbf{Q}}(m|a) = \frac{\mathbb{P}\_{\mathbf{Q}}(a|m)\mathbb{P}(m)}{\mathbb{P}(a)}.$$

In a sensory state, auditory responses in motor neurons are driven via synapses **V** according to the probabilistic model:

$$\mathbf{P}\_{\mathbf{V}}(m|a) = \frac{e^{m^T \mathbf{V}\_a}}{Z\mathbf{v}(a)}$$

with partition function

$$Z\mathbf{v}(a) = \sum\_{m} e^{m^T \mathbf{V} a}.$$

The error function in Equation 2 is replaced by the Kullbach-Leibler (KL) divergence between P**Q***(m*|*a)* and P**V***(m*|*a)*:

$$\begin{aligned} \mathrm{D\_{KL}}\left(\mathrm{P\_{Q}}(m|a), \mathrm{P\_{V}}(m|a)\right) &= \sum\_{m} \mathrm{P\_{Q}}(m|a) \ln \left(\frac{\mathrm{P\_{Q}}(m|a)}{\mathrm{P\_{V}}(m|a)}\right) \\ &= \sum\_{m} \mathrm{P\_{Q}}(m|a) \left[\ln \left(\mathrm{P\_{Q}}(m|a)\right)\right] \\ &\quad - \ln \left(\mathrm{P\_{V}}(m|a)\right) \left[\\ &= \sum\_{m} \mathrm{P\_{Q}}(m|a) \left[\ln \left(\mathrm{P\_{Q}}(m|a)\right)\right] \\ &\quad + \ln \left(Z\mathrm{y}\_{l}(a)\right) - m^{T}\mathrm{V}a \right] \quad \text{(A)} \end{aligned}$$

Before taking the derivative of DKL we compute the derivative of the partition function:

$$\frac{\partial}{\partial \mathbf{V}\_{\vec{\boldsymbol{w}}}} Z\_{\mathbf{V}}(a) = \sum\_{m} e^{m^{T} \mathbf{V} a} \frac{\partial}{\partial \mathbf{V}\_{\vec{\boldsymbol{w}}}} m^{T} \mathbf{V} a = \sum\_{m} e^{m^{T} \mathbf{V} a} m\_{i} a\_{j}$$

$$= Z\_{\mathbf{V}}(a) a\_{j} \sum\_{m} \mathbf{P}\_{\mathbf{V}}(m|a) m\_{i}$$

$$= Z\_{\mathbf{V}}(a) a\_{j} \langle m\_{i} | a \rangle\_{m},$$

based on which it follows that

$$\frac{\partial}{\partial \mathcal{V}\_{\vec{\eta}}} \ln(Z\mathsf{v}(a)) = \frac{1}{Z\mathsf{v}(a)} \frac{\partial}{\partial \mathcal{V}\_{\vec{\eta}}} Z\mathsf{v}(a) = a\_{\vec{\eta}} \langle m\_{i}|a\rangle\_{m} \dots$$

Using this relationship and Equation A3 we can calculate the derivative of the KL-divergence with respect to V*ij*:

$$\begin{split} &\frac{\partial}{\partial \mathbf{V}\_{ij}} \mathrm{D\_{KL}}(\mathbf{P}\_{\mathbf{Q}}(m|a), \mathbf{P}\_{\mathbf{V}}(m|a)) \\ &= \frac{\partial}{\partial \mathbf{V}\_{ij}} \sum\_{m} \mathrm{P}\_{\mathbf{Q}}(m|a) \left[ \ln \left( \mathrm{P}\_{\mathbf{Q}}(m|a) \right) + \ln \left( Z\_{\mathbf{V}}(a) \right) - m^{T} \mathbf{V} a \right] \\ &= \sum\_{m} \mathrm{P}\_{\mathbf{Q}}(m|a) \left( \frac{\partial}{\partial \mathbf{V}\_{ij}} \left( \ln \left( Z\_{\mathbf{V}}(a) \right) \right) - m\_{i} a\_{j} \right) \\ &= \sum\_{m} \frac{\mathrm{P}\_{\mathbf{Q}}(a|m) \mathrm{P}(m)}{\mathrm{P}(a)} \left( \frac{\partial}{\partial \mathbf{V}\_{ij}} \left( \ln \left( Z\_{\mathbf{V}}(a) \right) \right) - m\_{i} a\_{j} \right) \\ &\Rightarrow \left\langle -\frac{\partial}{\partial \mathbf{V}\_{ij}} \mathrm{D}\_{\mathrm{KL}} \left( \mathrm{P}\_{\mathbf{Q}}(m|a), \mathrm{P}\_{\mathbf{V}}(m|a) \right) \right\rangle\_{a} \\ &= \sum\_{a,m} \mathrm{P}\_{\mathbf{Q}}(a|m) \mathrm{P}(m) \left( m\_{i} \ a\_{j} - \frac{\partial}{\partial \mathbf{V}\_{ij}} \left( \ln \left( Z\_{\mathbf{V}}(a) \right) \right) \right) \end{split}$$

Thus, the gradient decent leads to

$$0 \Rightarrow \delta \mathbf{V}\_{i\bar{j}} = m\_i a\_{\bar{j}} - \frac{\partial}{\partial \mathbf{V}\_{i\bar{j}}} (\ln \left( Z\_{\mathbf{V}}(a) \right)) = \left[ m\_i - \langle m\_i | a \rangle\_m \right] a\_{\bar{j}}$$

This is the probabilistic analog of Equation A2, in which the silent postdictive motor activity *m*ˆ *i* is replaced by the conditional expectation *mi*|*a*m of activity in motor neuron *i* given the sensory response *a*.

#### **A2. CORRELATION OF MOTOR ACTIVITY DETERMINES AVERAGE SYNAPTIC CHANGE**

The average synaptic change under learning rule Equation A2 satisfies

$$
\left< \mathbf{V}\_{ij} \right> = \int\_0^\infty \left< m\_i(t-s)a\_j(t) \right> \mathbf{e}(s) \mathbf{d}s - \left< \sum\_k \mathbf{V}\_{ik} a\_k(t) a\_j(t) \right>,
$$

$$
= \int\_0^\infty \left< \sum\_k m\_i(t') \, \mathbf{Q}\_{jk} \, m\_k(t'+s-\tau) \mathbf{e}(s) \right> \mathbf{d}s
$$

$$
$$

where we have substituted *t* = *t* − *s*. We can write this equation as

$$
\langle \delta \mathbf{V} \rangle = \left[ \int\_0^\infty e(s) \mathbf{C}(s - \tau) \mathrm{d}s - \mathbf{VQC}(0) \right] \mathbf{Q}^T,\qquad \text{(A4)}
$$

where C*ij(s)* = *mi (t) mj(t* + *s)* is the cross-correlation matrix of motor activity at time lag *s*. In the following we assume without loss of generality that the delay *τs* of synaptic transmission between auditory and motor neurons is negligibly small.

## **A3. VARIABLE MOTOR CODE**

#### *A3.1. V is a causal inverse*

We assume a motor code with a narrow correlation function that is non-zero only for small *t*0.

$$\mathbf{C}(\mathbf{s}) = \begin{cases} 1C\_0 \text{ for } |\mathbf{s}| < t\_0/2 \\ 0 \text{ otherwise} \end{cases}$$

Where **1** is the unity matrix and *C*0 is a positive constant. The steady state solution δ**V** = 0 of Equation A4 leads to

$$\int\_0^\infty e(s)\mathbf{C}(s-\tau)\,\mathrm{ds}-\mathbf{VQC}(0) = 0.$$

Assuming that the eligibility trace is constant over short time intervals of duration *t*0 (over which the correlation function is non-zero) yields

$$\begin{aligned} \Rightarrow e\left(\tau\right)t\_0 \mathbf{1}C\_0 - \mathbf{V} \mathbf{Q} \mathbf{1}C\_0 &= 0 \\\\ \iff \mathbf{V} \mathbf{Q} &= e(\tau)t\_0 \mathbf{1} \end{aligned}$$

This implies that the auditory to motor mapping **V** is proportional to the inverse of **Q** weighted by the eligibility at time lag *τ* :

$$\mathbf{V} = e(\boldsymbol{\tau}) \boldsymbol{t}\_0 \mathbf{Q}^{-1} \tag{A5}$$

Thus, **V** is a causal inverse, at least when restricted to the image of **Q**.

## *A3.2. Variable motor codes are associated with large mirroring offsets*

We simulate a mirroring experiment in which we cross correlate in a given neuron the motor activity *mi(t)* and the activity *ma i (t)* that results from observation of the motor act (achieved in birds by song playback though a loudspeaker).

The auditory response in motor neuron *i* is given by

$$\begin{aligned} m\_i^a(t) &= \sum\_j \mathbf{V}\_{ij} a\_j(t) \\ &= \sum\_j (\mathbf{V}\mathbf{Q})\_{ij} m\_j(t-\tau) \\ &= e(\tau) t\_0 \sum\_j \delta\_{i,j} m\_j(t-\tau) .\end{aligned}$$

$$\begin{aligned} &= e(\tau) t\_0 \, m\_i(t-\tau), \end{aligned}$$

Where *δi,j* is the Kronecker-Delta, (*δi,j* = 1 for *i* = *j* and *δi,j* = 0 otherwise). Note that relative to song (either produced by the bird or played through the loudspeaker) the playback-evoked activity *ma i (t)* is shifted with respect to the motor activity *mi(t* − *τ)* by a time shift *τ* (as illustrated in **Figure 4**). The cross correlation Corr*(s)* between sensory-evoked and motor generated activity is defined as,

$$\text{Corr}(\mathbf{s}) = \frac{1}{T'} \int\_0^{T'} m\_i(t) m\_i^a(t+s) \, \text{d}t = \left< m\_i(t) m\_i^a(t+s) \right> \text{(A6)}$$

Where *T* is the duration of the motor behavior (e.g., the song motif or song). Inserting the expression for the motor activity evoked by the auditory response into Equation A6 yields,

$$\begin{aligned} \text{Corr}(\mathbf{s}) &= \langle m\_i(t)\mathbf{e}(\mathbf{r})t\_0 \, m\_i(t-\tau+s) \rangle \\ &= \begin{cases} e(\mathbf{r})t\_0 \mathbf{C}\_0 & \text{for } |\mathbf{s}-\tau| < t\_0/2 \\ 0 & \text{otherwise} \end{cases} \end{aligned} \tag{A7}$$

Thus, the cross correlation is nonzero in a small time centered around *τ* . The peak cross-correlation value is given by the eligibility trace at time lag *τ* .

$$\text{CorrPeak} = t\_0 C\_0 e(\mathbf{r})$$

Note that our calculations are valid in principle for meansubtracted *mi(t)*. In case of non-mean subtracted *mi* we should replace the cross correlation in Equation A6 by the cross covariance to obtain the same findings. However, in practice, mean subtraction is not necessary because the peak location (the mirroring offset) is independent of the mean.

#### **A4. STEREOTYPED MOTOR CODES**

#### *A4.1. V is a predictive inverse*

We describe the motor activity by a traveling pulse *mi(t)* = η *i* ω − *t* with speed ω, where

$$\eta(t) = \begin{cases} 1 & \text{for } |t| < t\_0/2 \\ 0 & \text{otherwise} \end{cases}$$

and *t*0 = 1*/*ω. The cross-correlation matrix for such a traveling pulse is a triangular pulse of height *t*0*/T* and width 2*t*0 which we approximate by a square pulse of width *t*0*,*C*ij(s) t*0 *T* η *i* − *j* ω + *s* in the following (illustrated in **Figure 3D**). To facilitate comparison of inverses associated with stereotyped and variable motor codes, we assume their peak correlations are identical, i.e.,

$$C\_0 = \frac{t\_0}{T}.$$

At a steady state δ**V** = 0, Equation A4 yields

$$\int\_0^\infty e\left(s\mathbf{C}(s-\tau)\mathbf{ds} - \mathbf{VQC}(0) = 0\right)$$

Assuming again that the eligibility trace is constant over short time intervals of width *t*0*,* we find

$$\Rightarrow e\left(\mathbf{r} - \frac{i - j}{\alpha}\right)t\_0 = \mathbf{VQ}.$$

Given sensory input *aj(t)* and the inverse map **V**, the auditoryevoked activity *ma j (t)* is proportional to the eligibility trace:

$$m\_i^a(t) = \sum\_j \mathbf{V}\_{ij} a\_j(t)$$

$$= \sum\_j (\mathbf{V}\mathbf{Q})\_{ij} m\_j(t - \tau)$$

$$= e\left(t - \frac{i}{\alpha}\right) t\_0,\tag{A8}$$

defined for *t* ≥ *i* ω .

By approximating the eligibility trace only by its maximum value *e t* − *i* ω = *e(*0*)*η *t* − *i* ω we have that approximately

$$m\_i^d(t) \cong e(0)t\_0 \, m\_i(t),\tag{A9}$$

and so the playback-evoked activity *ma i (t)* is roughly identical to motor activity *mi(t)* (as illustrated in **Figure 4**).

To compute the matrix **V** we use the same approximation for the eligibility trace to obtain

$$\frac{1}{t\_0}(\mathbf{V}\mathbf{Q})\_{\vec{\eta}} \simeq e(0)\eta(i-j-\tau) = e(0)\mathbf{H}^t\_{\vec{\eta}},$$

where we have set ω = 1, and where **H** is a shifter matrix (also called cyclic permutation or circulant matrix), e.g., for *n* = 4:

$$\mathbf{H} = \begin{pmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 1 & 0 & 0 & 0 \end{pmatrix}.$$

With this approximation we find that the synaptic mapping

$$\mathbf{V} \simeq e(0)t\_0 \mathbf{H}^\mathbf{f} \mathbf{Q}^{-1} \tag{A10}$$

is the inverse of the motor map shifted in time by *τ* , i.e., **V** maps sensory activity evoked by motor activity at time *t* onto motor activity at time *t* + *τ* . In other words, the sensory lag is compensated and sensory-evoked motor activity at time *t* predicts motor activity at time *t*. Hence, **V** is a predictive inverse.

#### *A4.2. Stereotyped motor codes are associated with small mirroring offsets*

Based on the sensory-evoked activity *ma i (t)* derived in Equation A8 we find for the cross-correlation function Corr*(s)* between motor activity *mi(t)* and sensory-evoked activity *ma i (t)*:

$$\text{Corr}\left(s\right) = \frac{1}{T} \int\_{0}^{T} m\_{i}(t) m\_{i}^{a}(t+s)dt$$

$$= \frac{t\_{0}}{T} \int\_{0}^{T} \eta\left(\frac{i}{\omega} - t\right) \varepsilon\left(t - \frac{i}{\omega} + s\right) dt$$

$$= \frac{t\_{0}^{2}}{T} \varepsilon(s) = t\_{0}C\_{0}\varepsilon(s). \tag{A11}$$

Thus, the cross-correlation function is proportional to the eligibility trace. If the eligibility trace is monotonically decaying we find that the peak of Corr*(s)* occurs at *s* = 0 and is given by

$$\text{CorrPeak} = t\_0 C\_0 e(s).$$

In other words, stereotyped neural codes are associated with zero mirroring offsets.

The ratio of peak cross correlation for variable (A6) and stereotyped (A10) motor codes is given by

$$r = \frac{e(\mathbf{r})}{e(0)},\tag{A12}$$

implying that the more stereotyped a neural code, the stronger is the observed mirroring effect.

#### **A5. MIRRORED RESPONSE STRENGTH IN A PROBABILISTIC MODEL**

In the following we define a probabilistic model of a motor neuron that allows us to compute the mirroring strength, i.e., the correlation between motor activity and sensory-evoked activity. We assume a minimal model in which a neuron has only two states *R* = 1 (active) and *R* = 0 (inactive). In addition, we assume two behavioral states *B* = 1 (behavioral feature present), and *B* = 0 (behavioral feature absent). During motor production, the degeneracy of the motor code quantified by the conditional probability of neural activity given that the feature of interest is present during the behavior (e.g., the finger is extended or the song pitch is high) is

$$P\_M(R=1|B=1) = p\_1$$

and the probability that the neuron is active while the behavioral feature is absent (intrinsic noise) is

$$p\_M(\mathbb{R} = 1 | B = 0) = p\_2 \dots$$

Hence, the average motor response [for prior *P(B* = 1*)* = 1*/*2] is given by

$$P(R)\_{\text{motor}} = \sum\_{i} 1 \times P\_{\text{M}}(R=1|B=i)P(B=i) = \frac{1}{2} \left(p\_1 + p\_2\right) = p.c.$$

In the sensory state (during observation of the behavior), the reliability of a response quantified by the conditional probability of triggering a sensory response given presence of the behavioral feature in the stimulus is given by,

$$P\_S(R=1|B=1) = q\_1$$

and the probability of a sensory response without the behavioral feature (intrinsic noise) is:

$$P\_S(R=1|B=0) = q\_2.$$

The parameters *p*1*, p*2*, q*1, and *q*2 can be freely chosen in this minimal model, for example *q*2 = *p*2 if intrinsic noise in sensory and motor states are assumed to be equal.

The average response in the sensory state is given by

$$
\langle R \rangle\_{\text{sensor}} = \frac{1}{2} \left( q\_1 + q\_2 \right) = q.
$$

The correlation between sensory and motor responses in this cell is

$$
\left\langle R\_{\text{motor}} R\_{\text{sensor}} \right\rangle = \sum\_{i} 1 \times P\_M(R=1|B=i) P\_S(R=1|B=i)
$$

$$
\times P(B=i) = \frac{1}{2} \left( p\_1 q\_1 + p\_2 q\_2 \right).
$$

And, the correlation coefficient between motor- and sensoryevoked response is

$$\text{CorrCoeff} = \frac{\langle R\_{\text{motor}} R\_{\text{sensor}} \rangle - \langle R\_{\text{motor}} \rangle \langle R\_{\text{sensor}} \rangle}{\sqrt{\langle R\_{\text{sensor}} R\_{\text{sensor}} \rangle \langle R\_{\text{motor}} R\_{\text{motor}}}},$$

$$= \frac{\frac{1}{2} (p\_1 q\_1 + p\_2 q\_2) - pq}{\sqrt{p(1 - p)q(1 - q)}}.$$

We can discuss the following special cases:

• perfect sensory tuning (*q*1 = 1*, q*2 = 0, no instrinsic noise in sensory state):

$$\text{CorrCoeff} = \frac{p\_1 - p\_2}{2\sqrt{p(1-p)}}$$

• same tuning and same intrinsic noise in motor and in sensory states (*q*1 = *p*1*, q*2 = *p*2):

$$\text{CorrCoeff} = \frac{\left(p\_1 - p\_2\right)^2}{4p\left(1 - p\right)}$$

In summary, the strength of mirrored responses scales linearly or quadratically with the contrastive probability that neural responses are locked to the behavioral feature vs. spontaneously driven.

## *Mike Skocik1 and Alexay Kozhevnikov1,2\**

*1 Department of Physics, Pennsylvania State University, University Park, PA, USA*

*2 Department of Psychology, Pennsylvania State University, University Park, PA, USA*

#### *Edited by:*

*Ahmed El Hady, Max Planck Institute for Dynamics and Self Organization, Germany*

#### *Reviewed by:*

*Ahmed El Hady, Max Planck Institute for Dynamics and Self Organization, Germany Noah Cowan, Johns Hopkins University, USA*

#### *\*Correspondence:*

*Alexay Kozhevnikov, Department of Physics, Pennsylvania State University, University Park, PA, USA. Department of Psychology, Pennsylvania State University, University Park, PA, USA. e-mail<EMAIL_ADDRESS>
Studies of behavioral and neural responses to distorted auditory feedback (DAF) can help shed light on the neural mechanisms of animal vocalizations. We describe an apparatus for generating real-time acoustic feedback. The system can very rapidly detect acoustic features in a song and output acoustic signals if the detected features match the desired acoustic template. The system uses spectrogram-based detection of acoustic elements. It is low-cost and can be programmed for a variety of behavioral experiments requiring acoustic feedback or neural stimulation. We use the system to study the effects of acoustic feedback on birds' vocalizations and demonstrate that such an acoustic feedback can cause both immediate and long-term changes to birds' songs.

**Keywords: acoustic feedback, animal vocalizations, behavioral neuroscience, sensory feedback, real-time data processing**

## **INTRODUCTION**

Distorted auditory feedback (DAF) is used for assessing the effects of auditory input on vocal production. Presenting DAF and assessing its effects on the song and on the neural activity have been used in songbirds to study the mechanisms of song production and learning (Leonardo and Konishi, 1999; Sakata and Brainard, 2006; Andalman and Fee, 2009; Keller and Hahnloser, 2009; Tschida and Mooney, 2012). Human speech is sensitive to certain types of DAF (Lee, 1950; Houde and Jordan, 1998), and DAF is used to study speech mechanisms. It is often desirable to have real-time DAF, i.e., to rapidly (in a few milliseconds or faster) detect the occurrence of specific acoustic elements in vocalization and present an auditory stimulus once the target acoustic element is detected.

In this paper, we describe an automated system for real-time DAF and demonstrate its use to study both the immediate and the long-term effects of DAF on the song of Bengalese finches. The system uses open-source software and, therefore, is extremely flexible and customizable by the user. It has a significantly lower cost than commercial systems.

Songbirds use auditory feedback to learn to sing when they are young and to maintain their songs in adulthood (Konishi, 1965; Brainard and Doupe, 2000). Long-term exposure to DAF has been shown to cause song degradation in songbirds (Okanoya and Yamaguchi, 1997; Woolley and Rubel, 1997; Leonardo and Konishi, 1999). Some bird species' songs exhibit immediate sensitivity to acoustic input. For these birds, DAF can have an immediate effect on the timing and acoustic structure of the song (Sakata and Brainard, 2006). Analyzing the effects of DAF can yield new understanding of the neural organization of the song and the mechanisms of song learning (Brainard and Doupe, 2000). To study the questions about the effects of time-localized DAF on birdsong, it is important to be able to deliver DAF with high temporal precision in relation to vocalization. To do this, it is necessary to rapidly and reliably detect the specific acoustic elements of the bird's song and, after detection of the acoustic element, generate an acoustic output.

It is a challenging technical task for an acoustic feedback system to be real-time. Real-time performance is most easily achieved with analog systems (Cynx and Von Rad, 2001), but digital systems offer significant advantages in terms of convenience and flexibility. However, the advantages of a digital system are accompanied by the difficulties of making a digital system have small and constant processing delays. The system has to perform analog-to-digital conversion, fast analysis of the recently acquired data and digital-to-analog conversion, and these operations have to take place with reliable timing and concurrently with saving the acquired data. Custom-made DAF systems have been developed and used in behavioral studies (Leonardo and Konishi, 1999; Kao et al., 2005), but their real-time processing characteristics have not been reported. Oftentimes, custom-made systems have significant and not well-controlled delays, especially for systems based on PC's running Windows. Commercial systems for real-time acoustic processing are available but are expensive.

We developed a real-time DAF system based on a PC running Linux and the Real-Time eXperiment Interface (RTXI) software (Lin et al., 2010) and a National Instruments A/D card. The system is low-cost (the cost is only the cost of the hardware, the software is free). The system is capable of A/D bandwidth of over 30 kHz with real-time processing of acoustic signals.

## **METHODS**

A PC with an Intel i7 six-core processor (2.66 GHz) and 4 GB or RAM running Ubuntu Linux <IP_ADDRESS>-rtai with RTXI version 1.1.2 and a National Instruments PCIe-6251 A/D card is used. The A/D card receives audio input from a microphone (AudioTechnika PRO-44, used with Behringer Shark DSP110 microphone amplifier). The output is sent to a speaker amplifier (SLA-1, Applied Research and Technology); the output of the amplifier is connected to a speaker.

The Data Recorder software within RTXI is custom modified. The simplified diagram of signal processing is shown in **Figure 1**. The system has two modes of operation—a non-triggered (idle) mode and a triggered (active) mode. At the core of the modified software is the circular buffer that takes data points one-by-one from the data acquisition engine once they become available. In the non-triggered mode [**Figure 1 (top)**], the system continuously (every 1 ms) computes the rms of the last 10 ms of the input signal. If the signal rms exceeds the threshold, the system is switched into triggered mode.

**FIGURE 1 | Block diagram of the acoustic feedback system.** When not triggered **(top)**, the system computes the rms of the input signal. When the rms exceeds the threshold, the system is triggered. When triggered **(bottom)**, the system computes the spectrogram of the most recent 20 ms of signal and computes the correlation coefficient of this spectrogram with the spectrogram of the template sound (e.g., song syllable). The template sound is detected when the correlation coefficient exceeds a threshold value; in this case, acoustic feedback can be generated. Both the input and the acoustic output are saved to the computer hard drive.

In triggered mode, the system does real-time processing of auditory data. In **Figure 1 (bottom)**, we show the processing done for recognizing the song syllable of a Bengalese finch. The FFT of the past 256 data points (∼8.4 ms) is computed every 1 ms and stored in an FFT circular buffer. Every 1 ms, the spectrogram of the most recent 40 ms of the input signal is obtained from the FFT circular buffer. A correlation coefficient between the input signal spectrogram and the spectrogram of the template is computed. If the correlation coefficient exceeds the threshold, the system detects the occurrence of the target song syllable, and acoustic feedback can be generated, or further processing can be done. While generating acoustic feedback, the system keeps going through all of the above steps, but is disallowed from registering another detection to prevent it from triggering on its own output.

The presence of the data circular buffer allows very fast access to chunks of the most recent data for processing without affecting the timing of the data acquisition process. The FFT circular buffer also allows extremely fast computations of the spectrograms of the sound (computing the spectrogram is a computationally-intensive task). This enables the system to recognize complex vocal elements based on their spectrogram (e.g., frequency sweeps) without compromising the timing. While triggered, every 1 s, the system computes the rms of the previous 200 ms of the input signal to check if the acoustic input is still present. If the rms is below a threshold (no signal), the system goes into the idle mode. While triggered, the system continually saves all the data acquired in a separate array and saves the data to the hard drive once it is switched back to idle mode. A more detailed description of this system, along with the source code, is available at http://www.phys.psu.edu/∼akozhevn/ac\_feedback/.

## **RESULTS**

We tested the performance of our DAF system in several tasks which are often needed in behavioral experiments using acoustic feedback. We also used the system to assess the effects of acoustic feedback on the song of Bengalese finch. All animal procedures were carried out in accordance with the locally approved IACUC protocol.

## **DELAY BETWEEN INPUT AND OUTPUT TEST**

A simple task is generating acoustic feedback when the input level exceeds a certain threshold. Although this task may be too simple for most behavioral experiments, the delay in the system between detecting the crossing of the threshold and producing the output is a useful figure for indicating how fast the system can be when it is solely converting A/D and D/A and saving data without any complex data processing.

The system was programmed so that, once the input exceeded a fixed threshold, the acquired input signal was sent to the D/A output with no extra processing. A square wave with the amplitude exceeding the threshold was applied to the input; the delay between the input and the output was measured with the digital oscilloscope. The measured delay between the output and the input was 27 ± 9µs (mean ± SD, min = 9 µs, max = 43µs). The sample rate was 30.3 kHz, so the observed delays corresponded to a delay of 1 data point between the input and the output. The observed variations of the delay are due to the difference in timing between the external input signal and the timing of the A/D events. In all cases, however, the delay between the input and the output does not exceed 1 datapoint. Therefore, the system has real-time capability.

#### **DETECTION OF SPECIFIC VOCAL ELEMENTS IN THE BIRD'S SONG**

A typical task in experiments using acoustic feedback is detection of a certain "template" sound. The template can be either a sound of a certain frequency or a more complex combination of frequencies, frequency sweeps, etc. Once the template is detected, the acoustic feedback can be played back to the animal. This task is computationally intensive because one needs to compute the characteristics of the recently acquired input signal, then compare these with the characteristics of the template and, if the input is sufficiently similar to the template, decide that the detection has occurred and generate acoustic output. The computation has to be done fast enough to enable realtime performance and not interfere with the data acquisition process.

Common techniques that have been used for detecting acoustic elements are spectrogram-based techniques (Leonardo and Fee, 2005) and feature-based techniques (Tchernichovski et al., 2000). In a spectrogram-based approach, the spectrogram of the recently acquired signal is computed and compared to the template spectrogram. A common way to accomplish this is to compute the correlation coefficient between the two spectrograms. Detection of the template sound occurs if the correlation coefficient exceeds a threshold value.

We tested the performance of the system for detection of specific syllables in the song of a Bengalese finch. The Bengalese finch song consists of a sequence of syllables separated by silences (inter-syllable gaps) (**Figure 2**). The acoustic structure of the song syllables is fairly stable; the main source of variability from one song to another is the sequence of syllables in each song (Honda and Okanoya, 1999).

To detect a specific song syllable, the system continuously computes the correlation coefficient of the spectrogram of the most recent 20 ms segment of the acquired signal with the spectrogram of a 20-ms syllable template (see Methods, **Figure 1**). The target syllable is detected by the system when the correlation coefficient exceeds the threshold value of 0.8. The value of the threshold was chosen by examining the target syllable detections by the DAF system in a set of about 20 songs and comparing the detected syllable occurrences with the actual occurrences of the target syllables determined by visual examination of the song spectrograms. If the threshold is set too high, the probability of missing the target syllable increases. Setting the threshold too low increases the probability of false positive detections. After the syllable detection, acoustic feedback (either white noise or the song syllable) can be played back to the bird.

Typical performance of the system on the real-time syllable recognition task is shown in **Figure 2**. The top spectrogram shows "detection only" mode—the system detects the target syllable in real time, but no playback is generated. The bottom spectrogram shows detection and playback generation—after detecting the target syllable, the system plays back another song syllable

**FIGURE 2 | Top:** spectrogram of the song of a Bengalese finch and the times of occurrence of one of the song syllables. The system was programmed to only detect the occurrences of the target syllable in real time, no acoustic feedback was generated. The detection times are shown as vertical red lines. **Bottom:** the system is detecting the target syllables (vertical red lines) and is generating acoustic feedback after detection. The acoustic feedback waveform is shown below. The feedback signal is one of the birdsong syllables; the acoustic feedback pickup by the microphone is visible on the spectrogram. The zoomed-in spectrogram of the template is shown on the right.

to the bird. The vertical red lines indicate the detection times of the target syllable. The zoomed-in spectrogram of the template is shown on the right. The template contains part of the intersyllable interval and the first 20 ms of the target syllable, so the end of the template (detection time) is approximately in the middle of the 40-ms long syllable.

Performance of the system was checked by comparing the results of automatic detections of the system with the manual identification of the target song syllables carried out by off-line examination of the spectrograms. Out of 659 target syllables, 610 were correctly detected and 49 were missed. There were zero false positives. Thus, the system shows robust performance with the real-time syllable recognition task: over 92% of the target syllables were correctly identified.

This demonstrates that the system is capable of real-time detection of target syllables in the song. Note that the syllables occurring after the target syllables in **Figure 2** are frequency sweeps that overlap with the template's frequencies. The system discriminates them from the target syllables because they have a different frequency profile. Such discrimination is an advantage of the spectrogram-based detection; this would not be possible if only instantaneous frequencies were detected.

Additionally, we tested the system on the detection of syllables in the song of a zebra finch—another bird species. We used our dataset of zebra finch songs with known syllable sequences obtained in a previous study (Kozhevnikov and Fee, 2007). Zebra finch songs were played back through the speaker, and the results of the real-time detection by our DAF system were compared to the known occurrences of the target syllable.

A small subset of songs (10 songs) was used as a test set: the threshold value for the syllable detection was adjusted to optimize the percentage of correctly detected syllables in this small test set. After this, the threshold kept was fixed, and the performance of the system was tested on the whole dataset (about 100 songs). Out of 756 target syllables in the dataset, 728 were correctly detected, 28 were missed; there were 4 false positives. The system correctly detected over 96% of the target syllables in the dataset; the probability of a false positive detection was less than 1%.

#### **EFFECTS OF AUDITORY FEEDBACK ON THE TIMING OF THE BIRDSONG**

Auditory feedback has been shown to have immediate effects on some animal vocalizations. For Bengalese finches, DAF has been shown to affect the timing of song syllables. DAF played after the song syllable increases the time interval between that syllable and the next syllable in the song (Sakata and Brainard, 2006). We tested whether our feedback system is effective in causing real-time changes to the Bengalese finch song. The system was programmed to detect one of the song syllables and, once the syllable was detected, to play back the same song syllable with a probability of 0.05. This ensured that the feedback was sufficiently sparse so almost in all cases there was only one playback during each song. The feedback and control trials were randomly interleaved. This simplified the analysis of the syllable timing and eliminated any confounding effects from playbacks being too close to one another. The delay between the syllable sung by the bird and the syllable playback was 40 ms.

The playback causes some pickup on the input channel, which can cause difficulties in precise determination of the timing of the syllable that is occurring during the playback. Therefore, the time interval between the target syllable and the following syllable (which is partially overlapped with playback) was computed as one half of the difference between the detection time of the target syllable and the detection time of the second syllable after the target syllable. The same procedure was performed in control trials to ensure consistency in data analysis. Since the distributions of time intervals may not be Gaussian, we use a non-parametric statistical test—two-way Kolmogorov–Smirnov test—to assess the statistical significance of DAF effects on the song timing.

**Figure 3** shows the distributions of the time intervals between the target syllable and the following syllable when the feedback is present (blue histogram) and when there is no feedback (red histogram). The widths of the distributions are due to the natural variability of the song timing. In the presence of feedback, the time intervals between the syllables become longer. Without DAF, the mean interval is 74.8 ms (*N* = 637 syllables); in the presence of DAF, the mean interval is 75.7 ms (*N* = 97 syllables). Although the change of the mean duration is small compared to the widths of the distributions, the effect is highly statistically significant (*p* = 0*.*001, two-way Kolmogorov–Smirnov test). The observed lengthening of the time interval between the song syllables is consistent with previous observations (Sakata and Brainard, 2006). Thus, our acoustic feedback has an immediate effect on the song: DAF immediately and reversibly affects song timing.

#### **LONG-TERM EFFECTS OF DAF ON ACOUSTIC STRUCTURE OF THE SONG**

DAF has been shown to cause long-term changes to animal vocalizations. For songbirds, prolonged repeated presentation of DAF

**Bengalese finch song syllables.** Shown above are the histograms of the time intervals between two subsequent syllables in the song in the presence of DAF (blue) and without DAF (red). The means are: *t*mean = 74*.*8 ms (control, *N* = 637 syllables) and *t*mean = 75*.*7 ms (feedback, *N* = 97 syllables), the difference is statistically significant (*p* = 0*.*001, two-way Kolmogorov–Smirnov test).

can cause gradual change of the song (Leonardo and Konishi, 1999; Warren et al., 2011). We tested out system on the task of causing long-term changes of the frequency of one of the song syllables.

The system is programmed to detect the fundamental frequency of one of the song syllables. After the target song syllable is detected, the temporal profile of the pitch (defined as the largest peak in the FFT of the latest 256 points) was computed. The lowest value in the pitch profile in the time window between 3 and 12 ms after the detection time was taken to be the pitch of the syllable. The feedback (white noise) is conditional on the detected pitch of the song syllable. For example, the feedback can be generated if the detected syllable pitch is smaller than a threshold value. Continuous exposure to such feedback has been shown to cause the bird to gradually shift the mean pitch of the syllable so that the feedback is generated less often—the bird adapts its song to avoid hearing DAF (Warren et al., 2011).

We tested whether such conditional DAF could shift the mean syllable pitch in both directions. **Figure 4** shows the long-term effects of DAF which is conditional on the syllable pitch. During days 1–6, the feedback was played back if the pitch was less than 3530 Hz; during days 7–12, the feedback was played back if the pitch was greater than 3530 Hz; during days 13–18, the feedback was played back if the pitch was less than 3510 Hz.

Our feedback is effective in causing gradual changes in the syllable pitch. Playing back DAF when the frequency is lower than the threshold value causes an upward drift of the mean pitch (days 1–6 and 13–18). Presenting feedback when the pitch is higher than the threshold causes a downward drift in pitch (days 7–12). This shows that the system is suitable for studies of long-term effects of DAF on animal vocalizations.

**FIGURE 4 | Prolonged exposure of the bird to DAF causes gradual changes in the song.** The mean pitch of one of the target song syllable was manipulated by feedback conditional on pitch. Each datapoint is the pitch of the target syllable averaged over all the renditions of the target syllable sung by the bird on a given day. The number of renditions varies from day to day (mean = 208, min = 126, max = 288). Error bars are s.e.m. Dashed lines are the values of the threshold. Arrows indicate the direction in which the syllable frequency was expected to change in response to DAF. Larger error bars on day 12 are mostly due to the smallest number of target syllable renditions (*n* = 126) sung on that day.

### **DISCUSSION**

The described real-time acoustic feedback system is a versatile tool for studies of the effects of auditory feedback on ongoing animal vocalization. The advantage of the system is the flexibility of the processing than can be realized. The circular buffer allows real-time acquisition of arbitrary-length segments of the most recent data without affecting the timing of data acquisition. In addition, having a separate FFT buffer facilitates real-time spectral processing of acquired signals. This feature enables a very quick creation of the spectrograms of long (tens or even hundreds of milliseconds) segments of signals at a high rate (the spectrogram is updated every 1 ms).

This capability is very useful for the detection of complex vocal signals. Often, it is not just a certain frequency that needs to be detected, but rather a certain spectrogram pattern, like multiple frequencies or the frequency sweeps frequently seen in birdsong syllables. Since the same frequency can occur in many syllables, it is the whole pattern of the spectrogram that allows real-time detection of the syllable. Our system is very well-suited for rapid spectrogram-based detection of acoustic elements.

The performance of the system will vary depending on the type of animal vocalization and the nature of the acoustic element being detected for two main reasons. First, there is always a natural variability in the acoustic structure of a vocal element, and the degree of this variability may be different for different vocal elements; this will affect the reliability of detection. For example, a birdsong syllable can possess a more or less stereotyped spectrogram; the detection will be easier for a more stereotypical song syllable. Second, a given vocal element can be more or less similar in its acoustic structure to other vocal elements; reliable detection of a target vocal element will be easier if it is spectrally more dissimilar to other vocal elements. To achieve optimum performance, adjustments to the threshold or detection algorithm may be needed; thus, it is important to have a highly customizable system.

It is worth mentioning that, when the song syllables are detected, data processing is not a time-limiting step, and significantly more complex processing can be done without decreasing the A/D rate. We tested the system with longer templates (60 and 100 ms); they did not affect performance. We also tested the simultaneous detection of two templates, so that, every 1 ms, the system computed the correlation coefficient of the sound spectrogram with two template spectrograms, and that also did not affect the A/D rate. For a template 60 ms long, the computation of the correlation coefficient takes 16µs; this time scales linearly with the length of the template. The computation of the FFT (to fill the column in the spectrogram buffer) takes 7–8µs. FFT and correlation coefficient computations are the slowest signal processing steps; all other steps combined take less than 1µs. Thus—If the spectrogram update rate is kept at 1 ms—the system should be capable of simultaneously detecting of over 10 different song syllables. Therefore, fairly complex real-time analysis and detection of multiple vocal elements can be done without compromising the speed of the system.