A dynamic network model of temporal receptive fields in primary auditory cortex — TEI-tagged view

RahmanpersName, WillmorepersName, KingpersName, HarperpersName (2019date). “A dynamic network model of temporal receptive fields in primary auditory cortex.” PLoS Computational Biologytitle·j.

Object of study — why cortical auditory responses appear to depend on hundreds of milliseconds of sound history when the delays available to the circuit are far shorter

Research question — how much of the long temporal receptive field can be produced by neurons' own decaying memory rather than by long delay linesterm

Contribution — shows a network whose units each carry their own exponential decay matches a 200 ms delay-line model using only 25 ms of delay, which is within the biologically available range

Full text — tagged

The reference list is not part of the text we hold for this paper. Our copy was extracted from the publisher's page and the bibliography did not survive that extraction, so the article below ends before its references. Nothing has been removed by us.

A dynamic network model of temporal receptive fields in primary auditory cortex | PLOS Computational Biology Article Authors Metrics Comments Media Coverage Reader Comments Figures Figures Abstract Auditory neurons encode stimulus historyterm, which is often modelled using a span of time-delays in a spectro-temporal receptive field (STRFterm). We propose an alternative model for the encoding of stimulus historyterm, which we apply to extracellular recordings of neurons in the primary auditory cortex of anaesthetized ferrets. For a linear-non-linear STRFterm model (LN modelterm) to achieve a high level of performance in predicting single unit neural responses to natural sounds in the primary auditory cortex, we found that it is necessary to include time delays going back at least 200 ms in the past. This is an unrealistic time span for biological delay linesterm. We therefore asked how much of this dependence on stimulus historyterm can instead be explained by dynamical aspects of neurons. We constructed a neural-network model whose output is the weighted sum of units whose responses are determined by a dynamic firing-rate equationterm·def. The dynamic aspect performs low-pass filtering on each unit’s response, providing an exponentially decaying memory whose time constantterm is individual to each unit. We find that this dynamic network (DNet) modelterm·def, when fitted to the neural data using STRFsterm of only 25 ms duration, can achieve prediction performance on a held-out dataset comparable to the best performing LN modelterm with STRFsterm of 200 ms duration. These findings suggest that integration due to the membrane time constantsterm or other exponentially-decaying memory processes may underlie linear temporal receptive fields of neurons beyond 25 ms.

Author summary The responses of neurons in the primary auditory cortex depend on the recent history of sounds over seconds or less. Typically, this dependence on the past has been modelled by applying a wide span of time delays to the input, although this is likely to be biologically unrealistic. Real neurons integrate the history of their activity due to the dynamical properties of their cell membranes and other components. We show that a network with a realistically narrow span of delays and with units having dynamic characteristics like those found in neurons, succinctly models neural responses recorded from ferret primary auditory cortex. Because these integrative properties are widespread, our dynamic network provides a basis for modelling responses in other neural systems.

Citation: RahmanpersName M, WillmorepersName BDB, KingpersName AJ, HarperpersName NS (2019) A dynamic network model of temporal receptive fields in primary auditory cortex. PLoS Comput Biol 15(5): e1006618.

https://doi.org/10.1371/journal.pcbi.1006618 Editor: Frédéric E. Theunissen, University of California at Berkeley, UNITED STATES Received: November 4, 2018; Accepted: April 13, 2019; Published: May 6, 2019 Copyright: © 2019 RahmanpersName et al. This is an open access article distributed under the terms of the Creative Commons Attribution License , which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: The experimental data underlying the results presented in the study were previously published by HarperpersName et al. (2016) and are available from https://osf.io/ayw2p/ . All codes are available from https://github.com/monzilur/DNetterm .

Funding: BDBW, NSH, and AJK were supported by Wellcome Trust funding (WT108369/Z/2015/Z). MR was supported by the Clarendon Fund Scholarship. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Introduction The tuning properties of auditory neurons are commonly described by a spectro-temporal receptive field (STRF) modelterm·def, which characterizes the linear dependence of the neural response on the sound spectrum at a range of latencies [ 1 – 14 ]. A static non-linearity is often applied to the output from the linear STRFterm—this linear non-linear (LN) model estimates the responses of neurons significantly better than the linear estimate [ 15 , 16 ]. While STRFterm models are somewhat successful in explaining the dependence of neural responses on the past few hundred milliseconds of stimulus historyterm [ 17 ], these models do not show how this temporal aspect of receptive fields might be implemented biologically. The most direct biological interpretation of the STRFterm model would involve auditory cortical neurons receiving inputs at a range of simple delays spanning out to 200 ms; however, this is biologically unrealistic. The onset latenciesterm of neurons in the ventral division of the medial geniculate body (vMGBterm), which provides the primary ascending input to primary auditory cortex (A1term), are typically less than ~14 ms in the ferret [ 18 ], the animal species that provided the data that were modelled in this study. In the cat vMGBterm response latencies are typically less than ~20 ms [ 19 ], in the mouse less than ~10 ms [ 20 ], and in the marmoset less than ~40 ms [ 21 ]. These values are far less than the duration of stimulus historyterm which influences the responses of A1term neurons. Although a few vMGBterm units do show onset latenciesterm beyond this range (in cat up to 70 ms [ 19 ], in mouse up to 60 ms [ 20 ], and in marmoset up to 300 ms [ 21 ]), this small fraction is unlikely to be the main determinant of the dependence of the neural response in primary auditory cortex on stimulus historyterm beyond a few tens of milliseconds. This is because there are other plausible mechanisms, such as the dynamic temporally-integrative properties intrinsic to neurons.

