## Abstract

Single cells in the motor and somatosensory cortex of rats were stimulated *in vivo* with broadband fluctuating currents applied juxtacellularly. Unlike the DC current steps used previously, fluctuating stimulation currents reliably evoked spike trains with precise timing of individual spikes. Fluctuating currents resulted in strong cellular responses at stimulation frequencies beyond the inverse membrane time constant and the mean firing rate of the neuron. Neuronal firing was associated with high rates of information transmission, even for the high-frequency components of the stimulus. Such response characteristics were also revealed in additional experiments with sinusoidal juxtacellular stimulation. For selected cells, we could reproduce these statistics with compartmental models of varying complexity. We also developed a method to generate Gaussian stimuli that evoke spike trains with prescribed spike times (under the constraint of a certain rate and coefficient of variation) and exemplify its ability to achieve precise and reliable spiking in cortical neurons *in vivo*. Our results demonstrate a novel method for precise control of spike timing by juxtacellular stimulation, confirm and extend earlier conclusions from *ex vivo* work about the capacity of cortical neurons to generate precise discharges, and contribute to the understanding of the biophysics of information transfer of single neurons *in vivo* at high frequencies.

**SIGNIFICANCE STATEMENT** Nanostimulation of single identified neurons *in vivo* can control spike frequency parametrically and, surprisingly, can even bias the animal's behavioral response. Here, we extend this stimulation protocol to time-dependent broadband noise stimulation in sensory and motor cortices of rat. In response to such stimuli, we found increased temporal spike-time reliability. The information transmission properties reveal, both experimentally and theoretically, that the neurons support high-frequency stimulation beyond the inverse membrane time. Generating a stimulus using the neuron's response properties, we could evoke prescribed spike times with high precision. Our work helps to establish a novel method for precise temporal control of single-cell spiking and provides a simplified biophysical description of single-neuron spiking under time-dependent *in vivo*-like stimulation.

## Introduction

The spiking of single neurons in response to time-dependent signals is an important research topic for several reasons. Information transmission in the brain presumably does not rely exclusively on a summed population activity, but rather is also mediated by the precise spike times of individual cells (Strong et al., 1998; Zuo et al., 2015). Measuring and modeling the single-neuron response to time-dependent stimuli thus remains a crucial task in understanding the neural code (Rieke et al., 1996). Here, we explore three related issues: (1) biophysical underpinnings, (2) implications for information processing, and (3) experimental control of precise spiking *in vivo*.

The essential biophysical underpinnings of single-neuron activity can often be understood through simplified biophysical descriptions (Herz et al., 2006). Furthermore, such reduced quantitative descriptions are needed in theoretical studies of neural networks. Integrate-and-fire (IF) models can capture the subthreshold membrane voltage and the spike times of some pyramidal cells under current injection *in vitro* fairly well (Badel et al., 2008). This model class can be used in analytical approaches (Brunel and Hakim, 2008), but also in large-scale network simulations (Gerstner et al., 2014). In such studies, using a faithful model for the single-neuron activity is likely as important as using the correct connectivity. Computational studies have shown, for instance, that the details of the single-neuron model determine whether an autonomous network can maintain low firing rates (Latham et al., 2000; Kriener et al., 2014) and if the neurons synchronize or remain in an asynchronous state (Abbott and van Vreeswijk, 1993; Brunel, 2000).

The implications of precise spiking for information processing are an unfathomed issue in neuroscience. The debate over the significance of discharge timing gained traction some 20 years ago (Softky and Koch, 1993; Shadlen and Newsome, 1994; König et al., 1995), but the contradicting views on the topic could not be resolved back then. Some of the most direct evidence for a role of precise discharge timing comes from single-neuron-stimulation experiments. In such “reverse physiology” paradigms, it was found that a few extra spikes may affect the behavior of an animal: single-neuron nanostimulation in rat motor cortex (Brecht et al., 2004) and somatosensory cortex (Houweling and Brecht, 2008, Doron et al., 2014) could be detected by the animal; a similar hypersensitivity with respect to single-neuron activity was observed in the visual cortex of mice (Li et al., 2009).

To address the debate about the role of precise discharge timing, we need improved experimental control of precise spiking *in vivo*. Juxtacellular stimulation of single neurons is an effective stimulation technique, but so far it has limitations with respect to control of spike timing. Specifically, it was shown that, by applying combinations of rectangular current pulses, only the control of spike number, frequency, and regularity (Houweling et al., 2010; Doron et al., 2014; Doron and Brecht, 2015), but not exact timing, could be achieved.

Addressing the issues raised above, we first investigated whether juxtacellular stimulation with broadband (white-noise-like) currents can gain control over the precise spike times of the evoked train of action potentials (APs) in both somatosensory and motor cortices. The reliability in our *in vivo* experiments was lower compared with classical *in vitro* studies (Mainen and Sejnowski, 1995), probably due to the synaptic background noise, but is similar for *in vivo* juxtacellular and intracellular stimulation.

We then analyzed the spike train statistics and signal-processing capabilities of single cells using broadband and pure periodic cosine stimuli. Using estimates of the transfer function of an individual cell, we were also able to prescribe a spike train (with constraints on its statistics, but not on single spike times) that this neuron will fire.

Finally, we used a multicompartment model (Eyal et al., 2014) and a two-compartment model (Ostojic et al., 2015) to reproduce quantitatively the statistics measured *in vivo* and compared the performance of the latter model with a one-compartment model.

## Materials and Methods

