Oscillations and variability in neuronal systems: interplay of autonomous transient dynamics and fast deterministic fluctuations

Neuronal systems are subject to rapid fluctuations both intrinsically and externally. These fluctuations can be disruptive or constructive. We investigate the dynamic mechanisms underlying the interactions between rapidly fluctuating signals and the intrinsic properties of the target cells to produce variable and/or coherent responses. We use linearized and non-linear conductance-based models and piecewise constant (PWC) inputs with short duration pieces. The amplitude distributions of the constant pieces consist of arbitrary permutations of a baseline PWC function. In each trial within a given protocol we use one of these permutations and each protocol consists of a subset of all possible permutations, which is the only source of uncertainty in the protocol. We show that sustained oscillatory behavior can be generated in response to various forms of PWC inputs independently of whether the stable equilibria of the corresponding unperturbed systems are foci or nodes. The oscillatory voltage responses are amplified by the model nonlinearities and attenuated for conductance-based PWC inputs as compared to current-based PWC inputs, consistent with previous theoretical and experimental work. In addition, the voltage responses to PWC inputs exhibited variability across trials, which is reminiscent of the variability generated by stochastic noise (e.g., Gaussian white noise). Our analysis demonstrates that both oscillations and variability are the result of the interaction between the PWC input and the target cell’s autonomous transient dynamics with little to no contribution from the dynamics in vicinities of the steady-state, and do not require input stochasticity.


Introduction
Variability in neuronal activity has been observed at all levels of neuronal organization Calvin and Stevens (1967); Shadlen and Newsome (1998); Ito et al. (2020); Churchland et al. (2006); Faisal et al. (2008); Romo et al. (2003); Goris (2014); ; Fox et al. (2006); ; Britten et al. (1992); Softky and Koch (1993); Arieli et al. (1996); Cohen and Kohn (2011). However, the mechanisms underlying the generation of variable voltage responses to fluctuating inputs are not well understood. In particular, it is not well understood which aspects of the variable voltage responses of neuronal systems are due to their intrinsic properties (e.g., ionic currents), which aspects are due to properties of the external inputs, and how the two interact.
The goal of this paper is to address these issues in the context of single cells receiving external current inputs. Through a series of case studies, we systematically investigate the role of the intrinsic properties of the target cells in the generation of oscillatory and variable responses to fluctuating inputs, and how the response variability is controlled by these intrinsic properties.
The cellular intrinsic and dynamic properties can be uncovered and characterized by the use of constant inputs. The voltage responses to constant input currents consist of a transient phase followed by a stationary phase Day et al. (2009). Examples of transient behavior are monotonic changes towards equilibrium, overshoots (or sags), and damped oscillations.