In this study, using recordings from anesthetised ferrets, we asked how much of the dependence on recent stimulus historyterm—the temporal receptive fields, of primary auditory cortical neurons can instead be explained by certain of these simple dynamical aspects of neurons–an approach more consistent with the known biology. To model a neuron’s response, we constructed a neural network model whose output is the weighted sum of the responses of multiple units, each of which resembles an LN modelterm. However, the response of each unit of the network was modified in accordance with a dynamic firing-rate equationterm·def. This low-pass filters the unit’s response (by convolving with an exponential decay impulse response), providing a simple exponentially decaying memory. This integrative characteristic can be related to the capacitance and resistance of nerve cell membranes, and the individual time constantterm of each unit can be related to the membrane time constantterm of real neurons [ 22 ]. The dynamic aspect of our network can alternatively be interpreted in terms of other neurobiological dynamical phenomena, such as channel-based neural adaptation, short term synaptic plasticityterm, or recurrent network properties. We found that our biologically-motivated dynamical model can accurately capture the temporal receptive fields of neurons without the need for a wide span of latenciesterm·var of input.

Materials and methods Ethics statement In vivo electrophysiological recordings in ferrets were performed under ketamine (5 mg/kg/h) and medetomidine (0.022 mg/kg/h) anaesthesia. All animal procedures were performed under license from the UK Home Office and were approved by the University of Oxford Committee on Animal Care and Ethical Review.

Stimuli Models were fitted to single-unit neural responses of anesthetized ferrets to natural sound stimuli. Altogether, 20 sound clips were presented containing human speech in different languages, other animal vocalizations (e.g. ferret) and environmental sounds (e.g. wind and water). All clips were 5 s long with a sampling rate of 48,828.125 Hz and with root mean square intensity ranging from 75 to 82 dB SPL.

Experimental setup and neural responses The electrophysiological data used in this study were taken from a series of experiments used in a previous study by HarperpersName et al. [ 23 ]. Briefly, recordings were made from the primary auditory cortical areas, A1term and the anterior auditory field (AAF) of 6 adult pigmented ferrets (5 females and 1 male) under ketamine (5 mg/kg/h) and medetomidine (0.022 mg/kg/h) anaesthesia for 20 repeats of the 20 natural sound clips played in random order. All animal procedures were performed under license from the United Kingdom Home Office and were approved by the local ethical review committee. In total, 56 penetrations resulted in 549 single and multi-units, of which 284 were single units. A single unit was taken for analysis only if its activity was driven by the stimulus according to the noise ratioterm (see below). For each unit, the number of spikes was counted in each 5 ms time bin and averaged over repeats, to provide a response profile y n ( t ), where t is time, and n is the clip number. The total number of time bins in a clip is T . For simpler notation, we will drop the subscript n , unless we note otherwise.

Noise ratioterm The noise ratioterm was calculated as a measure of how much the response of each unit is dependent on the stimuli [ 9 , 24 ]. The noise ratioterm was measured over all 20 stimuli and 20 repeats. Any unit with a noise ratioterm > 40 were excluded from the study. Also, only putative single units were used. This resulted in 73 neurons for the study.

Cochleagramsterm The models receive the input as a cochleagramterm, a spectrogram-like transformation of the sound waveform. The cochleagramterm approximates the spectral filtering performed by the auditory periphery and was calculated as follows [ 23 , 25 ]. For each sound clip, the amplitude spectrum was measured using 10-ms Hanning windows, overlapping by 5ms. The number of frequency channels was then reduced by weighted summation using overlapping triangular windows to provide 34 log-spaced channels (500 Hz to 22,627 Hz center frequencies, adapted from melbank.m ( http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html )). Next, a log function was applied to the value in each time-frequency bin, and any values below a low threshold was set to the threshold. The cochleagramterm was normalized to zero mean and unit variance over the training set (see Cross-validationterm·stat and model testing). The input x fq ( t ) at time t was the chunk of cochleagramterm preceding t , given by x fq ( t ) = K f ( t − q +1), where f is the frequency channel and q is the latency and K f ( t ) is the cochleagramterm.

The LN modelterm The two stage LN modelterm consists of a linear model followed by a sigmoid output nonlinearity.

Linear stage : (1) where a(t) is the model neuron's activation, w fq is the input weight and b is the baseline activation. Both w and b are free parameters of the model and were estimated by regressing neural response y(t) against input x ( t ) using glmnetterm·stat [ 26 ]. To overcome overfitting, the weights were regularized using an L1-norm (LASSO regularizationterm). Thus a(t) can be seen as the best linear estimate of y ( t ) from x ( t ). The regularization hyperparameter λ controlled the strength of the regularization. To select λ , the model was subjected to k-fold cross-validationterm·stat (see Cross-validationterm·stat and model testing).

Nonlinear stage : The nonlinear stage of the model was a parameterized sigmoid nonlinearity: (2) The four parameters ρ i of the function were fitted by minimizing the squared error between v(t) and y(t) .

The network receptive field (NRF) modelterm·def The NRF modelterm is an artificial neural network with a single hidden layer of 20 hidden unitsterm (HU) which converge onto an output unit (OU). Each unit of the model is similar to a single LN modelterm. The activation of the j -th HU, a j ( t ) is, (3) The output of the j -th HU, v j ( t ) is, (4) where g( a j ( t )) is a sigmoid non-linear activation function, 1/(1+exp(− a j ( t ))). The HUs feed these outputs to the OU. The activation function of the OU a o ( t ) is, (5) The output v o ( t ) of the OU is, (6) which is the model’s estimate of the neuron’s response at time t .