##### Electrophysiology.

Most experimental procedures were performed as described previously (Houweling et al., 2010). Briefly, we used standard surgical and electrophysiological techniques to prepare animals [Wistar rats, *n* = 46, postnatal day 20 (P20)–P50 at the day of surgery, of either sex] for acute experiments. Animals were anesthetized with urethane (1.5–2.0 g/kg, i.p.). Recordings were performed in the barrel cortex (P 2.5 mm, L 5.5 mm relative to bregma) and motor cortex (A 1.5 mm, L 1.5 mm relative to bregma). Glass pipettes for juxtacellular and whole-cell recordings were filled with intracellular solution containing the following (in mm): K-gluconate 135, HEPES 10, Na2-phosphocreatine 10, KCl 4, MgATP 4, and Na3GTP 0.3, pH 7.2. The juxtacellular signal was amplified and low-pass filtered at 3 kHz with a patch-clamp amplifier (Dagan) and sampled at 25 kHz with a Power1401 data acquisition interface under the control of Spike2 software (CED). For intracellular recordings, series resistance and capacitance were compensated. During single-cell juxtacellular stimulation trials, a 1 s square-wave current pulse either added with Gaussian noise or with a harmonic wave on top was injected into a neuron through a glass pipette and current strength was adjusted (range 1–14 nA for juxtacellular stimulation and 100–1000 pA for whole-cell stimulation) to elicit a maximal number of APs without damaging the neuron. A sketch of the setup is given in Figure 1*A*.

All experimental procedures were performed according to German guidelines on animal welfare under the supervision of local ethics committees.

##### Stimulation with bandpass and cosine stimuli.

In most experiments, we used a broadband signal as follows:
with mean value *I*_{0}α and SD *I*_{0}σ; η(*t*) is a band-pass-limited white Gaussian noise with unit variance and cutoff frequency *f*_{c} = 100 or 200 Hz.

In some experiments we used a pure cosine stimulus as follows: with frequencies in the range of 100–1000 Hz.

##### Signal analysis and statistics of interest.

Most of our statistics are based on Fourier transforms of time series defined by the following:
where *z*(*t*) is either the stimulus *I*_{s}(*t*) or the spike train *x*(*t*). The latter is given in terms of the spike times *t _{i}* as a sum of δ functions as follows:
We define the instantaneous firing rate as the average over

*N*

_{trial}spike trains

*x*

_{i}(

*t*) belonging to the same stimulus and smooth it with a normalized Gaussian of 2.5 ms width (asterisk denotes convolution) as follows: The power spectrum is defined as follows: Here, the asterisk denotes complex conjugation and <> averaging over trials and stochastic signals. To quantify the similarity between spike trains in Fourier space, we define the spike-train—spike-train cross-spectrum as follows: The average runs over all stimuli and over combinations of distinct spike trains that correspond to the same stimulus. Analogously, the input–output cross-spectrum is given by the following: According to the Wiener–Khinchin theorem, the inverse Fourier transform ℱ

^{−1}yields the corresponding cross-correlation functions as follows: A measure that quantifies the response to a pure cosine stimulus is the vector strength (Goldberg and Brown, 1969): Here

*f*

_{s}is the frequency of the input,

*t*

_{j}are the spike times, and <> denotes averaging over trials. The absolute value of the vector strength yields |

**(**

*r**f*)| = 0 if the spike times are unrelated to the periodic stimulus, whereas a perfect locking to one specific signal phase leads to |

**(**

*r**f*)| = 1. Furthermore,

*r*(

*f*) can be regarded as a normalized Fourier coefficient of the spike train at the driving frequency and is related to the cross-spectrum under white-noise stimulation through

*S*= <

_{sx}*Ĩ**

_{s}*x̃*>/

*T*∼

**(**

*r**f*). To infer the frequency-dependent information transfer, we use the spectral coherence function: which is bounded between zero (no linear correlation between input and output) and one (input and output are linearly related in a noiseless fashion). The coherence can be used to give a lower bound on the mutual information rate for Gaussian stimuli (DeWeese and Bialek, 1995; Gabbiani, 1996) as follows: The integrand in Equation 13 is a monotonic function of the spectral coherence and can be interpreted as an approximate mutual information rate density (Bernardi and Lindner, 2015). Therefore, the shape of the coherence mirrors the neuron's information transfer properties in a frequency-resolved manner. To smooth the data,

*S*(

_{sx}*f*),

*S*(

_{xx}*f*),

*S*(

_{xixj}*f*) and the coherence are convolved with a Gaussian of 3 Hz width in all figures, where spectral measures are shown.

##### Spike train similarity (intrinsic variability).

To quantify the similarity of two spike trains, *a* and *b*, we use the coincidence measure (Kistler et al., 1997) which is defined as follows:
Here *N*_{a}(*N*_{b}) is the number of spikes occurring in a spike train *a*(*b*) of length *T*; *N*_{coin} is the number of spikes that coincide within a temporal precision of ±Δ; and *N*_{chance} = 2Δ*N _{a}N_{b}/T* is the expected number of spikes that would coincide by chance if both spike trains would have homogeneous randomly distributed spike times. The normalization factor reads as 𝒩 = (1 − 2ΔN

_{a}/T) and bounds Γ

_{a,b}by 1 (equal to 1 for identical spike trains). Due to the subtraction of

*N*

_{chance,}it is expected to yield Γ

_{a,b}= 0 for independent spike trains. We use Δ = 2.5 ms as precision for coincident spikes and define Γ = 〈Γ