Action Editor: Pulin Gong
We collectively refer to them as the autonomous transient dynamics. These two phases reflect different ways in which the cellular intrinsic properties interact with the input and therefore different ways in which the voltage response encodes information about these cellular intrinsic properties. We argue that the autonomous transient dynamics plays a key role in determining the cell's variable response to external inputs, while the cell's stationary response to constant inputs plays at most a minor role. We additionally show that variability emerges in response to multiple presentations of (deterministic) piecewise constant (PWC) input functions having exactly the same constant pieces arranged in different, arbitrary orders across trials. This provides a way of elucidating the mechanistic role of both the cellular intrinsic properties through the corresponding autonomous transient dynamics in generating the response variability. This also supports the idea that variability, although it is present in response to stochastic inputs, is not inherently a stochastic phenomenon.
A prototypical example of the role of the transient dynamics of individual cells (autonomous transient dynamics) in determining their response to external inputs is given by subthreshold resonance Richardson et al. (2003); Rotstein (2014Rotstein ( , 2015. Subthreshold resonance refers to the ability of cells to exhibit a peak in the impedance amplitude profile (curve of the impedance amplitude as a function of the input frequency, Fig. 2-b3 and -c3) for a preferred (resonant) input frequency Hutcheon et al. (1996); Hutcheon and Yarom (2000) (see Appendix 1). Subthreshold resonance emerges in both cells having intrinsic oscillatory properties (e.g., two-dimensional linear cells with complex eigenvalues, Fig. 2-c3) and cells lacking intrinsic oscillatory dynamics Richardson et al. (2003); Rotstein and Nadim (2014) (e.g., two-dimensional linear cells having real eigenvalues, Fig. 2-b3), but exhibiting an overshoot (or sag). The effective time scale governing these autonomous transient dynamics is key for the generation of subthreshold resonance and the determination of the resonant frequency, which reflects balances between interacting processes. Richardson et al. (2003); Rotstein (2014Rotstein ( , 2015. Neurons and neuronal networks are subject to intrinsic random fluctuations Van Kampen (2011); Deco et al. (2009) ;Faisal et al. (2008) such as random opening/closing of ionic channels Shalinsky et al. (2002); White and Haas (2001); Dorval and White (2005); DeFelice (1981); Chow and White (1996); ; Schneidman et al. (1998); Sigworth (1980) and random synaptic inputs from other neurons in the network Calvin and Stevens (1967); Brunel et al. (2001); Destexhe et al. (2004); Fellous et al. (2003); Yarom and Hounsgaard (2011);Hong et al. (2012); Pena et al. (2018). Random fluctuations are also exerted by external factors. In mathematical models, these fluctuations are incorporated as noise interacting in various ways with the underlying deterministic dynamics Destexhe and Rudolph-Lilith (2012).
Common to these phenomena is the idea that noise repeatedly "kicks" trajectories and keeps them away from their stationary solutions (when they exist) to create alternative patterns (irregular versions of the underlying noiseless pattern or qualitatively different patterns). In practice, external noise interacts with the cell, which responds by integrating the intrinsic (deterministic) and input (stochastic) components (similarly to what was described above). A prototypical example of a stochastic input to a neuron is the Ornstein-Uhlenbeck (OU) process Uhlenbeck and Ornstein (1930); Gardiner (1985); Risken (1989) (Appendix 2). While the response's expected value coincides with the stationary solution of the noiseless system, the response's variance involves a combination of the (noiseless) model parameters and the variance of the Gaussian white noise input. In other words, the response variability is controlled by the autonomous transient dynamics. For higher-dimensional OU process, the determination of the dependence of the variability with the autonomous transient dynamics is more difficult than for one-dimensional OU processes given the complexity of the covariance formulas. This type of determination is analytically not possible for noise-driven nonlinear systems (including voltage-dependent nonlinearities or conductance-based synaptic inputs). These tasks require a different approach and the development of a conceptual framework that allows the investigation of the response properties in terms of the properties of the participating building blocks: the target cells and the fluctuating inputs. The autonomous transient dynamics are a key element of this framework.
Here, we use piecewise constant (PWC) input currents ( I ) with short duration pieces and variable amplitudes that allow for the autonomous transient dynamics to develop in response to each constant piece input. When the set of amplitudes ( ) are normally distributed and the durations are small enough, I is an approximation to Gaussian white noise Allen et al. (1998). However, in this paper we are not interested in this limit. We first show that an additive input I with normally distributed amplitudes evoke oscillatory voltage responses in cells whose stable equilibria are either foci or nodes (displaying damped oscillations and overshoots in response to constant inputs, respectively). We note that the former captures results described in Gardiner (1985) for Gaussian white noise (see also Thomas and Lindner (2019) for two-dimensional linear cells with damped oscillations). We then use the same input to investigate the response properties of a nonlinear (piecewise linear, PWL) model Rotstein (2014), which mimics the voltage-dependencies present in neurons, and a linear model receiving a (multiplicative) conductance-based synaptic-like current. We explain how oscillations are nonlinearly amplified or attenuated in these models. Next, we investigate how the voltage response variability is linked to the properties of the autonomous transient dynamics. All trials for a given protocol have the same set of linear piece amplitudes , but the order of the constant pieces is different for each trial. Each protocol consists of a subset of all possible order permutations of the elements of the set . This is the only source of uncertainty in the process. We compute the peak-and-troughs voltage response profiles P consisting of the set of peaks and troughs of the voltage response evoked by each constant piece, arranged in an appropriate way for comparison. We analyze the dependence of the variability of the P patterns across trials with the properties of the autonomous transient dynamics, and we use these results to understand how the response variability depends on the model parameters. More specifically, this variability results from the multiple different ways a target cell reacts to the same constant inputs (across trials) as it transitions from the response to the previous constant piece. Finally, we show how frequency-dependent inputs affect these processes.

Models
In this paper we use relatively simple biophysically plausible models describing the subthreshold dynamics of individual neurons subject to both additive and multiplicative inputs.

Linear model: additive input
For the individual neurons we use the following linearized biophysical (conductance-based) model Rotstein and Nadim (2014); Richardson et al. (2003) where v (mV) is the membrane potential relative to the voltage coordinate of the fixed-point (equilibrium potential) of the original model, w (mV) is the recovery (gating) variable relative to the gating variable coordinate of the fixed-point of the original model normalized by the derivative of the corresponding activation curve, C ( F/cm 2 ) is the specific membrane capacitance, g L (mS/cm 2 ) is the linearized leak conductance, g 1 (mS/cm 2 ) is the linearized ionic conductance, 1 (ms) is the linearized gating variable time constant and I(t) ( A/cm 2 ) is the time-dependent input current. In this paper we consider resonant gating variables ( g 1 > 0 ; providing a negative feedback effect). The linearization process for conductance-based models for single cells has been previously described in Richardson et al. (2003); Rotstein and Nadim (2014). We refer the reader to these references for details.

Piecewise linear (PWL) model: additive input
To account for nonlinear effects we extend the linear model (1)-(2) to include a piecewise linear function F PWL Rotstein (2014) where where v c is the cutting point (or breaking point) of the PWL model. This model has been used to investigate the dynamic mechanisms underlying the nonlinear amplification of the resonant voltage responses to sinusoidal inputs Rotstein (2014) and captures the nonlinear amplification effects of the resonant voltage responses off two-dimensional models with parabolic-like voltage nullclines Rotstein (2015). Note that the PWL model has linear dynamics for v < v c and therefore the PWL model becomes effectively linear if v c is large enough.

Conductance-based synaptic input model: multiplicative input
To account for the effects of conductance-based synaptic inputs we extend the model (3)-(5) to include a synaptic current where G syn (mS/cm 2 ) is the maximal synaptic conductance, E syn (mV) is the synaptic reversal potential and S(t) is the time-dependent synaptic input. We use E syn = 1 (excitatory synaptic current).

Piecewise constant (PWC) input functions with variable amplitudes
We use PWC functions with short-duration pieces as a tool to understand how the properties of the transient voltage response of individual cells to constant inputs (autonomous transient dynamics) control the response variability of these cells to time-dependent inputs with variable, abruptly changing amplitudes. White noise and related stochastic inputs (e.g. colored noise) satisfies these last two properties. PWC inputs with variable amplitudes and short durations provide a good balance between input variability and the ability of the autonomous transient dynamics to develop, and therefore they can be linked to the input amplitude that gave rise to them. In addition, PWC inputs with variable amplitudes allow us to use the same set of amplitudes in all trials for each protocol, arranged in different order for each trial, and therefore have a better control of the process. PWC functions with normally distributed amplitudes have been used to approximate white noise Allen et al. (1998) (see also Du and Zhang (2002)) and converge to white noise as the duration of all pieces approaches zero.

Piecewise constant input functions with normally distributed amplitudes
We partition the time interval [0, T max ] into N pieces of equal length Δ , and we define I (t) = k for t ∈ (t k−1 , t k ) for k = 1, … N . The values of the amplitudes { k } N k=1 of the constant pieces are normally distributed (zero mean and variance D) (Fig. 1, top, see also Fig. 3-a). For each protocol in our study we use multiple different input trial functions I (t) constrained to consist of different (random) permutations of the same set { k } N k=1 . As part of our analysis we use one specific rearranged version I ,step of I where the linear pieces { k } N k=1 are ordered in a monotonically increasing amplitude manner ( Fig. 1, middle). We use this as our notion of a reference ("ordered") input in the sense that I ,step has the minimum piece-to-piece variability (except possibly for the decreasing order).

Piecewise constant input functions with arbitrarily, but deterministically distributed amplitudes
The choice of normally distributed amplitudes ( Fig. 3-a) is motivated by the fact that in the limit of Δ → 0 , I approaches white Gaussian noise when the piece duration goes to zero. In order to decouple the effect of the autonomous transient dynamics evoked by abrupt transitions between constant pieces and the randomness of the signal, we use fully deterministic distributions of amplitudes within some range. More specifically, we use arbitrary order permutations of a number of constant pieces initially arranged in increasing order of amplitudes where the amplitudes are chosen according to deterministic rules. Each protocol (see Sect. 2.2.1) consists of a subset of all possible permutations (of the order of constant pieces), which is the only source of uncertainty in the protocol. We use two types of (deterministic) amplitude distributions for { k } N k=1 within the range [ min , max ] : equispaced ( Fig. 3-b) and bell-shaped-like ( Fig. 3-c). In the equispaced distribution the values of k cover the full range [ min , max ] and satisfy k+1 − k is equal for all k = 1, … , N − 1 . The "ordered" PWC function I ,step is linear (Fig. 3-b2). The bellshaped-like distribution was constructed from the random distribution (see the caption in Fig. 3). The corresponding "ordered" PWC function I ,step has an inflection point reflecting larger number of input amplitudes around zero ( Fig. 3-c3).

Impedance (amplitude) profile
The impedance (amplitude) profile is defined as the magnitude of the ratio of the output (voltage) and input (current) Fourier transforms : (-2,0,1,-1,2) as in row 2, but with a different (non-monotonic) distribution. c Passive cell. We used the following parameter values: C = 1 , g L = 0.25 . Red dots indicate the boundaries between constant pieces (there are no peaks and troughs). d 2D linear system exhibiting an overshoot in response to step-constant inputs. We used the following parameter values: C = 1 , g L = 0.25 , g 1 = 0.25 , 1 = 100 . e. 2D linear system exhibiting damped oscillations in response to step-constant inputs. We used the following parameter values: . In practice, we use the Fast Fourier Transform algorithm (FFT) to compute {x(t)} . Note that Z(f) is typically used as the complex impedance, a quantity that has amplitude and phase. For simplicity, here we used the notation Z(f) for the impedance amplitude. Since we deal with discontinuous inputs (PWC), we have checked if Gibbs phenomenon would be an issue. As a result, we did not identify any ringing effect in the Fourier transforms related to Gibbs phenomenon. The parameter values are the same as in Fig. 1. Row 1. Representative V traces (only 1000 ms are shown). The coral dots indicate peaks-and-troughs patterns computed as the maximum (minimum) of the V response if k > k−1 ( k < k−1 ). Row 2. Power spectra density (PSD) profiles for the sample V trace. Row 3. Impedance amplitude (Z) profiles for the sample V trace. a Passive cell ( f nat = f res = 0 ). We used the following parameter values: C = 1 , g L = 0.25 . b 2D linear system exhibiting an overshoot in response to step-constant inputs ( f nat = 0, f res ∼ 9Hz ). We used the following parameter values: C = 1 , g L = 0.25 , g 1 = 0.25 , 1 = 100 . c 2D linear system exhibiting damped oscillations in response to step-constant inputs ( f nat ∼ f res ∼ 8Hz ). We used the following parameter values: C = 1 , g L = 0.05 , g 1 = 0.3 , 1 = 100

Voltage and impedance (amplitude) envelope profiles
The upper and lower envelope profiles V +∕− ENV are curves joining the peaks and troughs of the steady state voltage response as a function of the input frequency f. The envelope impedance profile is defined as Rotstein (2014Rotstein ( , 2015 where A in is the input amplitude. For linear systems, Z ENV (f ) coincides with Z(f).

Voltage power spectral density
In the frequency-domain, we compute the power spectral density (PSD) of the voltage as the absolute value of its Fourier transform {v(t)} . We will refer to this measure as PSD or V PSD .

Peaks-and-troughs voltage response profiles
In each protocol, we compute the cells' voltage response to the multiple trials described above. For each trial, we The values of the amplitudes { k } N k=1 of the constant pieces are selected by following different rules (panels a to c). Each realization of I (t) corresponds to a permutation of the order of the set { k } N k=1 . I ,step (t) is the special case where k are organized in a non-decreasing order (used for reference). In all cases, Δ = 1 with a total time = 1,000,000 ms. a Random normal distribution with mean zero and variance one.  Fig. 1. The color dots indicate the peaks-and-troughs patterns for all trials reorganized so that the corresponding values of k from which they originate are ordered in a monotonically increasing manner. All dots for a given piece correspond to the same value of k . < P > (blue) is the mean value of the reordered peaks-and-troughs patterns for each linear piece. < P ,step > (green) is the peaks-andtroughs pattern corresponding to the (ordered) input function I ,step . b1 to d1. Variance (b1 to d1) of the peaks-and-troughs patterns in panel a. b1. Comparison between cells having a node and a focus. c1. Comparison between a cell having a node and the corresponding passive cell. d1. Comparison between a cell having a focus and the corresponding passive cell. b2 to d2. Normalized variance (b2 to d2) computed as the variance (b1 to d1) divided by the peak of the unforced cells' response to a step-constant input of amplitude 1. b2. Comparison between cells having a node and a focus. c2. Comparison between a cell having a node and the corresponding passive cell. d2. Comparison between a cell having a focus and the corresponding passive cell. We used the following parameter values: (i) C = 1 , g L = 0.25 , g 1 = 0.25 , 1 = 100 for the node, and (ii) C = 1 , g L = 0.05 , g 1 = 0.3 , 1 = 100 for the focus compute the sequence of the voltage response peaks and troughs consisting of the maximum value of the voltage v within the duration of a piece k if k > k−1 and the minimum of v within the range of a piece k if k < k−1 , respectively. We use these peaks-and-troughs patterns or profiles (e.g., Fig. 4 Fig 1 c-e) to characterize the voltage response patterns and their variability. We note that there are other possible metrics, including the whole v traces and the set of all peaks and troughs present in these traces (e.g., Fig. 4-c to -e). We found the peaks-and-troughs profiles as described above to be a simple and useful way to compare the voltage responses across trial (permutations).

Rearranged peaks-and-troughs voltage response profiles
In order to compare the voltage responses across trials we rearrange the peak-and-troughs voltage response profiles according to the increasing amplitude values of the elements of the set { k } N k=1 from which they were evoked. We use the notation P for the resulting rearranged voltage response patterns and P ,step for these voltage response profiles produced by I ,step (inputs with increasing order of constant piece amplitudes). In this way, we compare multiple different ways in which the response to a given constant piece is produced by transitions from pieces with different amplitudes (i.e., different initial conditions with respect to that constant piece). The differences in the voltage peak/troughs values across trials for each piece will be due to the differences in the values of V and the other variables at the arrival time of that piece (initial conditions) and therefore to the differences in the transient dynamics that they evoke.

Numerical simulations
We used the modified Euler method (Runge-Kutta, order 2) Burden and Faires (1980) with step size Δt = 0.01 ms. Preliminary simulations where lower step sizes than this have been tested did not change our results. All neural models and metrics, including phase-plane analysis, were implemented by self-developed Matlab routines (The Mathworks, Natick, MA) and are available in https:// github. com/ BioDa tanam ics-Lab/ imped ance_ input_ depen dent.  Fig. 1-e1 display damped oscillations before converging to the steady-state. We note that not all two-dimensional cells having a stable node or focus show prominent overshoots or damped oscillations. Fig. 1-c2 to -e2 show the voltage response of the same three cell types to PWC inputs with the same constant pieces as in Fig. 1-c1 to -e1, but the three intermediate pieces are arranged in a different order. The type of autonomous transient dynamics, which are activated at the transition points between input pieces, remains the same as in Fig. 1-c1 to -e1, but the response amplitude of the voltage response pattern increases for the N-and F-cells as compared to Fig. 1-d1 and -e1. This is due to a combination of factors that involve the differences in the jump sizes between linear pieces and the differences in initial conditions for each linear piece (e.g., Fig. 1-b). The values of the participating variables prior to the transition between linear pieces serve as initial conditions for the new regime. While for 1D system, the response amplitude is bounded by the steady-state for each linear piece regardless of the initial conditions, for 2D systems the amplitude response depends on the initial conditions not only of the main variable (v), but also the recovery variable that is not shown in the graph (e.g., Fig. 1-b).

Oscillatory voltage response properties to piecewise constant inputs with arbitrarily distributed amplitudes
The V patterns discussed above (Fig. 1, row 2) for a small number of input ( I ) partitions (large Δ = 200 ) display the fully developed autonomous transient dynamics of the three cell types to step-constant inputs. The extension of these results to the V patterns generated by inputs with a larger, more realistic number of partitions (smaller Δ ) is not straightforward since the smaller the partition size the less time the cell has to develop the autonomous transient behavior corresponding to each partition. Still, the transient dynamics reflect the model properties and the differences among models. Interestingly, not only the F-cell that shows damped voltage oscillations in response to step-constant inputs shows oscillatory activity in response to I , but also the N-cell that shows overshoots in response to step-constant inputs. This would be in principle not surprising since both cell types show resonance in response to sinusoidal inputs Richardson et al. (2003); Rotstein (2014), but the PSDs for the input patterns I with a smaller number of partitions ( Δ = 200 as in Fig. 1, row 3) shows no oscillations for the three responses (see insets in Fig. 2, row 2). The arbitrarily (randomly) distributed amplitudes allow the voltage response to explore a wide region of the phase-space and therefore capture information about the vector field governing the transient dynamics. The short duration of the pulses keeps the trajectory ("almost continuously") exploring the vector field by sequentially activating the autonomous transient dynamics, and this "permanent motion" prevents the interference of the regions of the vector field governing the dynamics in close vicinities of the steady-state. This mechanism is conceptually similar to the one governing the generation of resonance in response to oscillatory inputs in both linear and nonlinear systems Rotstein (2014Rotstein ( , 2015 where resonance can be observed in the absence of intrinsic oscillatory behavior Richardson et al. (2003); Rotstein and Nadim (2014); Rotstein (2014).

Uncovering the oscillatory dynamic properties of the target cells requires constant pieces with arbitrarily distributed amplitudes, but not amplitude randomness
Here we decouple the effects of the autonomous transient dynamics (in response to the sequence of input constant pieces) from the pieces' amplitude randomness on the ability of these PWC inputs to uncover the oscillatory (intrinsic) and resonant properties of the receiving cells.
To this end we use fully deterministic distributions of amplitudes within some range as described in Sect. 2.2.2. This leaves the choice of the subset of all possible permutations (number of trials) for each protocol as the only source of uncertainty in the input signal. Figure S1 ( Δ = 5 and Δ = 1 ) shows the response patterns to these PWC inputs (equispaced distributed, deterministic amplitudes) for the same parameter values as Fig. 2 (random amplitudes). The Z-profiles in these three figures are almost identical. The V PSD profiles in Figs. 2 and S1 are also very similar. The differences between the V PSD profiles in these figures and the ones for Δ = 5 are due to the larger value of the total time used there. Together, these results show that randomness is not needed to uncover the resonant and oscillatory properties of the receiving cells and these oscillatory properties emerge almost exclusively as the result of the sequential and fast activation of the cell's autonomous transient dynamics.

PWC inputs with arbitrarily distributed amplitudes capture the nonlinear amplification of the oscillatory voltage responses
The amplification of the voltage response oscillations to PWC inputs discussed so far was presented in the context of linear systems as the result of changes in the model parameter values (e.g, from N-to F-cells in Figs. 2-, S1-, -b and -c). Because the models used there are linear, the dependence of the voltage response amplitude on the input amplitude (D) is simple. Increasing values of D cause proportional increases in the voltage response amplitude so that the Z-profile remains unchanged.
Here we focus on the nonlinear amplifications produced in (nonlinear) models as the result of increasing values of the input amplitude D. To this end we use the piecewise linear (PWL) model (3)-(5), which is a continuous, nonlinear extension of the linear model (1)-(2) where the V-nullcline is "broken" (e.g., Fig. 5-b). It was shown in Rotstein (2014) that this type of models display nonlinear amplifications of the voltage response to sinusoidal inputs and capture similar phenomena observed in more realistic nonlinear models, in particular two-dimensional models having parabolic-like V-nullclines describing the subthreshold voltage dynamics Rotstein (2015).
We first introduce the ideas by examining the autonomous transient dynamics of the PWL model in response to a constant input I and then use these results to understand the nonlinear amplification of this PWL model in response to PWC inputs with increasing values of D. Figure 5-a and -b show the superimposed phase-plane diagrams for a PWL model (b) and the linear model (a) from which it originates for I = 0 (solid-red, baseline) and I = 1 (dashed-red). The w-nullcline is unaffected by changes in I. The trajectory (blue), initially at the fixed-point for I = 0 , evolves towards the fixed-point for I = 1 . In both cases (Fig. 5-a and -b) the voltage response exhibits an overshoot. The peaks (inset) occur when the trajectories cross the V-nullcline. For low enough values of I (lower than in Fig. 5-a and -b) the trajectory remains within the linear region (the trajectory does not reach the V-nullcline's "breaking point" value of V) and therefore the increasing values of I produce a linear voltage response amplification (no differences between the responses of the linear and PWL model; not shown). The nonlinear amplification is apparent for values of I for which the v-nullcline is high enough so the trajectory is able to cross from one linear regime (determined by the left piece of the dashed-red V-nullcline) where the trajectory is at t=0, to the other (determined by the right piece of the dashedred V-nullcline). The virtual fixed-point moved from the position where I = 0 and converged to the position where nullclines cross for I = 1 . Because the V-nullcline's "right piece" has a smaller slope than the "left piece" (the slope it would have if it were linear), the trajectory is able to reach larger values of V for the nonlinear system than for the corresponding linear one before turning around, and therefore the peak for the nonlinear system is higher than for the linear system. Nonlinear voltage response amplifications in this type of system are dependent on the time scale separation between the participating variables. For smaller values of 1 this nonlinear amplification is reduced and although the system is nonlinear, it behaves quasi-linearly Rotstein (2014Rotstein ( , 2015. The nonlinear amplification discussed above is particularly stronger for the transient dynamics (initial upstroke) than for the steady-state response, and therefore it is expected to have consequences for the nonlinear responses of nonlinear models to PWC inputs with large enough values of D (Fig. 6). We use a PWC input with a (non-random) equispaced distribution of constant piece amplitudes in the range [0, 2] multiplied by D. For small enough values of D ( Fig. 6-a1) the trajectory remains within the linear regime and therefore the responses to the linear and the PWL models are almost identical (blue and red). As D increases, the nonlinear amplification becomes stronger (Fig. 6-b1 and -c1, blue and red). This is accompanied by similar changes in the Z profile (not shown). As discussed above, these mechanisms are dependent on the time scale separation between the participating variables, determined by 1 and the nonlinear voltage response amplification is reduced for smaller values of 1 (see Fig 6-a2, -b2, and -c2 for 1 =10). Lowering the parameter 1 makes the system more linear and reduces the amplification.

Oscillation amplification and attenuation: currentversus conductance-based responses to synaptic-like PWC inputs
The oscillatory dynamics considered above emerge in response to additive PWC current inputs. However, the synaptic inputs received by neurons are multiplicative and conductance-based as described by the model (6)-(7).  (7) and (5) with S substituted by I. We used the following parameter values: C = 1 , g L = 0.25 , g 1 = 0.25 , 1 = 100 (same as in Figs. 1 and 2-b), v c = 1 and g c = 0.1 From the phase-plane diagram in Fig. 5-c we see that increasing values of I (replacing S in the model) reduces the nonlinearity of the V-nullcline (dashed-red) and increases (in absolute value) its slope. Both phenomena oppose the voltage response amplification (blue) and the overshoot becomes much less prominent. The triangular region (bounded by the V-axis, the displaced V-nullcline (dashed-red) and the w-nullcline (green)) is reduced in size as compared to the current-based inputs (panel b) and therefore the response is reduced in amplitude. Moreover, because the displaced V-nullcline in panel c is more vertical than the baseline V-nullcline, the size of the voltage overshoot in response to constant inputs is reduced and, in this sense, the voltage responses become quasi-1D. As a consequence, the initial portion of the transient responses to abrupt changes in input is reduced in size. This translates into the oscillatory response to PWC inputs, which is also attenuated and the resonant peak disappears or is significantly reduced (Fig. 6-c1, green).

Emergence of variability in response to piecewise constant inputs with normally distributed amplitudes
Here we address the relationship between the transient dynamic properties of individual cells (autonomous transient dynamics) and the variability of their responses to piecewise constant (PWC) input functions with pieces of the same duration and randomly distributed amplitudes. We primarily use Gaussian distributions and a relatively large number of pieces (with a small duration each). This allows us to separate the effects of the autonomous transient dynamics activated by the transition between pieces with variable size-steps from the variability of the piece duration, which is chosen to be in a range much lower than the cell's intrinsic (natural) and resonant frequencies.
In the next Section, we show that amplitude randomness is not needed, but only the arbitrary order of the amplitude of the constant pieces.

Variability in response to piecewise constant inputs emerges from the transient response properties of the autonomous system to step constant inputs
The voltage response of cells to PWC inputs consists of a sequence of transient behaviors (initiated immediately after the transition between two constant pieces), reflecting the autonomous transient dynamics in response to a set of initial conditions for the participating variables, followed by an approximation to the steady-state for each input piece. The latter and part of the former may be absent if the duration Δ of each constant piece is not long enough. As discussed above, Fig. 1-c1 to -e1 illustrate three qualitatively different ways in which cells respond to a sequence of increasing step function inputs according to whether they have a stable node (N, 2D), a stable focus (F, 2D) or they are passive (P, 1D). For illustrative purposes, the value of Δ was chosen to be large enough so as to show both dynamic components of the response to each linear piece. In subsequent simulations, Δ will be chosen to be much smaller so as to capture realistic situations (as discussed above), and therefore the voltage response will capture only the initial portions of the autonomous transient dynamics.
Passive (1D) cells exhibit a monotonic behavior towards the steady-state in response to each input piece. For the parameters chosen, the transient increase is relatively fast (Fig. 1-c1). 2D cells, in addition, can display overshoots ( Fig. 1-d1, N) and damped oscillations ( Fig. 1-e1, F) as they approach the steady-state in response to each input piece. The peak amplitudes of V in response to each input piece correspond to the transient component of the response and depend on the initial conditions in that regime ( Fig. 1-b), which in turn depend on the values both V and particularly w reach at the end of the previous regime. Because of this sensitivity of the transient responses to initial conditions, an order rearrangement of the constant pieces produces different response patterns ( Fig. 1-c2 to -e2). This variability, which is reflected in the peak and trough values of the "disordered" patterns (c2, d2, e2) as compared to the "ordered" ones (c1, d1, e1) is due to the differences in the initial conditions in each regime as v transitions between pieces with different amplitudes. The variability among patterns generated by different permutations of the set { k } N k=1 is inherited from this principle since the initial conditions for a given input piece depends on the "last" values of the participating variables in response to the previous piece. The variability is stronger for the 2D cells ( Fig. 1-d and 1-e) than to for the 1D cell ( Fig. 1-c1) since the sensitivity to initial conditions is weaker for the former than for the latter. The addition of a second variable (dimension) typically adds sensitivity to initial conditions. , N = 200 , Δ = 5 ) for a passive cell (P), a cell having a node (N) and a cell having a focus (F). The average < P > and the P ,step curves, superimposed to P , show that the average across trials and the response to the reference pattern (with minimal variability) are not identical. In this figure they are relatively close.

The response variability to piecewise constant inputs with normally distributed amplitudes depends on the cell's autonomous transient dynamic complexity
The variability of P across trials, Var(P ) is larger for the F-cell than for the N-cell and there seems to be no significant difference between the variability of N-and P-cells.
(The latter two are comparable since they have the same value of g L .) However, the qualitative differences between N-and F-cells are accompanied by different strengths in the responses to step-constant inputs (e.g., Fig. 1), and one has to account for these differences in order to isolate the effects of the type of transient dynamics. We note that these parameter values were chosen so that the two cells (N and F) have resonance in the same frequency band. Figure 4-b1 to -d1 show the variances of P (Var(P ) ) and Fig. 4-b2 to -d2 show the normalized variances (VarN(P ) ) computed as these variances divided by the peak of the unforced cells' response to a step-constant input of amplitude 1. Figure 4-b show that both Var(P ) and VarN(P ) are larger for the F-cell than for the N-cell. Figure 4-c shows that both Var(P ) and VarN(P ) are slightly larger for the N-cell than for the corresponding P-cell (obtained by making g 1 = 0 ). In contrast to this, Fig. 4-d shows that Var(P ) and VarN(P ) have different relative values for the F-cell and the corresponding P-cell; Var(P ) is larger for the P-than for the F-cell, but VarN(P ) is larger for the F-than for the P-cell. The values of these quantities in both cases are significantly larger for the F-/P-cells than for the N-/P-cells. Together, these results suggest that under the constraints imposed by the two cells having resonance in the same frequency band, the voltage response of the F-cell is more variable than the voltage response of the N-cell, and the larger variability of the F-/P-cells as compared to the N-/P-cells is due to a smaller value of g L , which in turn indicates a stronger amplification of the transient voltage response to step-constant inputs. However, the larger variabilities cannot be attributed to these differences in the voltage response amplitudes since they persist after the variances have been normalized.
A similar result is obtained when one lifts the resonance constraint on the N-and F-cells. Figure 7-b shows the effects of changes in g 1 on Var(P ) and VarN(P ) for fixed values of g L . As g 1 increases the transient dynamics of the cell transitions from P (green) to N (blue, f res = 7 ) to F (red, f res = 14 , f nat = 12.3 ) (Fig. 7-a). Var(P ) is slightly larger for the Fthan for the N-cell and larger than these two for the P-cell. In contrast, VarN(P ) is much larger for the F-than for the N-cell, and VarN(P ) for the P-cell is comparable to that for the N-cell.

The response variability along time and across trials depends on the levels of complexity of the autonomous transient dynamics
If the voltage response variability to PWC inputs with randomly distributed amplitudes primarily depends on the receiving cell's autonomous transient dynamics, then one expects the variability to be higher, the faster the response of the individual cells to step-constant inputs. As this response becomes faster, then it is easier for V to reach values further away from the mean. The parameter g L is the ideal one to test these ideas since it determines the time constant of the V equation. The smaller G L , the faster the response. Figure 7-c shows that Var(P ) is higher for g L = 0.1 than for g L = 0.2 and this remains true for various representative values of g 1 . The dependence of Var(P ) with g 1 is more complex and less clear. For g L = 0.2 , Var(P ) increases with g 1 , while for g L = 0.1 , Var(P ) first decreases and then increases with increasing values of g 1 . For Δ = 1 these dependences are less well separated ( Fig. 7-c2) when compared with Δ = 5 (Fig 7-c1).
For fixed values of g L and g 1 , the time constant 1 associated to the recovery variable w controls the time separation between the variables v and w. The larger 1 (the slower w), the faster the autonomous transient response. Because of this stronger sensitivity to initial conditions, the variability is expected to be larger, which is confirmed by Fig. 7-d. Together, these results show how the variability of the responses to PWC input functions with normally distributed amplitudes are controlled by the transient dynamics of the autonomous responses to piecewise constant inputs. This variability emerges as the input pieces within the same set are permuted for different trials. The differences among the responses to the same piece (piece with the same amplitude) across trials are due to differences in the initial conditions relative to these pieces caused primarily by the varying amplitudes of the preceding pieces. . The piecewise constant inputs I have Δ = 5 , N = 200 . a V traces for the node, focus and passive cells in response to a constant input of amplitude 1. b Effects of changes of the linearized resonant conductance g 1 on the variances Var(P ) and normalized variances VarN(P ) of the peaks-and-troughs patterns. Normalized variance computed as the variance (b1) divided by the peak of the unforced cells' response to a step-constant input of amplitude 1 (a). b1. Var(P ). b2. VarN(P ). We used the following parameter values: C = 1 , g L = 0.1 , 1 = 100 , g 1 = 0.2 (node), g 1 = 0.8 (focus) and g 1 = 0 (passive). c Effects of changes of the linearized conductances g L and g 1 on the variances Var(P ). c1. Δ = 5 . c2. Δ = 1 . We used the following additional parameter values: C = 1 and 1 = 100 . d Effects of changes of 1 on the variances Var(P ). We used Δ = 1 . We used the following additional parameter values: C = 1 , g L = 0.1 and g 1 = 2 . The piecewise constant inputs I have Δ = 1 , N = 1000 . The piecewise linear model is given by (3)-(4). Left and middle. The color dots indicate the peaks-and-troughs patterns for all trials reorganized so that the corresponding values of k from which they originate are ordered in a monotonically increasing manner. All dots for a given piece correspond to the same value of k . < P > (blue) is the mean value of the reordered peaks-and-troughs patterns for each linear piece. < P ,step > (green) is the peaks-and-troughs pattern corresponding to the (ordered) input function I ,step . Right. Comparison of Var(P ) between the linear (red) and piecewise linear model (blue). a We used the following parameter values: C = 1 , g L = 0.5 , g 1 = 0 , g c = 0.1 and v c = 0.5 . b We used the following parameter values: C = 1 , g L = 0.5 ,

The peak-and-trough profiles are able to capture the nonlinear properties of the target cells
In our previous discussion, we have used linear models as the receiving cells to the PWC inputs. We have also shown that the presence of certain types of nonlinearities affects not only the vector field but also the effective time constants of the cell, which in turn affect their autonomous transient dynamics. We reasoned that these types of nonlinearities may also affect the variability of the cells' response to PWC inputs. To test these ideas we used the piecewise-linear (PWL) model (3)- (5) where the v-nullcline (right hand side of the first equation equal to zero) is broken into two linear pieces at v = v c . We chose v c > 0 to be within the range of values of the v response to the PWC input and g L > g c . Effectively, the cell's membrane time constant transitions from g −1 L to g −1 c at v = v c . Therefore one should expect the variability pattern to be larger. Since the PWL model is asymmetric with respect to the equilibrium of the autonomous cell, then one should expect the variability pattern to be asymmetric with respect to the responses mean. Figure 8 (rows 1 and 2) illustrate these ideas. Figure 8 (row 3) summarize these results and also shows that the variability for the nonlinear models shows a strong dependence on the amplitude of the input constant pieces, while the variability pattern for the linear model is flatter. This is also expected since the larger the input amplitude, the more likely the response to reach values beyond v c . However, the presence of nonlinearities in the model does not necessarily guarantee a nonlinear response to external inputs. The latter strongly depends on the time scale separation between the participating variables as it occurs in other types of responses (e.g., to oscillatory inputs) where the mechanisms depend on the cell's autonomous transient dynamics Rotstein (2015). The stronger the time scale separation (the larger 1 ), the stronger the nonlinear amplification of the voltage response to time-dependent inputs. Consistent with this, Fig. 8-c1 shows that the peaksand-troughs profiles for the nonlinear system and 1 = 10 is Fig. 9 Piecewise constant inputs with arbitrarily ordered, but equispaced (non-randomly) distributed amplitudes capture the oscillatory properties of the target cells and the variability resulting from the autonomous transient dynamics. The piecewise constant inputs I and I ,step have constant pieces with equispaced distributed amplitudes (equal amplitude differences) in the interval [−2, 2] with Δ = 1 (total time = 1,000,000 ms for row 1 and total time = 1000 ms for each trial for row 2. Similar results (with less resolution) are obtained for a smaller number of pieces for row 1 (total time = 10,000 ms). Trials consist of different permutations (1000) of the same set of constant pieces { k } N k=1 . a g 1 = 0 . b g 1 = 0.2 . c g 1 = 1 . We used the additional parameter values: C = 1 , g L = 0.1 and 1 = 100 . Top. Impedance amplitude (Z) and (rescaled) power spectra density (PSD) profiles for the sample V trace for the responses to I (random) and I ,step (ordered). The PSD profiles were rescaled so that the maxima of the PSD and Z profiles coincide. Bottom. The color dots indicate the peaks-and-troughs patterns for all trials reorganized so that the corresponding values of k from which they originate are ordered in a monotonically increasing manner. All dots for a given piece correspond to the same value of k . < P > (blue) is the mean value of the reordered peaks-and-troughs patterns for each linear piece. < P ,step > (green) is the peaks-and-troughs pattern corresponding to the (ordered) input function I ,step more symmetric and spans a smaller range than the profile for 1 = 100 (Fig. 8-b1), while the peaks-and-troughs profiles for the corresponding linear systems are similar (Fig. 8-b2 and -c2). In addition, Var(P ) is smaller and flatter for 1 = 10 ( Fig. 8-c1) than for 1 = 100 (Fig. 8-c2).

The variability in the response to piecewise constant inputs requires constant pieces with arbitrarily distributed amplitudes, but not amplitude randomness
So far we have used PWC input functions with normally distributed amplitudes (Fig. 3-a) motivated by the fact that in the limit of Δ → 0 , I approaches white noise. Earlier in the paper we have shown that the ability of PWC inputs to capture the resonant properties of a given cell depends on the multiple ways the receiving cell's autonomous transient dynamics are activated by the different initial conditions resulting from the transition between different linear pieces, while randomness does not play any significant role in this process. Here we extend these ideas to the cells' response variability. To this end, we use fully deterministic distributions of amplitudes within some range as described in Sect. 2.2.2 leaving the choice of the subset of permutations (number of trials) for each protocol as the only source of uncertainty in the input signal. Figures 9, S2 and S3 show that for the three types of deterministic distributions presented in Fig. 3, the PWC inputs uncover the voltage oscillatory properties of receiving cells (top panels) and show the same type of variability as the responses to the normally distributed PWC inputs discussed above (bottom panels). The differences are in the details. The main observed differences (by inspection) are in the responses to the P cells (panels a2) where the VarP profiles are denser for the equispaced than for the bell-shapedlike PWCs and differences in the values of D for the latter are reflected in the peaks-and-troughs distributions in the VarP profiles. A more detailed analysis is beyond the scope of this paper.

The activation of the autonomous transient dynamics by the piecewise constant inputs with small 1 evokes the steady-state oscillatory properties of the receiving cells, and not their transient oscillatory properties
Along the previous sections, we have shown a number of examples where both the impedance amplitude Z(f) and the voltage PSD profiles of cells receiving PWC inputs with small durations Δ exhibit resonance, independently of whether the corresponding autonomous systems are a node (no intrinsic oscillations) or a focus (intrinsic damped oscillations). The impedance amplitude profile Z(f) of a system captures its steady-state response to oscillatory inputs. For linear systems, it is the magnitude of the complex-valued coefficient of the particular solution to the system (see Appendix 2) forced by a sinusoidal function. In this calculation, the transient component of the solution to the forced system (the solution to the corresponding homogeneous system) is ignored. The resonance frequency f res , the peak frequency of Z(f), is different from the natural frequency f nat (the frequency of the transient damped oscillations). The latter is computed from the eigenvalues for the homogeneous system (see Appendices 1 and 2, and also Richardson et al. (2003); Rotstein and Nadim (2014)), which controls the autonomous transient dynamics.
On the other hand, the dynamic mechanisms of generation of resonance involve the autonomous transient dynamics as uncovered by using dynamical systems tools Rotstein (2014Rotstein ( , 2015. Most significantly, the interplay of the input frequency (whose inverse is a measure of the input time scale) and the cell's intrinsic time scale determines the direction of motion of the response limit cycle trajectory in the phase-space. For the resonance frequency, this interaction of time scales produces the limit cycle trajectory with the maximal amplitude in the v direction. This mechanism does not rely on the cell's eigenvalues and eigenvectors (autonomous steady-state dynamics) and it underlies the resonant responses of both N-and F-cells.
Because white noise has a constant PSD, both the Z and the V PSD profiles have the same frequency-dependent properties. This extends to PWC inputs with small enough values of Δ (e.g., Figs. 2a2 to -c2, 9-and S2-to S3-a1 to -b1), but not necessarily for large values of Δ (e.g., Fig. 2a2 to -c2, insets). Figure 9-b1 (see also Figs. S2-to S3-b1) shows that both the Z and the V PSD profiles peak at f res when f nat = 0 (the receiving cell is a node). For these parameter values there is a clear separation between the cell's response to the oscillatory input and the autonomous intrinsic dynamics (captured by f nat ). In Fig. 9-c1 (see also Figs. S2-c1 to S3-c1) the V PSD profile seems to be superimposed to the Z profile, but because f res and f nat are very close, it is not clear whether the V PSD profile peaks at f res , f nat or in between. For the parameter values in Fig. 10-a1 the separation between f res and f nat is bigger and the V PSD profile clearly peaks around f res and far away from f nat for Δ = 1.
In contrast, Fig. 10-b show more complex responses for Δ = 50 . The duration of the constant pieces Δ determines an input time scale that interacts with the cell's intrinsic time scale to produce the V response. This is true for all values of Δ , but for small values of Δ the corresponding frequencies ( 1000∕Δ ) are very large, away from the resonant frequency band, and the response time to each constant input piece very small. As Δ increases, the interaction between time scales is stronger and felt for the lower frequencies.
For example, in Fig. 10-b for Δ = 50 , the V PSDs show resonances occurring at multiples of 1000∕Δ = 20 . The shape of these resonant patterns reflects the properties of the autonomous transient dynamics rather than the properties of the corresponding Z profiles. In Fig. 10-b1, the sequence of PSD peaks decreases, while in Fig. 10-b2, the sequence of PSD peaks increases and the maximum of these peaks (save the maximum at f = 0 ) occurs roughly at f nat . This reflects the stronger oscillatory responses evoked by each constant piece due to their autonomous transient dynamics.
In summary, the autonomous transient dynamics plays a role in shaping the V response PSD patterns, but these patterns transition from reflecting the stationary properties of the V response to oscillatory inputs, captured by the Z profile, to the intrinsic properties of the receiving cell as Δ increases.

Discussion
Neuronal systems are subject to fluctuations either intrinsically or externally Deco et al. (2009);Faisal et al. (2008); Shalinsky et al. (2002); White and Haas (2001); Dorval and White (2005); DeFelice (1981); Chow and White (1996); ; Schneidman et al. (1998);Sigworth (1980); Middleton et al. (2003); Richardson (2008), which have been modeled as random Gaussian white or colored noise Destexhe and Rudolph-Lilith (2012). Cells subject to variable inputs have been shown to exhibit a number of non-expected behaviors (see Introduction for more details and references), including oscillatory voltage responses and resonances, and they also exhibit variability across trials. These phenomena result from the dynamic, often non-trivial interplay of the cell's intrinsic properties and the properties of the noise. While variable inputs are ubiquitous, in particular in neuronal systems Laing and Lord (2010); Destexhe and Rudolph-Lilith (2012);Faisal et al. (2008), several aspects of the dynamic mechanisms governing these interactions are not well understood. In particular, it is not well understood what aspects of the cell's intrinsic dynamics play a role in each of these interactions, what properties of the variability are necessary, if at all, and how the two are integrated to produce these behaviors. Linear systems receiving Gaussian white noise can be solved analytically Uhlenbeck and Ornstein (1930); Gardiner (1985); Risken (1989), but, in general, the resulting formulas provide limited mechanistic information. Nonlinear systems are generally not amenable to analytical solutions. The investigation of dynamical systems primarily focuses on stationary solutions Strogatz (1994); Izhikevich (2006); Guckenheimer and Holmes (1983) and the properties of the transient behavior are less explored or even overlooked. Attractor networks are a key concept in neural computation Hopfield (1982); Amit (1989); Knierim and Zhang (2012); Samsonovich and McNaughton (1997); Redish et al. (1996); Boucheny et al. (2005); Zhang (1996) and have been proposed to capture a variety of cognitive process (e.g., working memory, navigation). Recently, transient dynamics have been argued to be important for neuronal processing and the encoding of information Rabinovich et al. (2008); Nachstedt and Tetzlaff (2016); Bondanelli and Ostojic (2020); Robinson and Harsch (2002); Rabinovich and Varona (2011);Mazor and Laurent (2005);Stopfer et al. (1997).
In previous work we showed that the cell's autonomous transient dynamics (the cell's transient voltage response to constant inputs, which is uncovered by the voltage response to abrupt changes between two constant inputs) play a significant role in the mechanistic explanation of the generation of resonance in linear and non-linear systems in response to oscillatory inputs Rotstein (2014Rotstein ( , 2015. Even for linear systems for which the impedance amplitude and phase profiles can be computed analytically, the resulting formulas do not explain why and under what circumstances resonance emerges in systems that do not display intrinsic oscillations (i.e., systems having only real eigenvalues). We developed a dynamical systems approach to show that the cell's response to oscillatory inputs can be thought of as a continuous sequence of transient voltage responses to constant inputs equal to the value of the input at a given time. The response trajectory for each input frequency tracks the cyclic motion of the voltage nullcline in the extended phaseplane diagram (technically, the projection of a surface in the three-dimensional space) and this evolution is governed by the properties of the autonomous transient dynamics. The resonant frequency corresponds to the optimal balance between the effective time scale at which the system transiently reacts to constant inputs and the input frequency. This could be captured by approximating the sinusoidal inputs by piecewise constant (PWC) functions with short-duration pieces and computing the sequence of (property joined) voltage response trajectories to these PWC inputs. Because in general these trajectories do not reach a close enough vicinity of the equilibrium for each constant piece, knowledge of the steady-state properties of the unforced target cell (eigenvalues and eigenvectors) does not inform the generation of resonance. We reasoned that similar ideas could shed light on the issues raised in the previous paragraph.
In this paper we set out to investigate the mechanisms of generation of oscillations and variability in response to PWC inputs I with short-duration constant pieces and variable amplitudes. For each constant piece, the autonomous transient dynamics are able to develop, but the voltage response does not reach a close vicinity of the corresponding steady-state. For randomly distributed amplitudes , the PWC inputs provide an approximation to Gaussian white noise Allen et al. (1998). However, although our results have implications for systems subject to random Gaussian noise, we are not making precise statements about the transition from the responses to the deterministic PWC inputs to the responses to stochastic Gaussian noise inputs as the duration of the constant pieces approaches zero. This requires more research and is beyond the scope of this paper.
We showed that oscillatory behavior can be generated in response to additive PWC inputs I in both linear and nonlinear systems, independently of whether the stable equilibria of the unperturbed systems are foci or nodes. These results persisted when the amplitudes of the constant pieces were not randomly distributed, but chosen following a deterministic rule. For linear systems, this captures earlier results described in Gardiner (1985) for Gaussian white noise (see also Thomas and Lindner (2019) for two-dimensional linear cells with foci). However, the dynamic mechanisms of generation of oscillations and their dependence on the autonomous transient dynamics, in particular for cells having a node were not known. To understand the effects of the nonlinearities on the oscillatory voltage response we used a piecewise linear model (PWL) developed in Rotstein (2014) that mimics the voltage-dependences present in neuronal models in the voltage nullclines (for the current-balance equation). We showed that the oscillatory voltage response is amplified by the nonlinearities. This nonlinear amplification is consistent with previous work on resonance (for sinusoidal inputs) Rotstein (2014) and can be also explained in terms of the interaction between autonomous transient dynamics and the input properties. The PWL model has two regimes, the linear regime, in the vicinity of the fixed-point, and the nonlinear regime for voltage values higher than the "breaking point" determining the boundary between the two PWL input pieces (see Fig. 5-b). The nonlinear regime is in itself linear with respect to a virtual fixed-point, which is different from the actual fixed-point (the origin in Fig. 5b). For each constant input piece for which the response trajectory crosses to the nonlinear regime, this trajectory is able to reach larger voltage values as compared to the linear model for the same input. This is because the slope of the voltage nullcline in the nonlinear regime is larger (less negative) than the slope in the linear regime, and therefore the trajectory for transient voltage response to the corresponding constant input is able to jump to higher values as compared to the linear regime. This clearly depends on the properties of the PWL model, in particular on the slope of the "broken" linear piece, which we chose to be larger (less negative) than the linear piece in the linear regime (see Fig. 5-b). Other slopes may lead to attenuation of the voltage response. Following similar ideas, we showed that the voltage response to multiplicative conductance-based synaptic inputs is attenuated as compared to the underlying model with additive noise. In this case, the "broken" piece is more negative than the linear one (see Fig. 5-c). This has mechanistic implications for the experimental results discussed in Fernandez and White (2008), comparing the voltage responses of current-and conductance-based synaptic fluctuations. A more detailed explanation requires the use of high-frequency Poisson-distributed synaptic-like inputs. We analyze this in more detail in a companion paper Pena and Rotstein (2022).
The effect of the autonomous transient dynamics on the response variability to I is already apparent in the simple illustrative example in Fig. 1-c to -e. The voltage response variability emerges as the result of the transition between regimes corresponding to different constant pieces. The variability is minimal for the response to I ,step ( Fig. 1-c1 to -e1) where the constant pieces are arranged by increasing amplitude and the variability increases as the constant pieces are permuted ( Fig. 1-c2 to -e2) since the transient voltage response for each constant piece regime depends on the initial conditions for that regime. These, in turn, depend on the value of the participating variables at the end of the previous regime, which is different for different trials (different permutations of the order of the constant pieces). To clarify these ideas we designed protocols where the set of amplitudes was the same for all trials and each trial used a different permutation of this set. We then rearranged the voltage response peak-and-trough profiles P in increasing order of input amplitudes . In this way, we were able to compare the different voltage responses to the same input amplitude across trials. We used both randomly distributed amplitudes and deterministic distribution of amplitudes. For each input amplitude k , the variability in P k across trials was the result of the multiple different ways the target cell reacts to the input k . Again, this is due to the differences in the initial conditions for that regime, which in turn depend on the values of the participating variables at the end of the previous regime. This variability was strongly affected by the properties of the target cells and their ability to produce overshoots and damped oscillations in response to abrupt changes in constant inputs, which in turn depends on the model parameters, which constrain these responses. This interpretation is strengthened by the fact that the variability patterns remained qualitatively the same when we relaxed the requirement that the amplitudes are randomly distributed and we used permutations of the same deterministic distribution for each protocol. Although it was not the full overshoot or the damped oscillations that affected the voltage response, but rather the initial portion to these responses, the particular dependence on the model parameters remains. In other words, cells that exhibit stronger autonomous transient dynamics show more variability. In this sense, this is consistent with the results for one-dimensional OU processes and it is intuitive for higher-dimensional OU processes, however, it is not directly clear from the complex covariance formulas. For nonlinear processes, analytic computations are not possible and therefore the method we developed allows to make predictions based on the knowledge of the autonomous transient dynamics and the structure of the input. Our results further predict that the variability would be reduced if the input functions involve gradual rather than abrupt transitions.
These predictions and the predictions about the oscillatory properties of cells in response to PWC inputs could be tested experimentally in vitro using current clamp techniques or in vivo using optogenetic tools Zhang et al. (2007); Deisseroth (2011); Bernstein and Boyden (2012). In the first case, variability in response to PWC inputs is expected to be well correlated with appropriate metrics for the autonomous intrinsic dynamics measured in response to constant inputs. Obtaining these correlations in the second case could be more challenging, but still possible by using the appropriate tools to measure the response membrane potentials.
Time-dependent inputs can be approximated by a discrete sequence of short-duration constant pieces (e.g., Fig. 1-a) and therefore the sequence of the properly joined transient solutions to these constant pieces (Fig. 1-c to -e) provides an approximation to the system's voltage response to the original (time-dependent) inputs. The short duration of the constant pieces does not allow for the corresponding steady-state solutions to develop and therefore the information about the input is encoded by the transient trajectories. Complex odors (different types of stimuli at different concentrations) and their representations have been proposed to operate in a similar way Rabinovich et al. (2008);Stopfer et al. (1997). A potentially more involved experimental test or application of the ideas presented in this paper is the presentation of complex odors consisting of sequences of different permutations of a number of "basic odors" Rabinovich et al. (2008);Stopfer et al. (1997). Variability in the responses is expected across trials and an increase in the variability patterns is expected as the "intensity" of these odors increase. However, odor representations in the olfactory research is a network phenomenon and therefore additional research is needed to extend our results to include network effects.
Our results have also implications for the understanding of the emergence of variability in neuronal systems  and dynamical systems in general, and in the absence of noise. In the simple systems we used here, variability emerges as the result of the interaction between the system's autonomous transient dynamics and a set of deterministic inputs where the only source of uncertainty is the choice of the subset of all possible PWC functions with the same set of constant piece amplitudes. Key for this variability to emerge are the abrupt transitions between constant piece regimes and the sensitivity to initial conditions of the autonomous transient dynamics. In higher-order networks, this can be evoked by the network internal dynamics without the need of external inputs. Therefore our results have implications for the encoding of information in general and for these resulting from expansions or quenching of variability Churchland et al. (2010); Ito et al. (2020). More research is needed to establish in these networks how variability is affected by the node and connectivity parameters. Finally, if variability encodes the autonomous transient dynamics of the target cells, then decoding methods should be able to infer these dynamics. These also require more research.

Intrinsic and resonant oscillatory properties of 2D linear systems
Consider where a, b, c and d are constants, = 2 f ∕1000 > 0 is the input frequency and A in ≥ 0 is the input amplitude. The prime sign represents the derivative with respect to t. The units of t are ms and the units of f are Hz.

Intrinsic oscillations
The characteristic polynomial for the corresponding homogeneous system ( A in = 0 ) is given by The eigenvalues are given by and the natural (intrinsic) frequency of the (damped) oscillations (in Hz if t has units of ms) is given by assuming (a − d) 2 + 4bc < 0.

Resonance and the impedance amplitude profile
The impedance amplitude profile Z( ) for system (10)-(11) is the magnitude (10) x � = a x + b y + A in e i t , y � = c x + d y (11) r 2 − (a + d) r + (a d − b c) = 0.
(12) r 1,2 = a + d ± √ (a − d) 2 + 4bc 2 , (13) f nat = √ −(a − d) 2 − 4bc 4 1000 of the complex valued coefficient of the particular solution to the system For 1D system, these quantities are given, respectively, by and The resonance frequency f res (in Hz if t has units of ms) is the frequency at which Z reaches its maximum

Response to constant inputs
The equilibrium solution to system (10) for a constant input A in (i.e., = 0 ) is given by

One-dimensional OU process
The 1D OU process Uhlenbeck and Ornstein (1930) is described by the following linear stochastic differential equation where a > 0 , I and are constants and (t) is zero-mean and -correlated Gaussian white noise. The parameter a is the inverse of the time constant and measure the strength by which the system reacts to perturbations. The parameter measures the intensity of the noise. The quotient I∕ is the asymptotic mean.
Using standard methods Risken (1989); Gardiner (1985) one can compute the solution satisfying X(0) = x 0 , which is the sum of a deterministic function with the form (24) and an integral of a deterministic function with respect to a Wiener process. The solution is normally distributed with mean and variance given, respectively by and

Higher-dimensional OU process
The multivariate OU process Uhlenbeck and Ornstein (1930) is described by the following linear stochastic differential equation where X is an n-dimensional vector, A is an n × n matrix, B is an n-dimensional vector, Σ is an n × m matrix and H is a vector of independent zero-mean and -correlated Gaussian white noise components. Using standard methods Risken (1989); Gardiner (1985) one can compute the solution satisfying X(0) = x 0 . The solution is normally distributed. The mean is given by (1 − e −2at ).