The dynamic network (DNet) modelterm·def A neuron’s response depends on the activity of its presynaptic neurons through the total synaptic current they evoke. However, a neuron’s membrane potential, and consequently its firing rate, is not an instantaneous function of this current. Due to the capacitance and resistance of the neural cell membrane, the membrane potential and hence the firing rate is approximately a low-pass filtered version of the synaptic current. A common simple dynamic model of neurons that captures this phenomenon is the differential equation ( Eq 7 ) known as the firing-rate equationterm·def [ 22 ].

(7) Here, τ cont is the membrane time constantterm, v is post-synaptic firing rate, I ( t cont ) is the synaptic current at time t cont ,F() is a nonlinear function that relates synaptic current to the post-synaptic firing rate, and t cont is time. Here we use t cont and τ cont with the superscript ‘cont’ to denote that for this equation these variables operate in continuous time rather than discrete time. τ cont determines how rapidly firing rate approaches steady state for a constant synaptic current, and can be related to the membrane time constantterm [ 22 ]. Although the relationship is not strict, we shall call τ cont the membrane time constantterm. This is equivalent to convolving the instantaneous output of function F() with an exponential decay impulse response.

Eq 7 can be discretized by the Euler method, where c = 1 time-bin = 5 ms. We redefine v and I to depend on time-bin number, t , rather than continuous time ( t = t cont / c , where t is an integer). Hence, the discretized form of the firing rate equationterm·def is: (8) The dynamic network model incorporates units that behave like Eq 8 . In recognition that the characteristics of units in our model may not exactly correspond to biophysical properties of neurons described by the firing rate equationterm·def, and to be consistent with the NRF modelterm, we replace F() with g() and the current I ( t −1) with activation a ( t ).

As a result, activation a j ( t ) of the j -th HU of the DNetterm model is the same as the NRF modelterm, (9) However, the output v j ( t ) is different from that of the NRF modelterm, incorporating the discretized firing-rate equationterm·def.

(10) Similarly, the activation a o ( t ) of the OU is, (11) and the output v o ( t ) of the OU is, (12) Here, τ j and τ o are the membrane time constantsterm of j -th hidden unitterm and output unit respectively.

The sDNet modelterm·def The firing-rate equationterm·def ( Eq 7 ) assumes that the time constantterm that governs the relationship between the current and the firing rate is substantially larger than the time constantterm that governs the decay of post-synaptic potentials. However, if we assume the reverse, i.e., that the synaptic time constantterm is substantially larger, a different simple model of the relationship between the firing rate and the current becomes appropriate [ 22 ]. This model is given by: (13) (14) Here, I is synaptic current and v is post-synaptic firing rate, z ( t cont ) is the immediate influence of the presynaptic spikes on the synaptic current, F() is the nonlinear function of synaptic current to post-synaptic firing rate, is synaptic time constantterm, and t cont is time.

Similar discretization, like that used for the DNetterm model (see above), can be used for Eq 13 to use this for the units of a neural network model. Again, replacing the synaptic current, I , with the activation a, the activation of hidden unitsterm of the sDNet modelterm·def can be modelled as the following: (15) And the hidden unitterm output, replacing F() with g(), is given by: (16) Hence, in the sDNet modelterm·def, the activation of each HU, rather than the output of the nonlinearity, is convolved with an exponential decay impulse response.

Similarly, for the output units, (17) And the output, (18) Here, τ sj and τ so are the synaptic time constantsterm of j -th hidden unitterm and output unit respectively.

Objective function of the NRF and DNetterm models The free parameters w jfq , w j , b j , and b o (also τ j and τ o for the DNetterm and τ sj and τ so for the sDNetterm) were optimized by minimizing the squared error between v o ( t ) and y(t) subject to L1-regularization of the weights. Thus, the objective function is given by, (19) Here, n is included to indicate clip number, but was left out of other equations for simplicity. N is the number of clips used in training.

Parameter estimation All models except for the LN modelterm are fitted by minimizing the objective function with respect to the free parameters using the sum-of-function optimizer algorithm [ 27 ]. In using this algorithm, we take one clip to be one minibatch. The optimization algorithm requires calculation of the error gradients in respect to each of the parameters. For the NRF modelterm, error gradients are calculated using standard chain rule (backpropagation). For the dynamic models, the process is similar (see below).

For the network models, before training, the weights were initialized by modified Glorot initialization from a uniform distribution ranging from— to + , where fq is the number of input weights to a HU and l = 1 is the number of output weights from a HU. The biases were initialized similarly [ 28 ].

For the DNetterm model, to prevent the term running into mathematical error (when τ = 0), we substitute with . Parameters, d j and d o (for hidden and output units) were each independently initialized with the square root of a random variable from an exponential distribution with a mean of 1. Similar, substitution is used for the sDNet modelterm·def.

Error-gradient calculation for the DNetterm model To estimate the parameters, the gradients of the error with respect to various parameters are calculated using a method similar to the one used for real time recurrent learning [ 29 ].

Substituting , the dynamic equations for HUs ( Eq 10 ) and OU ( Eq 12 ) of the DNetterm model become respectively: (20) (21) Using the chain rule, the gradient of the error term, E( t ) at time t in respect to the weights of the output layer, w j , (22) (23) The prime denotes the derivative of the function. Similarly, the gradient of the error term, E( t ) in respect to the bias of the output layer, b o , (24) (25) Finally, the gradient of the error term, E( t ) with respect to the parameter, d o , (26) (27) Hence, the gradient is passed forward from the current time step to the next time step to be used in calculation of the gradient at that new time step. At t = 1, the values of and are undetermined. So, at time t = 1, the values of these terms are set to zero.