_{a,b}〉 as the average over all stimuli and combinations of distinct spike trains corresponding to the same stimulus realization. When comparing experiment and simulation in terms of the coincidence factor, the best match is achieved when two criteria are fulfilled. First, simulation and experiment should have the same intrinsic variability. This is the case when the coincidence factor calculated for pairs of simulated spike trains Γ

_{ss}is as large as it is for pairs of experimentally measured spike trains Γ

_{ee}. Second, the model should reproduce the experimental spike times with a desired accuracy. This is achieved when the coincidence factor calculated for mixed pairs of experiment and simulation Γ

_{se}is as large as possible. In particular, Γ

_{se}should be as large as the intrinsic variability Γ

_{ee}. In other words, the simulation has to be as similar to the experiment as the experiment is to itself.

##### Spike extraction.

To extract the spike times, we filter the measured voltage by removing slow transients and frequency components shared with the input current. We typically remove the frequencies up to 60 Hz for the cases with constant input and frequencies up to 300 or 400 Hz for cases with fluctuating input. In the filtered trace, spike times are identified using a simple threshold criterion. Each time the voltage crosses the threshold from below is stored as spike time. To avoid artifacts arising from transients, we exclude the first and last 20 ms from our analysis. We defined a signal-to-noise ratio by SNR = spike height/SD(*V*), that is, the ratio of the height of putative spikes and the SD of the voltage, to quantify the quality of a recording. To be conservative, we accept data only if the average spike height exceeds the fourfold SD of the filtered voltage, corresponding to a signal-to-noise ratio >4. Because the filtering procedure excludes the power below the stimulus cutoff frequency, the latter could not be taken arbitrarily high and was limited to maximum 200 Hz for the broadband stimulation. For the experiments with pure cosine stimulation, this problem could be circumvented by excluding only the power in a narrow frequency band around the periodic stimulus and, therefore, we could go up to stimulation frequencies of *f*_{s} = 10^{3}Hz. Note that this spike extraction technique relies on the limited frequency range of the stimulus power. Therefore, one of the simplest stimuli, a sequence of narrow rectangular current pulses, would not allow for reliable spike extraction in this setting.

##### Multicompartment model.

For simulations, we use a multicompartment model that has been proposed by Eyal et al. (2014). This is a ball (soma) and two sticks (dendrite and axon) model implemented in the Neuron software package. It consists of 107 compartments [soma (1), axon (21 + 25) and dendrite (60)]. In particular, we refer to their case of an intermediate conductance load (ρ = (*G*_{dend} + *G*_{soma})/*G*_{axon} = 190), where *G*_{dend}, *G*_{soma}, and *G*_{axon} are the inverse input resistances of the isolated compartments, respectively. This dendrite has a length of 3000 μm and a diameter of 5 μm (parameters as in Eyal et al., 2014; Fig. 2, gray curves). For comparison with the experimentally studied cell shown in Figure 8, we inject a rescaled input current and additionally white Gaussian noise in the model to account for the intrinsic variability. The model parameters itself have not been modified.

##### Two-compartment model.

We also use a two-compartment model that contains one somatic compartment and one dendritic compartment (Ostojic et al., 2015). To allow for electrotonic coupling, both compartments are connected via an Ohmic conductance. *g*_{c}. The somatic compartment, characterized by the membrane voltage, *V*_{s}(*t*), is equipped with an active spike-generating mechanism, whereas the dendritic compartment, described by the voltage *V*_{d}(*t*), is purely passive (see the equivalent circuit in Fig. 1*B*):
Here the parameter α > 1 accounts for the fact that, in the juxtacellular situation, only a fraction of the injected current enters the cell. We assume that the soma and the dendrite have the same resting potentials and shifted the voltages to set it to zero. The extracellular potential is assumed to be constant at reference potential. Dividing Equation 15 by *g*_{s}Δ_{T} and Equation 16 by *g*_{d}Δ_{T} and measuring the voltages in units of the spike slope factor, *V*/Δ_{T} → *V*, yields the following dimensionless equation:
Here τ* _{s}* =

*C*and τ

_{s}/g_{s}*=*

_{d}*C*are the somatic and dendritic membrane time constants of the isolated compartments, respectively, and

_{d}/g_{d}*V*

_{Th}/Δ

*determines the threshold for*

_{T}*V*

_{s}, at which point the exponential starts to rise rapidly. The somatic spiking mechanism obeys a fire-and-reset rule: Whenever

*V*

_{s}>

*X*, we register a spike time and reset

*V*

_{s}→ 0. The parameter

*g*(

_{c}/g_{s}*g*) is the ratio of the coupling conductance and the somatic (dendritic) membrane conductance. The somatic input current

_{c}/g_{d}*I*

_{in}=

*I*

_{s}+

*I*

_{noise}consists of two parts. First, the external stimulus injected via the juxtacellular electrode that is given by

*I*

_{s}(

*t*). Second, to account for the random fluctuations arising in the soma

*in vivo*, we add intrinsic white Gaussian noise,

*I*

_{noise}(

*t*), with zero mean value 〈

*I*

_{noise}(

*t*)〉 = 0 and correlation function 〈

*I*

_{noise}(

*t*)

*I*

_{noise}(

*t′*)〉 = 2

*D*δ(

_{s}*t*−

*t′*). The dendritic current,

*I*

_{in,d}=

*I*

_{base}+

*I*

_{noise,d}, which accounts for the synaptic input that the dendrite receives, consists of a constant base current,

*I*

_{base}, and Gaussian white noise with 〈

*I*

_{noise,d}(

*t*)

*I*

_{noise,d}(

*t′*)〉 = 2

*D*δ(

_{d}*t*−

*t′*). In total, the two-compartment model involves the following nine intrinsic parameters: For a given set of these parameters, the biophysical parameter values can be calculated, for example, when assigning a fixed value to the spike slope factor Δ

_{T}and the parameter α that gives the fraction of current that enters the cell. We chose Δ

_{T}= 0.5 mV and α such that we obtain an input resistance in the model of 45 MΩ, as we have measured on average in intracellular experiments.

##### One-compartment model.

We compare the performance of the two-compartment model with that of the one-compartment model (Fourcaud-Trocmé et al., 2003). The one-compartment model results from Equation 17 in the limit of vanishing coupling conductance *g*_{c}. Because the one-compartment model still has to account for the synaptic input of the surrounding neurons, the constant input *I*_{base} is shifted to the somatic compartment. This yields the following equation:
The spiking mechanism remains the same as in the case of the two-compartment model.

##### Model fitting.

For comparison between experiment and simulation, we define a cost function as follows:
The individual terms have the form ‖Δφ* _{i}*‖/‖φ

*‖ where φ*

_{i}^{exp}_{i}is some measurable quantity. The term ‖Δφ

*‖/‖φ*

_{i}*‖ is the difference between simulation and experiment (Δφ*

_{i}^{exp}_{i}) normalized by the value from the experiment (φ

*), which makes different functions φ*

_{i}^{exp}_{i}comparable. The first term describes the relative deviations in power spectrum (φ

_{1}=

*S*), followed by the spike-train—spike-train cross-spectrum (φ

_{xx}_{2}=

*S*) and by the stimulus—spike-train cross-spectrum (φ

_{xixj}_{3}=

*S*). The last term in Equation 20 contains the relative deviation of the mean firing rates

_{sx}*r*. Because φ

_{1}, φ

_{2}, and φ

_{3}are frequency-dependent quantities, we define a norm as the area under the absolute value of the function up to a frequency

*f*

_{F}(we chose

*f*

_{F}= 1.5

*f*) as follows: To minimize Equation 20, we use the CMA-ES algorithm (covariance matrix adaptation evolution strategy; Hansen and Ostermeier, 2001; Hansen et al., 2003). The essential CMA-ES iteration steps are as follows. It draws an ensemble of points (feasible solutions) in the parameter space from a multivariate normal distribution centered at the actual favorite solution. In each iteration step, the candidate solutions are ordered according to their cost function values and weighted to prefer the best ones (evolution). These weighted candidates are used to calculate the new mean and the covariance from which a new multivariate normal distribution is calculated to draw feasible solutions in the next step. We use the CMA-ES implementation from the Python “cma” package (Hansen, 2015).

_{c}##### Goodness-of-fit.

Because we deal with finite datasets, the spectral measure φ_{i} will be subject to measurement noise. To quantify deviations between simulation and experimental data that go beyond this measurement noise, we define two functions φ_{i}^{±} = 〈φ* _{i}^{sim}*〉 ± 3SD(φ

*) based on the statistics of the measure in the simulated model. Typically, the simulated measure will be within the region that is limited by φ*

_{i}^{sim}_{i}

^{±}. The deviation of the experimental data, φ

*, from this region provides us with a measure for the goodness-of-fit: where Θ denotes the Heaviside function. Intuitively, Λ can be regarded as a percental deviation between simulation and experiment.*

_{i}^{exp}##### Generating noise that evokes a prescribed spike train.

Using a linear ansatz, input and output of a system can be related according to the following:
which, after applying a Fourier transformation, turns into:
In Equations 23 and 24, *K*(*t*) is the system specific linear response kernel and χ(*f*) its Fourier transform, known as susceptibility or the rate-modulation factor. On the one hand, this equation determines the susceptibility according to the following:
On the other hand, given a certain desired spike train, *x*_{†}(*t*), and its Fourier transform *x̃*_{†}(*f*), we can estimate the most likely stimulus evoking *x̃*_{†}(*f*) to be the following:
To avoid stimulation with too large pulses that might damage or kill the cell, we constrain the stimulus to obey Gaussian statistics. To this end, we apply two somewhat ad hoc steps. First, instead of using the spike train *x*_{†}(*t*), we use a filtered version of it, *x̄*_{†}(*t*) = *N*(0,2.5 *ms*) * *x*_{†}(*t*) in Equation 26, which yields a continuous function in time. Second, for the resulting signal, *Ī*_{†}(*t*), we estimate the cumulative probability distribution, *P*(*x*) = Prob(*x* < *Ī*_{†}) and then pass it through a static nonlinearity as follows:
where *F* has the form:
Here
is the inverse of the cumulative probability distribution for a Gaussian variable with mean μ and variance σ^{2}. In this way, we have achieved that the resulting stimulus *I*_{D}(*t*) obeys a particular Gaussian probability density. However, due to the nonlinear transformation, *I*_{D}(*t*) no longer has the same cutoff frequency as *Ī*_{†}. Because a finite cutoff frequency is essential for the spike extraction procedure (see above), we have to further process the stimulus. To this end, we remove all power at frequencies beyond the desired cutoff frequency. This in turn changes again the profile of the probability distribution to a non-Gaussian shape so that the nonlinear transformation has to be applied again. When performing nonlinear transformation to a Gaussian and the constraint on high-frequency power in an iterative way, the stimulus converges to the desired statistics within 10 iterations (for a similar scheme, see Nichols et al. (2010)). Of course, compared with the stimulus resulting from the optimal linear reconstruction, we now deal with a modified variant of this stimulus and it has to be tested whether this stimulus is still effective in evoking the desired spike train *x*_{†}(*t*).

From Equation 26, it becomes clear that one has to know the susceptibility, χ, which is specific for the neuron under investigation, in advance before calculating the desired stimulus. In addition, we are not completely free in choosing the prescribed spike train and restrict ourselves to spike trains that have the same mean firing rate and coefficient of variation (*C*_{V}, SD over mean) of the interspike interval (ISI), which in turn requires us to measure these statistics for the specific neuron.

Because we have to know the characteristics of the neuron before calculating *I*_{D}(*t*), we divide the experiment in two phases. In phase 1, we stimulate the neuron with 10 different white Gaussian noise stimuli (Eq. 1). Each of these stimuli is applied 10 times so that we gain a dataset that consists of 100 trials from which, after spike extraction, the mean firing rate, the *C*_{V}, the cross-spectrum, and, thus, via Equation 25, the susceptibility can be calculated. Using the statistics from phase 1, we generate the desired spike train as a renewal process; practically, we simulate a white-noise-driven perfect IF model that obeys these statistics (Vilela and Lindner, 2009a).

With the desired spike train and the susceptibility of the neuron, we can calculate the stimulus (same mean, μ_{exp}, and variance, σ_{exp}, as applied in phase 1) that should evoke the desired spike train in this neuron. In phase 2, this stimulus is applied to the neuron several times to generate raster plots and to calculate the reliability in evoking the desired spike times.

## Results

### Reliable spike extraction in juxtacellular broadband stimulation

When spikes are evoked by long DC current steps in conventional juxtacellular stimulation, determining spike times is relatively straightforward because stimulation artifacts are limited to the onset and offset of current steps. Extraction of spikes becomes a nontrivial matter, however, when a broadband nanostimulus is applied juxtacellularly. Because our stimulus frequencies extended only up to a certain cutoff, we could filter out this frequency band. We first ensured that this filter procedure for spike extraction is reliable.

Figure 2 illustrates the filtering procedure for spike extraction for juxtacellular stimulation with strong input fluctuations (Fig. 2*A1–A3*) and for an intracellular recording with smaller input fluctuations (Fig. 2*B1–B3*). In the raw data, Figure 2*A1* APs are hidden in large fluctuations of the voltage trace. Figure 2*A2* displays the trace after applying a filter yielding a much cleaner sequence of spikes clearly distinguished from the baseline. Due to the variability in AP width and the filter procedure, the AP height extracted in this way fluctuated as is illustrated by the histogram on the right of Figure 2*A2*. To further test the filtering procedure, we applied it also to intracellularly measured voltage traces (raw trace in Fig. 2*B1*), resulting in a spike sequence of comparable variability in spike height (histogram at the right side of Fig. 2*B2*). In particular, the *C*_{V} (SD over mean) of the two height distributions was comparable (*C*_{V} = 0.27 for juxtacellular, *C*_{V} = 0.2 for intracellular recordings).

The quality of spike extraction and the effect of the filtering procedure are shown for the whole juxtacellular dataset in Figure 2*A3*. Here, we show the SNR before and after filtering. The height of the gray boxes in the two top rows denotes the SD of the particular trace (including the large excursions during the APs). For the raw traces, we used the spike height measured at constant stimulation (σ = 0) to calculate the SNR at larger σ. As a conservative criterion for spike extraction, we used SNR = 4 as the threshold (dashed line); that is, only spike trains with sufficiently high amplitudes were used in our analysis. Because the filtering procedure excludes power below the stimulus cutoff frequency, the latter could not be taken as arbitrarily high and was limited to the maximum 200 Hz for the broadband stimulation. For the experiments with pure cosine stimulation, this problem could be circumvented by excluding only the power in the narrow frequency band around the periodic stimulus and, thus, we could go up to stimulation frequencies of *f* = 1 kHz.

### Precise control of spike timing by juxtacellular stimulation

Next, we investigated the extent to which we have control of the neuron's behavior by stimulating it with a broadband stimulus (Eq. 1) and changing the parameter (σ, α). For each neuron, *I*_{0} was adjusted to yield a firing rate approximately in the range of 10–30 Hz and was kept constant while recording this neuron. In Figure 3*A*, a reconstructed juxtacellularly stimulated layer 2/3 pyramidal neuron in the motor cortex is shown. For a constant current, the raster plots in Figure 3*B* reveal randomly distributed spikes; their timing was determined by the synaptic input and intrinsic noise, which were different in each trial. In marked contrast, a fluctuating input leads to more reliable spike timing (Fig. 3*C*). The high spike-time reliability under noisy juxtacellular stimulation also becomes apparent in the correlation between stimulus and spike train *C*_{sx}(τ) (Fig. 3*D*) and in the correlation among different spike trains under the same (frozen) stimulus *C _{xixj}*(τ) (Fig. 3

*E*). In summary, as shown previously in the classical