The gradients of the error term, E( t ), with respect to the rest of the parameters of the DNetterm model and the sDNet modelterm·def are obtained similarly.

Cross-validationterm·stat and model testing From the data for the 20 sound stimuli, 4 were chosen as a test set which was not used during training and cross-validationterm·stat. The cross-validationterm·stat set (the remaining 16 stimuli) was used to fit the models using k -fold cross validationterm·stat, where k = 8. The cross-validationterm·stat set was randomly divided into a training set of 14 stimuli and a validation set of 2 stimuli. The model was trained on the training set for 18 different values of the hyperparameter λ. A log spaced range of lambda values was used, but with a somewhat lower density at the extremes. For the LN modelterm, the exact values of λ used were: 1.00 x 10 −1 , 2.00 x 10 −2 , 1.17 x 10 −2 , 6.84 x 10 −3 , 4.00 x 10 −3 , 2.34 x 10 −3 , 1.37 x 10 −3 , 8.00 x 10 −4 , 4.68 x 10 −4 , 2.74 x 10 −4 , 1.60 x 10 −4 , 9.36 x 10 −5 , 5.41 x 10 −5 , 3.20 x 10 −5 , 6.40 x 10 −6 , 1.28 x 10 −6 , 2.56 x 10 −7 , and 5.12 x 10 −8 . For the rest of the models, the values of λ used were: 1.00 x 10 −3 , 2.00 x 10 −4 , 1.17 x 10 −4 , 6.84 x 10 −5 , 4.00 x 10 −5 , 2.34 x 10 −5 , 1.37 x 10 −5 , 8.00 x 10 −6 , 4.68 x 10 −6 , 2.74 x 10 −6 , 1.60 x 10 −6 , 9.36 x 10 −7 , 5.41 x 10 −7 , 3.20 x 10 −7 , 6.40 x 10 −8 , 1.28 x 10 −8 , 2.56 x 10 −9 , and 5.12 x 10 −10 . For each of the fitted models, neural responses were then predicted for the validation set, and the correlation coefficient between the actual neural responses and the prediction was measuredterm·stat. This process was repeated 8 times for different non-overlapping validation sets. The model was then retrained with the whole cross-validationterm·stat set using the λ value that provided the highest mean correlation coefficientterm·stat over all 8 folds. Next, the retrained network was used to predict the neural responses to the test set. All the correlation coefficientsterm·stat and normalized correlation coefficientsterm·stat shown are for this held out test setterm·stat, and all the model parameters shown are for the retrained network.

To ensure that all models receive the same amount of training dataterm·var, models with different latency spans of STRFsterm are provided with exactly the same durations of neural response. This was done by clipping the necessary number of time points from the beginning of the neural response data.

Effective hidden unitsterm·def The hidden layer of each of the networks contained 20 HUs, but because their weights were L1-regularized, they tended to develop substantive weights only if they explained aspects of the neural response that were not explained by any other units. As a result, many HUs developed weights close to zero, and thus the network models tended to have only a small number of effective HUs. The ‘effectivenessterm’ of each unit j was calculated [ 23 ] as the variance over time of the unit’s weighted output w o v j ( t ). Any HU with variance greater than 5% of the sum of the variances of all 20 HUs is considered effective.

IE scoreterm The IE scoreterm [ 23 ] measures whether the degree to which a HU is excitatory or inhibitory, on a scale between -1 and 1. Most inhibitory is IE = -1, when all of the input weights are negative or zero, and the output weight is positive, or when all of the input weights are positive or zero and the output weight is negative. Most excitatory is IE = 1, when all of the input weights are positive or zero, and the output weight is positive, or when all of the input weights are negative or zero and the output weight is negative.

(28) Membrane time constantterm Because the membrane time constantterm is modelled as 1+ d 2 and the duration of each time step is 5 ms, the value of membrane time constantterm in ms is 5 (1+ d 2 ).

Knocking out large time constantsterm of the DNetterm model After fitting a network model, HUs with time constantterm greater than the length of the STRFterm are selected. The output weights connecting these units to the output unit are then set to zero. As a result, HUs with large time constantsterm become inactive during the test process. Everything else in the model remained unchanged.

Results The neural response dataset used in this study was recorded from the primary auditory cortical areas, A1term and AAF, of anesthetized ferrets in response to a diverse selection of natural sounds (20 stimuli of 5 s duration), including speech, animal vocalizations and environmental sounds [ 23 , 25 ]. In total, 73 single-unit neural recordings showed sensitivity to the sound stimuli (see Methods ), as measured by their noise ratiosterm, and these were analysed in this study.

To probe the computations underlying stimulus encoding, models were fitted to estimate the neural firing rate (averaged over stimulus repeats) of each neuron as a function of the preceding sound stimulation ( Fig 1 ). First, for all models, the sound stimuli were pre-processed to generate a spectrogram-like representation (a cochleagramterm) that approximates the processing that occurs in the cochlea and auditory nerve. Second, using responses to 16 of the natural stimuli, a parameterized model was trained to estimate the firing rate as a function of the preceding cochleagramterm. After parameters of the models were fitted to the data, their fit quality was tested on a held-out test setterm·stat composed of the responses to the remaining 4 natural stimuli.

Download: PNG larger image TIFF original image Fig 1. The models have two-stages.