*in vitro*study by Mainen and Sejnowski (1995) for intracellular stimulation, our data demonstrate that a high spike-time reliability can be achieved using noisy juxtacellular stimulation

*in vivo*.

To examine the roles of mean and SD of the noisy stimulation, we used various combinations of the stimulus parameters α and σ (Fig. 4). The spike number increased with the mean of the applied current. Therefore, positive (negative) currents depolarize (hyperpolarize) the cell (Houweling et al., 2010). More importantly, the timing reliability across trials grew gradually with increasing SD. It was largest for a large mean value (ensuring a high mean firing rate) and a high amplitude of the time-dependent part of the noise. Even in the case of the highest amplitude noise with the largest mean, however, locking to the stimulus was never perfect because external synaptic input and intrinsic noise were still present. Pooled population data for the mean firing rate versus the mean input parameter α are shown in Figure 5*A*; Figure 5*B* displays the firing rate versus the coincidence measure Γ. To better elucidate the relation between firing rate and coincidence measure, we present separately data for small stimulus fluctuations (Fig. 5*A1*,*B1*) and large stimulus fluctuations (Fig. 5*A2*,*B2*). In all cases, the firing rate increased monotonically with α, but this relationship became more linear for stronger stimulus fluctuations (i.e., a larger value of σ) due to the well known linearization effect of noise (Longtin, 2000; Vilela and Lindner, 2009a). From Figure 5*B1*, it becomes evident that stimulation with a small or intermediate amplitude σ resulted in little coincident evoked spiking despite a broad range of firing rates. When the mean firing rate was high enough, we obtained strong synchronization among trials only for large σ (used in Fig. 5*B2*). In both Fig. 5*B1* and Fig. 5*B2*, we also show data for the equivalent intracellular experiment that display similar dependence between firing rate and coincidence measure. Interestingly, both methods resulted in similar maximal Γ values. This means that, in our *in vivo* setting, the feasible reliability by the juxtacellular stimulation technique is similar to that of the classical intracellular stimulation.