A. The first stage involves pre-processing that mimics the auditory periphery, where the sound waveform is converted into a time-frequency representation (cochleagramterm). B. The output of the first stage is then used as input to a model that is fitted to the neural responses of auditory cortical neurons. C. After parameter estimation, the model was used to predict neural responses to arbitrary stimuli.

https://doi.org/10.1371/journal.pcbi.1006618.g001 Using standard models, access to 200 ms of recent stimulus historyterm is needed to achieve best prediction For each neuron, an LN modelterm was fitted to estimate its firing rate as a function of the preceding cochleagramterm. The LN modelterm consists of an STRFterm, followed by a sigmoid nonlinearity ( Fig 2A ). As mentioned in the Methods, the STRFterm is the weighted sum of the cochleagramterm over frequency components at a span of latenciesterm·var. The effect of different spans of latencies (lengths) of the STRFterm was investigated. For the held-out test setterm·stat, and averaged over the 73 neurons, the mean normalized correlation coefficientterm·stat, mean CC normterm·stat [ 30 , 31 ], between the actual neural responses and the prediction of LN modelsterm with 25, 50, 100, 200, and 400 ms long STRFsterm are respectively 0.51, 0.64, 0.69, 0.71, and 0.71 ( Fig 2B ), where 1.0 is the maximum achievable prediction given the estimated neuronal noise and 0.0 indicates that there is no correlation between actual and predicted neural responses. The means of non-normalized Pearson correlation coefficients (CC)term·stat are 0.40, 0.50, 0.54, 0.55, and 0.55 respectively. Fig 2C shows STRFsterm of some example neurons estimated by the LN modelterm.

Download: PNG larger image TIFF original image Fig 2. Prediction of neural responses by the linear non-linear (LN) model.

A. The LN modelterm, which consists of a spectrotemporal receptive field (STRFterm) followed by a sigmoid nonlinearity. B. Normalized correlation coefficientsterm·stat (CC normterm·stat ) between neural responses and the prediction of the LN modelterm with different latency spans. Each black dot represents a neuron, the curve indicates the mean over 73 neurons and error bars indicate standard error of the meanterm·stat. C. STRFsterm of some example neurons, for LN modelsterm with different latency spans. Red, excitatory fields, blue, inhibitory. Excitatory and inhibitory fields extend up to 100 ms in most cases and up to 200 ms in some cases.

https://doi.org/10.1371/journal.pcbi.1006618.g002 A network receptive field model (the NRF modelterm) was also developed and fitted to each neuron. This model consisted of the weighted sum of multiple LN-like units. This model is the same as a previously published version [ 23 ], except that it uses sigmoid nonlinearities instead of scaled tanh non-linearities. A sigmoid non-linearity was used to be consistent with the dynamic network model that we will come to, which requires non-negative outputs and hence also uses sigmoid non-linearities. The NRF modelterm we use here is essentially an artificial neural network with a single hidden layer, using a sigmoid non-linearity in its hidden unitsterm and its single output unit ( Fig 3A ). The NRF modelterm makes use of multiple component STRFsterm, which allows modelling of neuronal sensitivity to non-linear interactions of multiple features [ 23 ]. STRFsterm of at least 400 ms are needed for best performance in predicting neural responses ( Fig 3B ). For the test set, the mean CC normterm·stat between neural responses and the model predictions made by NRF modelsterm with 25, 50, 100, 200 and 400 ms long component STRFsterm are 0.51, 0.64, 0.70, 0.72 and 0.73, respectively (mean CC 0.40, 0.50, 0.54, 0.55 and 0.56, respectively). The NRF modelterm predicts the neural responses equal to or better than the LN modelterm over all durations examined, although by a more modest amount than seen in an earlier study [ 23 ]. This may be because the NRF modelterm in this study uses a sigmoid rather than tanh nonlinearity, or because the earlier study [ 23 ] used a different way of selecting the test set. Fig 3C shows the input weights (analogous to STRFsterm) for the ‘effective’ HUs (see Methods ) of the model-fit for one example neuron, obtained using the NRF modelterm. Fig 3D and Fig 3E show the ‘effectivenessterm’ (see Methods ) and IE scoreterm (see Methods ) of the units shown in Fig 3C .

Download: PNG larger image TIFF original image Fig 3. Prediction of neural responses by the network receptive field (NRF) modelterm·def.