### Neurons transmit high-frequency stimulus information

How well does the neuron transfer certain frequency components of the input? To address this question, we calculated the cross-spectra *S*_{sx}(*f*) between input *I*_{s}(*t*) and spike train *x*(*t*) for band-pass stimuli with σ = 1 and α = 1 and the vector strength for pure cosine stimulation (Fig. 6).

Focusing on experiments that yielded a minimal reliability of Γ > 0.15, we found that individual cross-spectra were either flat or even increased with frequency up to the cutoff frequency of *f*_{c} = 100 Hz and so does the average over all cross-spectra (colored, dashed line in Fig. 6*A*). This is surprising because the cells' firing rates barely exceeded 40 Hz.

To further test the high-frequency encoding capabilities of the cells, we used cosine stimuli that permit to explore much higher frequencies (100 Hz ≤ *f*_{s} ≤ 1000 Hz) and plot the vector strength normalized to its value at *f*_{s} = 100 Hz. Because there was a large variability in the frequency dependence, we sorted the data into curves that show a pronounced (initial) increase (Fig. 6*B*) or not (Fig. 6C). The different curves of the vector strength in Figure 6*B* showed shallow maxima around *f*_{s} = 600 Hz, whereas some of the curves in Figure 6*C* started declining between 200 and 400 Hz. This decline, however, was only slight and it was difficult to identify a clear cutoff frequency from the experimental data. Given the strong cell-to-cell variability, there was no principal difference with respect to high-frequency encoding for the data measured in motor cortex (blue) and somatosensory cortex (red).

The cross-spectra inspected above tell us how much signal power is transferred in each frequency bin, but do not quantify the flow of information. A frequency-resolved measure for information transmission is the coherence function (Borst and Theunissen, 1999) given in Equation 12 as the ratio of the squared cross-spectrum and the product of power spectra of stimulus and spike train, respectively. For most cells, the coherence function did not show any decrease, but rather attained a flat shape (Fig. 7*A*). This implies an equally strong coding of information for slow (low-frequent) and fast (high-frequent) components of the input stimulus and stands in marked contrast to the low-pass coding observed in simple model neurons (Vilela and Lindner, 2009b).

The coherence function can also be used to calculate a lower bound on the mutual information rate (MIR; Eq. 13). Remarkably, when plotted versus the coincidence measure Γ (Fig. 7*B*), the MIR approximately followed a linear relationship, MIR ∼ Γ. The values of the lower bound that were achieved for strong coincidence go beyond 100 bits/s and are thus comparable to information rates measured in various sensory neurons (cf. Table 2 in Borst and Theunissen (1999)). Note that for strong stimulation, resulting in only little trial-to-trial variability, the true mutual information rate can be significantly higher than the lower-bound estimate because of nonlinear encoding that is not described by the coherence function (Borst and Theunissen, 1999; Bernardi and Lindner, 2015).

### Two-compartment model can reproduce neuronal behavior

Theoretical studies of the neural response to time-dependent stimuli discussed at length how features such as temporal correlations in the background noise (Brunel et al., 2001; Fourcaud-Trocmé et al., 2003), noise coding (Lindner and Schimansky-Geier, 2001; Boucsein et al., 2009), or cooperativity among sodium channels (Ilin et al., 2013) may endow the neuron with a slow decay or even a saturation of its rate-modulation factor (susceptibility) at high frequencies. However, none of the suggested one-compartment models would reproduce the feature seen in our experiments, namely, that the modulation factor increases with frequency well >100 Hz. Indeed, we were unable to reproduce the spike-train statistics under nanostimulation with a simple one-compartment IF model that included colored background noise (as in the study by Brunel et al., 2001), high onset speed of the AP (by an appropriate choice of the slope parameter in the exponential IF model as in Fourcaud-Trocmé et al., 2003), or both.