A. The NRF modelterm, which involves a weighted sum of multiple LN-model like units, each with an STRFterm. B. Normalized correlation coefficientsterm·stat (CC normterm·stat ) between neural responses and the prediction of the NRF modelterm with different latency spans. Each black dot represents a neuron, the curve indicates the mean over 73 neurons and error bars indicate standard error of the meanterm·stat. C. Component STRFsterm of an example neuron (Neuron #1). Each column represents an NRF modelterm with different latency spans. D. ‘Effectivenessterm’ (see Methods ) and E. IE scoreterm (see Methods ) of the components STRFsterm shown in panel C.

https://doi.org/10.1371/journal.pcbi.1006618.g003 The DNetterm model can characterize neural responses using very short duration STRFsterm Although both the LN and NRF modelsterm can capture dependence on stimulus historyterm using a span of input latencies extending up to at least 200 ms, there is little biological evidence for such a wide span of latenciesterm·var in the auditory thalamic neurons that feed to auditory cortex (see Introduction ). With the aim of better understanding the biological underpinnings of stimulus historyterm dependence in auditory cortical responses, a model with a dynamic aspect was developed that integrates the output of its component units over time in a manner similar to those of real neurons (see Methods ). Each unit in the DNetterm model includes exponentially-decaying response integration, motivated by the integrative properties of real neuronal membranes ( Fig 4A ). This model is in essence the NRF modelterm but with dynamic units, each modelled by the dynamic firing-rate differential equation ( Eq 7 ) [ 22 ]. That is, each unit of the model applies an exponentially decaying impulse response to the output of their non-linear activation function. Each unit’s decay rate is governed by its individually fitted time constantterm, τ (see Methods ).

Download: PNG larger image TIFF original image Fig 4. Prediction of neural responses using a dynamic network (DNet) modelterm·def.

A. The DNetterm model. The cochleagramterm is passed through a set of linear-nonlinear filters, whose response is then integrated over time using an exponentially-decaying impulse response to produce the hidden unitsterm’ outputs. A weighted sum of these outputs is passed through a similar output unit to make a prediction of the neural response. B. CC normterm·stat for the DNetterm model with different latency spans. Each black dot is a neuron, the curve indicates the mean over 73 neurons and error bars are standard error of the meanterm·stat. C. The input weights (STRFsterm) for the ‘effective’ hidden unitsterm (see Methods ) of 8 example neurons (each column is a different neuron) for the DNetterm model. D. ‘Effectivenessterm’ (see Methods ) and E. IE scoreterm (see Methods ) of the hidden unitsterm shown in panel C.

https://doi.org/10.1371/journal.pcbi.1006618.g004 We found that the mean CC normterm·stat between the actual neural responses and the predictions of the DNetterm model with 25, 50, 100, 200 and 400 ms STRFsterm are respectively 0.71, 0.71, 0.70, 0.68 and 0.68 ( Fig 4B ) (mean CC 0.55, 0.55, 0.54, 0.52 and 0.52, respectively). Because the 25-ms model performs the best and the biological range of latencies are typically less than 30 ms[ 32 ], we focused on the 25-ms DNetterm model for further analysis. Fig 4C shows the input weights (analogous to STRFsterm) for the ‘effective’ hidden unitsterm (see Methods ) of the model fits for 8 example neurons, using a DNetterm model with 25-ms long STRFsterm. Fig 4D and 4E show the ‘effectivenessterm’ (see Methods ) and IE scoreterm (see Methods ) of the units shown in Fig 4C .

The 25-ms DNetterm model out performs the DNetterm with longer (200 ms and 400 ms) STRFsterm, as measured by prediction of neural responses on a held-out test setterm·stat averaged over all 73 neurons ( Fig 5A ). It also achieves performance similar to the best performing LN modelterm, but performs a little worse than the best performing NRF modelterm ( Fig 5A ). This can be seen qualitatively in the time course of the prediction for an example neuron and stimulus in Fig 5B . A neuron by neuron comparison shows that this result is consistent for most neurons ( Fig 5C and 5D ). Also, the 25-ms LN and NRF modelsterm are outperformed by a large margin by the 25-ms DNetterm model ( Fig 5A ). A neuron by neuron comparison shows that this result is also consistent for almost all neurons ( Fig 5E and 5F ).

Download: PNG larger image TIFF original image Fig 5. Comparison of the performance of the DNetterm model with other models.

A. Mean CC normterm·stat of linear-nonlinear (LN), network receptive field (NRF) and dynamic network (DNet) modelsterm·def. B. Time course of predictions by best performing versions of each of the LN, NRF and DNetterm models, compared with real neural responses for an example neuron and stimulus. C. Neuron by neuron comparison of CC normterm·stat values between the 25-ms DNetterm model and 200-ms LN modelterm. Each red dot is a neuron. D. Neuron by neuron comparison of CC normterm·stat values between the 25-ms DNetterm model and 400-ms NRF modelterm. E. Neuron by neuron comparison of CC normterm·stat values between the 25-ms DNetterm model and 25-ms LN modelterm. F. Neuron by neuron comparison of CC normterm·stat values between the 25-ms DNetterm model and 25-ms NRF modelterm.

https://doi.org/10.1371/journal.pcbi.1006618.g005 Note that the DNetterm model is equivalent to the NRF modelterm when the discretized time constantterm for every DNetterm model unit is at the smallest possible value ( τ = 1, d = 0). Hence, because all possible NRF modelsterm are a subset of all possible DNetterm models, the DNetterm should in theory always equal or outperform the NRF modelterm. The fact that this is not the case suggests that, at least in some cases, the very best fits are not being found for the DNetterm model, perhaps as a consequence of some overfitting. That the DNetterm model performs worse for large latency spans, where there are more variables, is again probably due to challenges in fitting the DNetterm model, likely due to some overfitting to the training set.

Network structure and response integration are both required for an effective DNetterm To understand how much of an effect HUs of the DNetterm model with long time constantsterm have on the model’s predictions, the fitted model was modified by ‘knocking out’ (see Methods ) HUs with time constantsterm longer than the duration of their STRFsterm. This was achieved by setting the output weights of those units to zero, while keeping all other parameters intact. We found that the 25-ms DNetterm model is the DNetterm whose CC normterm·stat is most reduced by this alteration ( Fig 6A ).

Download: PNG larger image TIFF original image Fig 6. Multiple time constantsterm in a network architecture are required to capture the dependence of neural response on stimulus historyterm.

A. DNetterm model performance when large time constantsterm are knocked out. B. Performance of the single-unit dynamic model. C. Number of effective hidden unitsterm·def (HUs) used to model each of the 73 neurons using the DNetterm model. D. Distribution of time constantsterm of HUs of the DNetterm model. E. Distribution of time constantsterm of the output units (OUs) of the DNetterm model. F. Relationship between time constantterm and HU effectivenessterm for the DNetterm model.

https://doi.org/10.1371/journal.pcbi.1006618.g006 A single-unit dynamic model was also developed to examine the effect of a time constantterm in an LN-like single STRFterm model. Prediction performance of this model was found to be worse than the LN modelterm ( Fig 6B ). This indicates that both long time constantsterm and a network architecture are needed for the 25-ms DNetterm model to achieve prediction performance similar to that of standard models.

Although the DNetterm model had 20 hidden unitsterm, the L1 regularization eliminated unnecessary units. Hence there were fewer effective HUs, and for the 25-ms DNetterm model this number was found to be between 2 and 8 ( Fig 6C ). We defined an effective HU as one whose effectivenessterm (variance of weighted output, see Methods ) exceeded a certain threshold.

In the 25-ms DNetterm model, the time constantsterm of the effective HUs ranged between 5 and 475 ms, often far greater than the duration of the STRFsterm ( Fig 6D ). The time constantsterm of the OUs of the DNetterm model ranged from 5 to 130 ms ( Fig 6E ). To investigate whether long time constantsterm are associated with effective HUs, for each neuron HUs were ranked in order of effectivenessterm and the mean of the time constantsterm of all HUs in the same rank (over all neurons) were calculated. This showed that high-ranked HUs tend to have long time constantsterm ( Fig 6F ).

Hidden unitsterm with long time constantsterm tend to be more inhibitory To test whether a relationship exists between the time constantterm and unit’s inhibitory/excitatory nature, we classified the effective HUs of the 25-ms DNetterm model into excitatory or inhibitory using an IE scoreterm (see Methods ). If the output weight is positive, an IE scoreterm of 1 means a unit with entirely excitatory influence and an IE scoreterm of -1 means a unit with entirely inhibitory influence. We found that almost all units having large time constantsterm are inhibitory in nature ( Fig 7 ).

Download: PNG larger image TIFF original image Fig 7. Relationship between inhibitory/excitatory (IE) score and time constantterm.

Almost all the HUs of the DNetterm model with time constantsterm greater than the duration of the STRFsterm (units right of the vertical dotted line, the duration of the 25 ms STRFsterm) are inhibitory in nature. IE scoreterm = inhibitory/excitatory score.

https://doi.org/10.1371/journal.pcbi.1006618.g007 Membrane time constantterm versus synaptic time constantterm In the DNetterm model, the integration governed by the time constantterm occurs after each unit’s nonlinearity. Although the time constantterm is conceived and implemented in a manner consistent with a membrane constant, within the network it may capture other dynamic neural phenomena such as synaptic time constantsterm, forms of a channel-based or synaptic adaptation, or certain forms of network recurrency. An alternative network model places the integration just before the nonlinearity, which is more consistent with the time constantterm being synaptic in nature [ 22 ]. To investigate how much of the neural activity can be explained by such a model, we modified the DNetterm model to form the sDNet modelterm·def (synaptic DNetterm model), with each unit governed by Eqs 13 and 14 [ 22 ].

The sDNet modelterm·def achieves best performance in predicting neural responses when the duration of the STRFsterm is 200 ms (CC normterm·stat = 0.71) ( Fig 8A ). Although the sDNetterm with 25-ms STRFsterm predicts neural responses better (CC normterm·stat = 0.67) than the 25-ms LN modelterm (CC normterm·stat = 0.51) and the 25-ms NRF modelterm (CC normterm·stat = 0.51), it is worse than the 25-ms DNetterm model (CC normterm·stat = 0.71). This suggests that a large percentage of the variance of neural responses can be explained by a time constantterm before non-linearity, but that having a time constantterm after the non-linearity further improves prediction.

Download: PNG larger image TIFF original image Fig 8. A modified DNetterm model with synaptic time constantsterm (sDNet modelterm·def).

A. Performance of the sDNet modelterm·def. B. Distribution of the values of time constantsterm of hidden unitsterm (HUs) in the sDNetterm and DNetterm model C. Distribution of the values of time constantsterm of the output units (OUs) of the sDNetterm and DNetterm model. D. Relationship between a HU’s time constantterm and its effectivenessterm for the sDNet modelterm·def. E. Relationship between a HU’s inhibitory/excitatory (IE) score and its time constantterm for the sDNet modelterm·def. Most of the HUs with time constantsterm greater than the duration of the STRFsterm (units right of the vertical dotted line) are inhibitory in nature.

https://doi.org/10.1371/journal.pcbi.1006618.g008 While the largest time constantsterm of the effective HUs of the sDNet modelterm·def can be as large as ~470 ms ( Fig 8B ), the distribution of time constantsterm of output units have smaller values compared to the DNetterm model ( Fig 8C ). As with the DNetterm model the more effective HUs of the sDNet modelterm·def tend to have longer time constantsterm ( Fig 8D ) and units with long time constantsterm tend to be inhibitory in nature ( Fig 8E ).

Discussion We find that a model of primary auditory cortical neurons which incorporates exponentially-decaying memory processes–the DNetterm model–can capture some of the temporal aspects of the receptive field using very short (25 ms) component STRFsterm in a network. Standard models such as the LN modelterm (and the NRF modelterm) represent temporal receptive fields of neurons using delay linesterm whose latencies extend up to a few hundred milliseconds. Our results suggest that, beyond 25ms, the temporal receptive field can be modelled more succinctly by exponentially-decaying memory processes than by delay linesterm.

When auditory neurons are fitted using a standard LN modelterm, the resulting STRFsterm often contain narrowband inhibitory fields that decay over approximately 200 ms. In contrast, we find that the STRFsterm of the DNetterm model have short inhibitory fields with long time constantsterm. This suggests that these temporally extended inhibitory fields in STRFsterm may reflect neural mechanisms that can be approximated by exponentially-decaying memory processes.

Although they have not been described in the ferret, a few vMGBterm units do show onset latenciesterm beyond a few tens of milliseconds: in cat latencies up to 70 ms have been described [ 19 ], in mouse up to 60 ms [ 20 ], and in marmoset up to 300 ms [ 21 ]. Areas other than vMGBterm that project to cortex may also include long latency responses. Hence, we cannot rule out the possibility that long-latency inputs at least partially explain the dependence of auditory cortical neurons on stimulus historyterm–indeed it is likely that they do to some extent. However, these long-latency responses are relatively rare and we show that it is possible to explain a lot of the dependence of the neural response in primary auditory cortex by exponentially-decaying memory properties consistent with those properties intrinsic to neurons, such as membrane dynamics.

There are many biological mechanisms which may underlie exponentially-decaying memory processes. Neurons have various intrinsic dynamic integrative and adaptive characteristics originating from basic membrane dynamics, channel-based adaptation, post-synaptic potentials and short-term synaptic plasticity [ 22 ]. Recurrent network effects could also underlie such memory processes. In a previous study, a model with synaptic depression has been shown to partially capture the integration of stimulus historyterm in A1term [ 33 ]. Here, we found that a better prediction of neural responses can be achieved by a model based on the integrative properties of neural membranes (the DNetterm model) than by the one based on the dynamics of synaptic potentials (the sDNet modelterm·def). We note, however, that although the DNetterm model was developed with membrane dynamics in mind, aspects of all the above listed dynamic processes could be captured by the same exponentially-decaying memory in the DNetterm. These different processes are not exclusive of each other, and combinations of them may contribute to the effects we find.

The firing rate of A1term neurons has been demonstrated to depend on stimulus historyterm up to ~4 s in the past [ 17 ]. This dependence cannot be fully captured by linear models such as STRFsterm (17). This indicates that there are some nonlinear dynamic processes that provide history dependence which cannot be captured either by a linear or LN modelsterm, nor by the DNetterm model. Perhaps this is because the DNetterm model has limited capacity to capture temporal dynamics that are not exponential in nature.

Much of the relationship between stimulus and neural responses results from neurons being embedded in a network structure and from the nonlinear properties of the units in this network. From the auditory periphery to the A1term, there are numerous feedforward, feedback and lateral connections. Many previous models of auditory cortical neurons have not explicitly represented this network character. However, some recent models have incorporated non-linear contextual effects [ 14 , 25 , 34 – 36 ]. Furthermore, a few models do explicitly incorporate multiple receptive fields [ 37 – 43 ] or more extensive network structure [ 23 ], which perform better in predicting neural responses of A1term/AAF neurons than STRFterm models [ 23 , 37 ].

When modelling the dependence of neural responses on stimulus historyterm, one approach is to include the history of a neuron’s spikes and those of its recorded neighbours in a generalized linear model (GLM) model [ 35 , 44 ]. However, unlike the DNetterm model, this approach does not model multi-layer network structure with hidden unitsterm that are inferred but not recorded. Modelling hidden network structure has been applied to the salamander retina, where in vitro responses to simple artificial stimuli are well predicted by a network model (LNFDSNF) with various non-linearities and response delays mimicking the retinal pathway [ 45 ]. The LNFDSNF model is particularly relevant to our study because it is a network model whose units have some memory of their past activity, though the form of this memory differs from that of the DNetterm model. The LNFDSNF model units have a feedback kernel that acts before the nonlinearity—in contrast, the DNetterm model units integrate their activity after the nonlinearity using an infinite impulse response that decays exponentially.

The succinct and biologically-motivated description of neural responses provided by the DNetterm model provides a basis from which to develop further models to account for the ~30% of the CC normterm·stat that remains to be explained. One variant of the model that can be developed involves a network with multiple hidden layers where the units each have just a single latency (i.e. each unit in the first hidden layer uses just a spectral receptive field at a single delay). This approach provides an opportunity to build deeper neural-network-like models with fewer parameters, potentially overcoming overfitting problems for deep network models of neural responses [ 46 ]. This network can provide a basis for new network models that incorporate a greater range of intrinsic dynamic processes more explicitly, for example modelling in each unit its basic membrane dynamics, synaptic depression and channel-based adaptation together, with different appropriate time constantsterm in each case. Such models may show more predictive power.

Neurons are intrinsically dynamic systems, each with memory of past events. Furthermore, neurons are embedded in networks. We show that models incorporating both these characteristics can be valuable for explaining moment to moment in vivo responses of neurons in the auditory cortex. This is an arguably more biologically plausible way of capturing the history dependence of neural responses than a set of delay linesterm. In artificial neural networks, memory is usually implemented through recurrent connections between units. However, real biological neurons also have intrinsic memory processes with particular characteristics. Including such single-unit memory has potential advantages for explaining the in vivo responses of biological neurons and may also be useful in artificial neural networks. In summary, the dynamic network model we present here provides a biologically-plausible alternative encoding model to the classic LN modelterm, as well as a basis for further development of encoding models with hidden unitsterm with diverse intrinsic dynamic processes.