A recent computational study by Eyal et al. (2014) and *in vitro* experiments and modeling on Purkinje cells by Ostojic et al. (2015) revealed that the presence of a purely passive dendritic tree can enhance the response at high frequencies. We tested whether similar models as used in these studies can reproduce, not only the rate modulation factor, but also all of the second-order statistics that were accessible in our *in vivo* experiments: spike train power spectrum, cross-spectra between stimulus and spike train and between two different spike trains (i.e., two different trials) for a frozen stimulus, and the coherence function.

We used the rescaled input current and additional white Gaussian noise to stimulate the multicompartment model by Eyal et al. (2014) without any further parameter tuning. Doing so, we observed for a specific cell from the motor cortex in Figure 8 a surprising agreement between simulation results of the multicompartment model (blue) and data (black). We also used a two-compartment exponential IF model (equivalent circuit in Fig. 1*B*) and a one-compartment exponential IF model, for which the parameters were estimated by an optimization procedure (see Materials and Methods). In particular, the two-compartment IF model (red curves) could reproduce quantitatively: (1) the increase of the cross-spectrum with frequency (Fig. 8*D*); (2) the flat shape of the coherence as a function of frequency (Fig. 8*E*); and (3) the nonlinear effect of the stimulus on the power spectrum (Fig. 8*B*) and the cross-spectrum between trials (Fig. 8*C*) even outside of the stimulus frequency band (i.e., here for *f* > 100 Hz). In contrast, the one-compartment model (green curves) could also reproduce the power spectrum (Fig. 8*B*) and the cross-spectrum between trials (Fig. 8*C*), but failed to reproduce the input–output statistics (Fig. 8*D*,*E*). The cross-spectrum did not reproduce the increase with frequency (Fig. 8*D*) and the coherence displayed a pronounced low-pass behavior in marked contrast to the experimental data.

A less abstract illustration of how well the two-compartment model captures the time-dependent response of the neuron can be gained by looking at the instantaneous firing rate, *r*(*t*); that is, the time-dependent probability density for observing a spike at a certain instant in time (Fig. 8*A*). The two-compartment model and the experiment show both rapid, signal-related changes and their firing rates agree well. A similar agreement in spike train statistics between the two-compartment model and experiment can also be achieved for other cells from motor cortex and somatosensory cortex and for a higher cutoff frequency of the stimulus (*f*_{c} = 200 Hz), whereas the one-compartment model systematically performs worse in reproducing all the spectral statistics.

To quantify the performance of both the one- and the two-compartment model, we plot Λ, a measure for the goodness-of-fit (Eq. 22), for 16 cells (Fig. 9*A*). This measures the deviation between experiment and the region one would expect for the simulations (red transparent region in Fig. 8). In Figure 9*A*, the cells were sorted according to the Λ values corresponding to the two-compartment model. Using this measure, the superiority of the two-compartment model became apparent. However, satisfying two-compartment fits could not be found for all cells; beyond cell index 10, we observed a distinct increase in Λ with cell index. Therefore, we divided the population into cells (blue symbols) below the corresponding threshold (dashed line in Fig. 9*A*) and cells (red symbols) above threshold. Note that the two-compartment model provided the better fit even for the cells for which statistics could not be reproduced well (red symbols).

Deviations between model and experimental data can be also described in terms of the spiking times and the associated coincidence measure Γ. Using the same color code as in Figure 9*A*, the performance of both models in terms of the Γ factors is given in Figure 9, *B* and *C*. Figure 9*B* shows the trial-to-trial reliability within simulations (Γ_{ss}) versus that of the experiment (Γ_{ee}). For both models, the data are aligned to the diagonal, which corresponds to the same degree of reliability in model and experiment. Clearly, the two-compartment model remains closer to the diagonal, indicating that it reproduces the trial-to-trial reliability better than the one-compartment model. Figure 9*C* shows the similarity between simulation and experiment, Γ_{se}, versus Γ_{ee}, revealing again superior performance of the two-compartment model, the values of which are close to the diagonal. When comparing the two-compartment model for different cells, those performing well in terms of the spectral statistics (blue symbols below threshold in Fig. 9*A*) tended to be distributed closer to the diagonal than those with suprathreshold values of Λ (red symbols above threshold in Fig. 9*A*). In Figure 9*D*, the different parameters that belong to the 10 cells below threshold are presented. For some of the parameters, specific values can be found in the literature (Rall, 1959; Koch et al., 1996; Badel et al., 2008; Bekkers, 2011; Harrison et al., 2015; Ostojic et al., 2015), whereas, for others (*D*_{s}, *D*_{d}, *I*_{base}), the order of magnitude can be estimated based on biophysical considerations. For all parameters, we found a reasonable agreement between the expected order of magnitude and the values shown in Figure 9*D*. However, because the two-compartment model is still a strongly simplified biophysical description, these values should be regarded as effective parameter values only.

### Neurons can be stimulated to generate prescribed spike times *in vivo*

As outlined in the Materials and Methods section, we could stimulate the neuron to fire a prescribed spike train with a high accuracy. In this two-phase experiment, we first estimated the mean firing rate, the *C*_{V} of the ISI, and the transfer function of an individual cell (phase 1). We then generated a target spike train that had about the same rate and *C*_{V}. Furthermore, by means of the transfer function and an iterative procedure, we computed a Gaussian stimulus with vanishing power beyond our cutoff frequency that approximately elicited this spike train. Finally, we applied the computed stimulus several times to the neuron under investigation (phase 2).

In total, we used five neurons from motor cortex from which we created seven datasets (two cells have been explored with two stimulus intensities). Figure 10*A* shows the stimulus and Figure 10*B* the corresponding raster plot for one example cell for which the spikes appeared with high precision close to the prescribed times. Population data for the coincidence factors measured in the two phases of the experiment (Fig. 10*C*) show that the trial-to-trial similarity in phase 2 was at least as large as in phase 1 and for most cells even substantially larger. This is due to the still somewhat pulsatile nature of the computed stimulus (Fig. 10*A*). The similarity of the phase 2 trials with the target spike train was even higher than the self-similarity in phase 2, which is plausible because the desired spike train lacks the intrinsic noise of the individual trial. In conclusion, this method demonstrates the ability to prescribe spikes in a precise and reliable manner using juxtacellular stimulation *in vivo* in cortical neurons.

## Discussion

In the present study, we investigated whether and under which conditions juxtacellular stimulation can evoke reliably spike timing in single sensory and motor pyramidal neurons *in vivo*. We found that nanostimulation in the form of strongly fluctuating inputs with Gaussian statistics can evoke reproducible spike timing comparable to what has been achieved previously by intracellular stimulation *in vitro* (Mainen and Sejnowski, 1995) and by visual stimuli *in vivo* (Buracas et al., 1998). This is especially important in experiments in which the relevance of single spikes for a behavioral response is tested (Brecht et al., 2004; Houweling and Brecht, 2008; Doron et al., 2014).

In addition, we investigated the ability of single neurons to track high-frequency modulations and to encode information on rapidly changing signals. Consistent with previous *in vitro* (Köndgen et al., 2008; Boucsein et al., 2009; Ilin et al., 2013) and *in vivo* (Tchumatchenko et al., 2011) results, we found that the neuronal response (i.e., modulation of the instantaneous firing rate) and information transmission (spectral coherence function) is not limited to frequency bands below the mean firing rate or inverse membrane time constant. The maximum in the vector strength observed in Figure 6*B* approximately agrees with recent *in vitro* observations for Purkinje cells (Ostojic et al., 2015).

Different explanations for the transmission of high frequencies have been suggested, among them the effect of temporal noise correlations (Brunel et al., 2001) and the sharpness of spike onset (Fourcaud-Trocmé et al., 2003; Naundorf et al., 2006; Boucsein et al., 2009). However, our attempts to capture the statistics of our data in a one-compartment model with intrinsic colored noise and/or with arbitrarily sharp onset dynamics failed. Only if a second compartment representing the dendrite were included in the model could we reproduce the instantaneous firing rate, power spectra, cross-spectra between trials, and input–output statistics with great accuracy. Our results are in qualitative agreement with the recent study by Ostojic et al. (2015), who investigated the high-frequency transmission of time-dependent stimuli in Purkinje cells.

Coming back to the issue of spike control, as demonstrated by stimulation with white Gaussian noise stimuli, we were able to generate reproducible spike trains when stimulating the neurons with fluctuating Gaussian stimuli. However, using this approach, the exact timing of the APs remains uncontrolled in advance because the applied stimuli are random (frozen) noise realizations. Here, we extended the feature of generating reproducible spike trains to spike trains of which the AP timing is prescribed in advance. Specifically, we calculated the stimulus to evoke a desired spike train by considering the susceptibility of the neuron under investigation and applying an iterative procedure to ensure Gaussianity and finite frequency support (cutoff frequency). The desired spike train had similar statistical characteristics (mean firing rate and ISI *C*_{V}) as the spike trains evoked by broadband nanostimulation. Repetitive application of this stimulus forced the neuron to generate spike trains that showed a larger intrinsic reliability (trial-to-trial similarity) than was the case for stimulation with frozen white Gaussian noise. In addition to the high intrinsic reliability, the generated spike trains were evoked precisely at the desired spike times. Remarkably, the similarity between desired and evoked spike train was therefore at least as good as the trial-to-trial reliability in the white noise experiments.

We limited our analysis to spike trains with a similar statistics (rate and *C*_{V}) and found that, although our procedure of the stimulus computation is somewhat approximate, such a target spike train can be evoked with high accuracy and temporal precision. It is an interesting task to investigate how far this can be extended to target spike trains that deviate in their statistics significantly from what the inspected neuron would do under a Gaussian white noise stimulus with a given mean and variance.

In summary, we showed here that juxtacellular stimulation *in vivo* allows for high control over the activity of individual stainable neurons. This might be of particular interest for reverse physiology experiments, in which the influence of single-neuron activity on neuronal network dynamics or the animals' behavioral response is investigated.

## Footnotes

This work was supported by Bundesministerium für Bildung (Forschung Grant 01GQ1001A), the Humboldt Universität zu Berlin, and the Deutsche Forschungsgemeinschaft (a Gottfried Wilhelm Leibniz Prize).

- Correspondence regarding theory and modeling should be addressed to Benjamin Lindner, Bernstein Center for Computational Neuroscience Berlin, Philippstr. 13, Haus 2, 10115 Berlin, Germany. benjamin.lindner{at}physik.hu-berlin.de
- Correspondence regarding experiments should be addressed to Guy Doron, NeuroCure Cluster of Excellence, Charité Cross Over - Campus Mitte, Charitéplatz 1 10117, Berlin, Germany. guy.doron{at}charite.de