{'fraction_non_alphanumeric': 0.03327168552364375, 'fraction_numerical': 0.009501015375689006, 'mean_word_length': 4.657878539187526, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Estimation of Synaptic Activity during Neuronal Oscillations

In the study of brain connectivity, an accessible and convenient way to unveil local functional structures is to infer the time trace of synaptic conductances received by a neuron by using exclusively information about its membrane potential (or voltage). Mathematically speaking, it constitutes a challenging inverse problem: it consists in inferring time-dependent parameters (synaptic conductances) departing from the solutions of a dynamical system that models the neuron's membrane voltage. Several solutions have been proposed to perform these estimations when the neuron uctuates mildly within the subthreshold regime, but very few methods exist for the spiking regime since large amplitude oscillations (revealing the activation of complex nonlinear dynamics) hinder the adaptability of subthreshold-based computational strategies (mostly linear). In a previous work, we presented a mathematical proof-of-concept that exploits the analytical knowledge of the period function of the model. Inspired by the relevance of the period function, in this paper we generalize it by providing a computational strategy that can potentially adapt to a variety of models as well as to experimental data. We base our proposal on the frequency versus synaptic conductance curve (f − gsyn), derived from an analytical study of a base model, to infer the actual synaptic conductance from the interspike intervals of the recorded voltage trace. Our results show that, when the conductances do not change abruptly on a time-scale smaller than the mean interspike interval, the time course of the synaptic conductances is well estimated. When no base model can be cast to the data, our strategy can be applied provided that a suitable f − gsyn table can be experimentally constructed. Altogether, this work opens new avenues to unveil local brain connectivity in spiking (nonlinear) regimes.


Introduction
Estimating non-measurable connectivity parameters from experimental measurements is an important challenge to understand brain function. Brain activity is mainly driven by the intrinsic dynamics of neurons and their mutual connections, which are essentially regulated by synapses. A thorough knowledge of individual membrane potentials and synaptic terms is, therefore, crucial to understand how the brain operates. However, we are far from this utopia: current recording methods can only retrieve partial information, such as the time evolution of membrane potential of one or a few neurons. Unfortunately, other quantities that shape the connectivity, like synaptic conductances, are impossible to be directly recorded in experiments.
The time course of synaptic conductances arriving at a single cell not only provides local information about its functional connectivity, but it often constitutes an indicator to validate more general theories such as obtaining information about the architecture of the cortex, discerning between phase insensitive or spatial phase selective coupling [27], unveiling whether the integrative properties of cortical neurons are altered under high-conductance states [13], as well as studying many neuronal problems related to the balance between excitation and inhibition [39]. Unfortunately, there are no direct methods to estimate synaptic conductances and, thus, they have to be deduced from experimentally observable quantities, which leads in turn tosince no direct methods are known to estimate synaptic conductances, they have to be deduced from experimentally observable quantities. This leads to the need of using inverse mathematical methods. As many inverse problems, it has been proved to be an arduous analytical and computational task for which no general procedures are yet available.
Even though all these strategies provide circumstantial solutions, it is also well known that they can only be applied in the subthreshold regime, where ionic currents are negligible (see [18] and [21]), a situation for which simplied models (close to linear) are representative enough to account for the voltage dynamics. However, it has also been demonstrated that these linear procedures can also provide wrong estimations when ionic currents are active in the subthreshold regime causing non-spikes in the membrane potential of the postsynaptic neuron. This fact motivates the introduction of new strategies to estimate conductances under these non-linear paradigms (see [37] and [36]). More recently, in [9], the authors identify that another artifact that can provide wrong estimations is to assume that the synaptic conductances are voltage-independent; to avoid this problem, they provide a new estimation procedure, also in the linear framework. Despite of the many eorts in the literature to get rid of artifacts that can contaminate the estimation of conductances (inuence of noise, activity of subthreshold ionic currents, variability across trials,. . . ), one of the most resilient troubles are the nonlinearities that the f − I curve exhibits in the presence of spiking mechanisms.
On the other hand, a line of research is devoted to look for possible strategies to estimate conductances in the spiking regime. For instance, in [8], the authors use a ring clamp technique where a description of the spike is given in terms of the membrane potential threshold and the amplitude of the spike. These traces are approximated by calibrating the membrane potential using dynamic clamp in in vitro experiments of isolated neurons and a quadratic polynomial approximation. More recently, in [17] we introduce a mathematical approach to estimate conductances by characterizing the period of the membrane potential, and in [24] estimations of the conductances of spike trains of the membrane potential are performed.
In [17] we present a rst proof-of-concept to estimate synaptic conductances in the spiking regime based on information about the f − I curve of the neuron. The idea was exemplied with the McKean model, a piecewise linear system that allows for a quasi-analytic expression of the f − I curve. More specically, we approximated the ring period, T , as a function of the parameters of the system and, in particular, of the synaptic conductance parameter, g syn . Therefore, for a given interspike interval (ISI) of length T * , one can numerically solve the implicit equation T (g syn ) = T * to obtain an estimation, g syn,est , of the synaptic conductances during that ISI. In order to obtain the time course of the synaptic conductances, the procedure is repeated for each ISI. By interpolation, a continuous estimation of g syn,est (t) can be obtained, which is proved to be an accurate estimation of the prescribed conductances provided that they have slower variation (in time) than the membrane potential. Although the results obtained in [17] show that using the f − I to estimate synaptic conductances in spiking regimes can be a suitable strategy, there were only applied to the McKean model, a very specic and non-realistic model, which is nothing but a caricature of the Fitzhugh-Nagumo system. An improvement of this estimation procedure is to consider a dierent underlying mathematical model that could better t experimental data. On the one hand, desirable candidates would be conductance-based models as per the Hodgkin-Huxley formalism but, even though some of them (pyramidal cells [38], stellate cells [33],. . . ) have been investigated in the context of estimation of conductances in subthreshold regimes, it is extremely dicult to devise a procedure for a generic conductance-based model. On the other hand, integrate-and-re (IF) models have also been considered as good approximations to the neuronal activity of dierent type of neurons (see [15,19,31,32,12]), and are simple enough. We found that the exponential integrate-and-re (EIF) model represents an optimal balance between adaptability to experimental data and computational f − I curves, and model simplicity.
In this paper, following the ideas introduced in [17], we propose a procedure to estimate conductances under a regular spiking regime. This procedure is based on a modication of the standard frequency-input curve, which consists of replacing the input constant current by a constant synaptic conductance (f − g syn curve). We aim to estimate the global synaptic conductance time course from a membrane potential time series of a spiking neuron, for which we can have an accurate f − g syn curve. For simplicity, here we assume that all synaptic input is excitatory.
The f − g syn curve can be provided in several ways: (I) from dynamic clamp experiments, using in vitro recording; (II) approximating the real data by a base model (IIa) for which we can obtain an analytic expression for the f − g syn curve, such as the leaky or the quadratic integrate-and-re models (see [3] and [2]) or (IIb) allowing an approximate expression for the f −g syn curve, such as the McKean model (see [17]) or the exponential integrate-and-re (see Section 2.1); nally, (III) using a more detailed model that can be numerically integrated in order to construct a f − g syn curve table, that is, providing a set of period-conductance pairs extracted from numerical experiments (see Section 2.2). Once the f − g syn relationship is obtained, we have an expression of the period T = 1/f depending on the conductances (or viceversa). In cases (I) and (III), where we have a discretization of this relationship, in order to obtain a continuous function T (g syn ), one can interpolate the data using, for instance, the piecewise cubic Hermite interpolating polynomial (PCHIP). Hence, given an interspike interval T * , the corresponding conductance within this interval can be obtained by solving analytically or numerically the equation T (g syn ) = T * .
Notice that for case (I), which involves dynamic clamp, the method can be considered to extract the f − g syn curve in in vitro experiments but not in in vivo since several injected currents would be required in order to have a precise f −g syn curve, which could add both variability across trials and noise measurement. Alternatively, we provide dierent computational strategies (II-III) to obtain the f − g syn curve. Even though strategies (IIa) and (IIb) provide good approximations to the f − g syn curve, they should not be applied in general to approximate the f − g syn curve for a wider range of neurons. For this reason, the computational approach presented in Section 2.2 focuses on the idea of case (III).

Neural signal without synaptic inputs
given a neural signal with synaptic inputs, estimate g syn (t) Fit the neural signal using a base model Compute the period when different contant conductances are applied Figure 1: Scheme of the estimation procedure performed by Algorithm 1. Given a neural signal from a neuron which is not receiving synaptic inputs (gsyn(t) = 0; upper left) but exhibiting regular spiking, we t this signal using a computational model where the synaptic term is added afterwards (bottom left, where the red color stands for the tting). Next, we use this model to compute the period (or frequency) for dierent constant synaptic conductances (T (gsyn), bottom right). We nally estimate the synaptic conductance time course from a neural signal subjected to synaptic inputs (upper right).
The paper is organized as follows. In Section 2.1, we present the neuronal model we use to generate data by means of prescribed conductances and its approximation to the EIF model. In Section 2.2, we present the estimation procedure giving a pseudo-code of its implementation sketched in Figure 1. Next, in Section 3, we show the results of the estimation procedures. Finally, in Sections 4 and 5 we discuss the results obtained and present our conclusions.

Methods
In this section we present the estimation procedure we use to estimate conductances from a neuronal signal that contains high activity (Subsection 2.2). In order to prove the eciency of the procedure, we need to generate in silico data by specifying the actual conductances. These data are the input of the procedure, whose output is compared with the actual conductances. Then, we also present here the neural models we use and the prescribed synaptic drive (Subsection 2.1).

2.1
Models and data treatment 2.1.1 Neuronal models Along the paper, we consider dierent neuronal models both for the purpose of generating membrane potential traces and to use them as base models for the estimations of synaptic conductances. We keep the choice of models minimal; as we mention in other parts of the paper, we can adapt our procedure to many other models.
Pyramidal cell model We rst consider a simplied version of a model for a pyramidal neuron, see [38], where only the soma compartment has been considered. It is a spiking model driven by the interaction of sodium (I N a ) and potassium (I K ) currents. The membrane potential V follows the dierential equation where C m is the capacitance, I syn is the synaptic current, I app is the applied current, I L is the leak current, and I N a and I K are the ionic currents, which are described by equations where V x and g x , with x ∈ {L, N a, K}, represent the specic reversal potentials and maximal conductances, respectively. The variables h and n are gating variables governed by rst-order kinetics of typė where w represents either h or n; w ∞ (V ) := α w (V )/(α w (V ) + β w (V )) is the w steady-state and τ w (V ) = (α w (V ) + β w (V )) −1 is the w time decay. In our version of this model, the m-type (fast) variable is considered to be at the steady-state m = m ∞ (V ).
The oscillatory regime for this model is approximately I app ∈ (0.225653, 108.285) µA/cm 2 . By default, we take I app = 0 µA/cm 2 . (see Appendix B in [18] for a justication of this choice).
For more details about the biophysical parameters of the model, see Appendix A.
We will also explore how the random synaptic input or measurement noise aects the estimation of conductances and so, we will add a noise term to the membrane potential equation (1). For the sake of reference we present the resulting system as a separate model. Noisy pyramidal cell model Considering a white noise process η(t) having zero mean and standard deviation σ V , system (1) writes as Stellate cell model The stellate cell model is taken from [33]. The membrane potential dynamics is given by where C is the membrane capacitance, I syn the synaptic current and I app the applied current. The leakage and ion currents are given by: where V x and g x , for x ∈ {L, N a, K, N aP, h}, denote the corresponding ion reversal potentials and maximal conductances, respectively.
The gating variables w, which can be either m, h, n, p, r f or r s , follow the same rst order dynamics (2), as for the pyramidal cell model.
In absence of synaptic inputs, this model presents spikes for an applied current greater that I T = −1.46µA/cm 2 . In our case, to have a high frequency of spikes, we inject a constant current equal to I app = 8 µA/cm 2 .
For more details about the biophysical parameters of the model, see Appendix A.
Exponential Integrate-and-Fire (EIF) model with constant adaptation We chose the EIF model as representative of base models due to well-known properties of being able to shape its parameters in order to t to diverse f − I, see [15]. Therefore, it will be used during the estimation procedure of the synaptic conductances. We write it as Here, I T is the largest input current for which the neuron does not spike, V T is the corresponding membrane potential on the V − I curve, ∆ T stands for the spike slope factor (which is related to the sharpness of the spike initialization), and the rest of terms are described as in the pyramidal cell model (1).
When the membrane potential reaches the numerical threshold ϑ th , which denes the spiking time t spike , the membrane potential is reset to V reset and the integration restarts at t spike + ∆ abs , where ∆ abs is the absolute refractory time (typically chosen between 0 and 5 ms, see [16]).
To t the pyramidal cell model (1), the specic EIF parameters have been set to be I T = 0.16 µA/cm 2 , V T = −59.9 mV , ∆ T = 2.97 mV , ϑ th = −51 mV , V reset = −71 mV , and ∆ abs = 1.25 ms (see Figure 2). a b Figure 2: Approximation of the pyramidal cell model by the EIF model. Panel a shows the membrane potential trace obtained from the pyramidal cell model (black curve) given in equation (1) when no synaptic input is considered, and the approximation of the membrane potential trace obtained by the EIF model (red curve) described in equation (5). Panel b shows their respective T − gsyn curve, where T is the period of oscillation (the black curve corresponds to the pyramidal case and the red curve to the EIF model).
In this case, the interspike period can be computed as the time of ight between the reset value (V reset ) and the threshold (ϑ th ) for the voltage, adding the absolute refractory time (∆ abs ). By solving the dierential equation (5) in this range of membrane potential values, if g syn is considered constant, the interspike period is given by We note that this expression provides an implicit T − g syn relationship since I syn depends on g syn . Hence, for a specic conductance value, we can numerically solve it using a numerical method for integrals (see 2.1.3).
Unfortunately, there is no method to further simplify this expression.
Note that equation (6) can only be used if we suppose that dV /dt is far from zero, that is, if we are not in the resting regime. However, for those specic time intervals where dV /dt ≈ 0, from (5), the synaptic conductance is given by where V = V (t).

Synaptic drive
In our simulations we assume that the target neuron is stimulated through two dierent currents: (1) an external current, I app , which is a constant value that places the cell in the excitable regime, ready to induce spiking activity within a biologically plausible frequency range; and (2) the synaptic input I syn , which comes from other neurons of the network and eventually leads to the variation of the spike frequency.
In the dierent neural models, the synaptic input is always supposed to be of the form I syn = g syn (V −V syn ), where g syn = g syn (t) is the synaptic conductance trace and V syn is the synaptic reversal potential; assumed to be excitatory in our simulations.
To generate the membrane potential of the neuron, we consider specic prescribed conductances. In our in silico experiments, we inject the corresponding synaptic currents, we extract the voltage response and, pretending not to know them, we nally apply our methods to estimate the synaptic conductances. For this purpose, we use two dierent types of conductance traces: a synthetic three-frequency oscillatory drive and a more realistic trace obtained from a computational network. Note that, due to the nature of the problem we tackle here, there are no actual conductance traces available unless dynamic clamp experiments are performed.
We dene the three-frequency oscillatory drive as being A 1 = 0.0022, A 2 = 0.002, A 3 = 0.001, A 4 = 0.025, T 1 = 150, T 2 = 320, and T 3 = 50. These values have been chosen to ensure a plausible spiking frequency for such kind of neurons. The resulting conductance trace is shown in Figure 3a (black curve); observe that it elicits a membrane potential completely immersed in the spiking regime (without silent periods interspersed), see Figure 3c.
On the other hand, to explore the synaptic regime interspersed by short silent periods, we chose a more realistic synaptic conductance input, consisting of a trace with a 1 ms resolution obtained from a computational network that models the layer 4Cα of the primary visual cortex (see [27] and [35]). A small amount of noise has been added to these synaptic conductances of 0 mean and 0.0001 deviation. The conductances trace is shown in Figure 8a (black curve).

Numerical methods
Since we are considering neural models described both by ODE and SDE systems, we need to use dierent numerical methods to integrate them.
To integrate the system of ordinary dierential equations of the pyramidal cell model without noise, and also to integrate the EIF model used in the estimation procedure, we used the 4-th order Runge-Kutta method with a 0.01 ms time step; the choice was made to sample our data with an accuracy similar to that obtained experimentally.
On the other hand, the noisy pyramidal cell model was integrated using the Euler-Maruyama method with a xed time step of 0.01 ms.
For the interpolation of the estimated conductances, we use the Matlab function pchip, which is an implementation of the piecewise cubic Hermite interpolation. For our data sets, the errors induced by this method showed to be lower than those induced by the interpolation by cubic splines.
Finally, to numerically solve integrals, we use the Romberg method with positive integer n = 15 to ensure an error of order 10 −8 or lower (see [7]). When implicit equations need to be solved, we use the secant method with error tolerance 10 −5 .

Estimation procedures
In order to estimate the synaptic conductances from the membrane potential trace we use similar ideas to those presented in [17]. In that paper, which was uniquely devoted to a piecewise linear system (the McKean model), the knowledge of an analytical approximation of the period T as a function of the synaptic conductance g syn was essential. However, as exposed in the Introduction, in general this analytical approximation cannot be obtained and, therefore, we next propose a computational adaptation.
We rst t the in silico data, namely the membrane potential, to a base model; here, we use the exponential (EIF) model, already described in Subsection 2.1.1.
Once we have determined the parameters of the base model, we integrate it for dierent constant synaptic conductances g syn and compute the period of oscillation, T . This results in a table of corresponding T − g syn pairs, which can further be interpolated, using for instance the PCHIP interpolation, to obtain a continuous function G syn (T ).
Finally, given any membrane potential trace, we compute, for each interspike time interval (ISI), the ISI duration, T * , and its corresponding G syn (T * ). Repeating the procedure for all the ISIs, we obtain a discrete trace of synaptic conductances, which can nally be interpolated using also the PCHIP interpolation.
A schematic representation of this estimation procedure is presented in Algorithm 1 and displayed in Figure  1.
Algorithm 1 Estimation procedure 1. Data collection. Let V 0 := {V (t), t ∈ J} be the membrane potential of a neuron N on some time interval J in absence of synaptic inputs.
2. Fitting data to a base model. Fit V 0 to a base model (for instance, an integrate-and-re model). Call P the set of parameters of the base model that gives the best tting.
3. The f − g syn curve for the base model. Build up the function G syn (T ) for the base model with the set P. Depending on the base model, one possibility is to nd either G syn analytically, or its inverse T (g syn ) and then invert it. Otherwise, we will proceed computationally: (a) Consider a discretization of N equally-spaced synaptic conductances, that is, assuming that the base model is regularly spiking for g syn ∈ [g syn,0 , g syn,N ].
(b) By numerical integration, for each g syn,i with i = 1 . . . N , obtain the periodT i of the corresponding membrane potential. Note that this correspondence could also be obtained using one of the three other options presented in the Introduction (page 3). In all cases, a table containing the relation betweenT and g syn is obtained. Then, interpolate it to obtain a continuous function G syn (T ) using, for instance, a piecewise cubic Hermite interpolating polynomial (PCHIP).
4. Estimation: computing ISIs. Let V = {V (t), t ∈ J } be the membrane potential trace, for some time interval J , of neuron N subject to synaptic input bombardments. Assume that N emits n spk spikes along the time span J and let T j , for j = 1, . . . , n spk , be the j-th interspike interval of V.
6. Estimation: interpolation. Interpolate the set {(t j ,ĝ syn,j )} n spk j=1 , where t j is a point on the j − th ISI.
Here, we take t j the endpoint of the j − th ISI and use a piecewise cubic Hermite interpolating polynomial (PCHIP). * Note that for those models such that a qualitative expression of the T −gsyn curve can be provided, the interpolation step does not proceed.
Note that, in step 3b of Algorithm 1, several dierent synaptic conductances can provide the same period T . This lack of uniqueness is a trouble for the following steps in the algorithm. Even though in all cases we have treated, this issue had not happened, a control condition needs to be considered in this step to ensure the monotonic behaviour in the T − g syn curve. Hence, we have overcome this drawback by choosing the closest value to the last computed conductance.

Results
In order to test the estimation procedure, we apply Algorithm 1 (see Section 2.2) to membrane potential traces, often referred to as actual membrane potential traces, obtained computationally from the pyramidal cell model (see Section 2.1.1). This model will be considered either with additive noise (3), or without it (1). Moreover, in our simulations we consider the model with dierent synaptic drives to simulate spiking regimes without (Section 3.1) and with (Section 3.2) short subthreshold regimes interspersed (see Section 2.1.2 for more details in the synaptic drive).
For the dierent synaptic drives, we perform the estimation of conductances using dierent base models in Algorithm 1 step 3b. First, since we know that, by construction, the actual membrane potential follows the pyramidal neuron model, we consider this model as the base model, so we can test how powerful the estimation procedure is if we know the exact neural model. In this case, we deepen our results by exploring the cases where the data has been obtained with dierent amounts of noise and where the T − g syn curve has been constituted using dierent discretizations.
However, in most of the cases, the specic model that describes the neural data we have in consideration is not always known; hence, an approximation of it is required. For this reason, secondly, we use the exponential integrate-and-re (EIF) model with constant adaptation as a base model and try to approximate the activity of the pyramidal cells to it (as in [15]; see also Section 2.1).
Finally, we explore the situation of having a small sample of (g syn , T (g syn )) pairs which are obtained with some additive noise rather than having a reference model to obtain a more precise T − g syn curve. This specic case could be given when trying to obtain the T − g syn curve experimentally, where data is obtained with measurement noise and variability across trials. In this case, we rst generate a small sample of (g syn , T (g syn )) pairs where we add a small bias in order to simulate possible errors in the recordings. Then, we use a proper tting to obtain the G syn (T ) curve to be able to proceed with the estimation procedure (see Section 3.1.4 for more details on how). We will refer to this curve as the T − g syn experimental curve.
In Table 1   In this section we present the results obtained from simulations of the noisy pyramidal cell model using the prescribed three-frequency oscillatory drive dened in (8). This conguration allows us to obtain a membrane potential trace in the spiking regime without interspersed silent periods. The estimation has been performed using, as base model, the three mentioned cases: the pyramidal cell model (Section 3.1.1), the EIF model (Section 3.1.3) and the T − g syn experimental curve (Section 3.1.4).

Estimation when the neural model is known
In this case, we apply the estimation procedure given in Algorithm 1 considering the pyramidal neural model as a base model in step 3b of the algorithm. a b c Figure 3: Results of the estimation procedure applied to a membrane potential generated through the pyramidal cell model with no additive noise, and using this same model as a base model in the estimation procedure. Panel a shows the prescribed conductances (black curve) and the estimated conductances (red curve). Panel b shows the membrane potential (black curve) obtained from the pyramidal cell model described in Section 2.1.1 and using the prescribed conductances of panel a, and the reconstructed membrane potential (red curve) obtained also from the pyramidal cell model described in Section 2.1.1 but using the estimated conductances. Panel c shows the scatter plot of the prescribed conductances versus the estimated ones. The red line represents the identity line.
Results show that the estimated synaptic conductance trace follows the actual one before and after interpolating the estimated values (see Figure 3a). After a short transient discrepancy observed in the rst 20 ms, due to the the fact that we cannot have an a priori guess of the actual initial conditions, we observe a strong agreement between the estimation and the actual values, as we can see in Figure 3c where the actual conductance trace versus the estimated one (after being interpolated) lies in the vicinity of the identity line, meaning a good approximation. In fact, computing the errors between the two traces, we have that the mean squared error is of order O(10 −9 ) and the relative error is, in mean, of order O(10 −3 ), indicating that the Algorithm 1 accurately estimates the synaptic conductance time course (see Table 2). Note that, in this case, since we are using the same model to generate the data and to estimate the conductances, the only source of errors could arise from the process of inversion, which is also related to the discretization provided by the spiking frequency of the neuron. Nevertheless, the magnitude of the errors obtained corroborates that they do not penalize the estimation (see Figure 3b and the comment about the eects caused by the discretization in the next paragraphs). Moreover, the mean of the conductances after interpolating the data is also well estimated (relative error of 5.547 × 10 −4 ). Hence, when we reconstruct the membrane potential by using the pyramidal cell model with the estimated synaptic conductances, neither the spike frequency, the spike amplitude nor the spike onset time are modied (see Figure 3b), thus reinforcing the fact that the small errors caused in the estimation are not substantially important.
Two dierent sources of errors can arise during the estimation procedure. First, errors coming from the discretization of the conductances used to build up the T − g syn table that relates spiking periods with constant synaptic conductance values and, second, the measurement noise, which is embedded in the membrane potential obtained from the recordings. Next, we discuss how these errors are reected in the estimation.
Eects caused by the discretization We want to analyze the inuence on the estimations of the discretization of synaptic conductances in step 3b of Algorithm 1, that is, the eect of parameter δ g used to build up the T − g syn table. For this reason, and to avoid errors coming from additive noise, the actual membrane potential is also generated without considering noise. Then, we explore the errors (mean square error, MSE) for dierent discretization steps (see Table 2). Since the range of the actual conductances is approximately 8 × 10 −3 µS/cm 2 (see Figure 3a), we start using a step size of 10 −3 . Results show that the discretization has no inuence on the goodness of the procedure since increasing the discretization step does not yield a substantial increase of the errors. These results show that, when the base model is exactly the one that the neuron follows, the error of the estimation is given basically by the discretization of the conductances. Therefore, since the error coming from the discretization step size does not improve for values smaller than 10 −3 and, for this value the precision of the estimation is already acceptable, from now on, we will consider this discretization step size.      For values of the standard deviation up to σ V = 0.1 mV (approximately), the estimation provides errors similar to the deterministic case (one order of magnitude higher in the worst case, see Table 3). On the other hand, for greater magnitudes of noise (σ V ≈ 1 mV ), several discrete conductances are overestimated leading to higher errors both before and after the interpolation.

Estimation of conductances for the stellate cell model
In order to test the robustness of our method, the algorithm has been also applied to a membrane potential trace coming from the stellate cell model, see [33], with the synaptic input described in Section 2.1.2. The dynamics of the system presents oscillations both in the spiking and in the subthreshold regime. The base model considered in this case is the same model used to generate the membrane potential, that is, the stellate cell model. We also obtain good approximations (MSE= 8.37 × 10 −3 ), thus reinforcing the goodness of the estimation method. Figure 4 portrays the results of the estimation procedure. a b c Figure 4: Results of the estimation procedure for a membrane potential from a stellate cell model using the same model as a base model in the estimation procedure. Panel a shows the prescribed conductances (black curve) and the estimated conductances (red curve). Panel b shows the membrane potential (black curve), obtained from the stellate cell model described in (4) when the prescribed conductances of Panel a has been used, while the red curve depicts for the reconstructed membrane potential, obtained also from the stellate cell model but using the estimated conductances. Panel c shows the scatter plot of the prescribed conductances versus the estimated ones. The red line represents the identity line.
As we can see in Figure 4 a, the shape of the prescribed synaptic conductances is well approximated by that of the stellate cell model, even though some small uctuations are missing and the estimation is not as good as the case of the pyramidal neuron. We can also see the goodness by plotting the actual conductances versus the estimated conductances, showing values in the vicinity of the identity line. Moreover, we obtain that the mean squared error is of approximately 7.4 × 10 −4 , while the mean conductance has been estimated with a relative error of 2.047 × 10 −1 .

Estimation when the neural model is approximated by a base model
In many experimental cases, the exact model describing the target neuron is not known. One plausible way to apply our methodology is by tting the measured data to a reference model. In this section we explore this paradigm and study the errors of the estimation procedure. For this purpose, we generate membrane potential traces from the noisy pyramidal cell model (3) (see also Figure 2 in Section 2) with the synaptic drive shown in Figure 5a, and use them to t an exponential integrate-and-re model (5). The underlying assumption is that the tted EIF model provides a good approximation of the original model (see Figure 2).
In order to challenge our methodology, we consider the worst scenario regarding data discretization and noise levels according to our previous results (see Tables 2 and 3): we assume a discretization step size for the conductances (δ g = 10 −3 ) and a high standard deviation for the white noise in the membrane potential (σ V = 0.1 mV ). Even so, results still show a good performance of the estimation procedure, see Figure 5, unexpectedly similar to the results obtained in Section 3.1.1. a b c Figure 5: Results of the conductance estimation applied to a membrane potential generated from the noisy pyramidal cell model, but using the EIF model as a base model in the estimation procedure.
Panel a shows the prescribed conductances (black curve) and the estimated conductances (red curve). Panel b shows the membrane potential (black curve) obtained from the noisy pyramidal cell model (3) with σV = 0.1 and using the prescribed conductances of panel a, and the reconstructed membrane potential (red curve) obtained also from the noisy pyramidal cell model but using the estimated conductances as input. Panel c shows the scatter plot of the prescribed conductances versus the estimated ones. The red line represents the identity line.
As we can see in Figure 5a, some of the small oscillations are tracked properly, yielding to some temporal shifts on the onset of spikes when reconstructing the membrane potential (see Figure 5b). However, despite of this loss, the estimation keeps capturing the spiking frequency, and the conductances are estimated with a relative error of 2.7% (mean over time) whereas the MSE is 7.038 × 10 −7 . After interpolating the data, the relative error drops to 1.45% and the MSE is equal to 7.348 × 10 −7 . Note that the dierence between these errors and those provided by the same case when we use the exact neural model in the estimation procedure does not dier qualitatively (see the rst two rows of Table 4). Therefore, even though we assume that the exact model is not known, we still obtain good estimates of the time course of the synaptic conductance extracting the T − g syn relationship from an approximated base model (EIF in this case).
In the estimation procedure we have built a conductances' table to be able to make the estimation. However, Note that for this specic case, and when the conductances are supposed to be constant, there exists a closed expression for the EIF period, T EIF (g syn ), see equation (6). This fact allows us to replace the table of the estimation procedure for the solution T EIF (g syn ) = T * , where T * is the specic interspike considered in each iteration of the procedure. Knowledge of this expression allows us to estimate the conductances more eciently and with less numerical errors (see third row in Table 4).

Estimation without a reference base model
To fully explore all possibilities of how to extract a relation between the oscillatory period and a constant value of the conductance, apart from considering an approximate model from which we can build up a table or nd an explicit expression, we also consider the possibility of having the T − g syn relationship experimentally. This relationship could be directly obtained by using the dynamic clamp technique or indirectly through the I − V curve by using other techniques (such as current clamp). Indeed, if the current I only contains synaptic currents then the T − g syn curve can be obtained by computing g syn = I V −Vsyn and extracting the period T from the membrane potential trance.
Note that, if the T − g syn curve can be obtained with high accuracy (with a suciently ne discretization of points), then the estimation will present similar errors than those presented in Section 3.1 (see Table 2). However, if the discretization is poor, errors could be magnied.
In this section, we investigate the case where only a small sample of (g syn , T (g syn )) pairs are obtained. To computationally simulate this specic situation, we assume that the data set is generated from the noisy pyramidal cell model. By using this same model, we generate a small sample of ten (g syn , T (g syn )) pairs to construct the T − g syn curve. Moreover, to consider possible errors in the recordings (measurement noise and variability across trials), we have added noise to the resulting set of points. This bias has been included by adding a random normal distributed value with 0 mean and 0.5 standard deviation to the period component. Finally, to be able to estimate the conductance for any given period and since we are supposing that no extra information is known about the actual data, we have interpolated these data using the PCHIP method (see Figure 6a). Actual conductances Estimated conductances PCHIP approximation Figure 6: Results of the estimation procedure applied to a membrane potential generated through the pyramidal cell model with additive noise, and using a simulated T − gsyn curve in the estimation procedure. Panel a shows the actual T − gsyn curve (blue curve), the T − gsyn curve obtained with the EIF method (black curve), a simulation of an experimental T − gsyn curve (red dots), and an exponential tting of the experimental curve (red curve). Panel b shows the scatter plot of the prescribed conductances versus the estimated ones; the red line represents the identity line. Panel c shows the prescribed conductances (black curve) and the estimated conductances (red curve). In panel a, the experimental curve (red dots) has been computed using 10 dierent trials where noise has been added to the period (a random normal distributed value with 0 mean and 0.5 standard deviation) and the interpolated using the PCHIP method.
Applying now the estimation procedure to the membrane potential shown in Figure 3b (black trace) using the simulated experimental curve instead of an approximate base model to generate the T − g syn curve, the results show how the estimated conductances are clearly underestimated (see Figure 6b,c) and, even though the frequency of oscillation of conductances is maintained, the amplitude of the oscillation is reduced. This fact is mainly due to misestimations that arise in the T − g syn experimentally tted curve (with respect to that obtained with the actual model), specially the underestimation obtained for those values of conductances between 0.02 and 0.04, which is in fact the interval where the actual conductances lie.
However, in view of the magnitude of the conductances, these underestimations are not substantially significant, providing similar errors to the cases where the T − g syn curve is numerically computed from a base model (see fourth row in Table 4).
Note that the errors in the T − g syn curve can alter consistently the overall estimation causing an important increase of the errors in the estimation. One source of these misestimations is the interpolation method. We also tested an exponential tting, instead of the PCHIP interpolation, which gave us higher misestimations of the T − g syn curve in the working interval, see Figure 7. Results show, as previously, that the frequency of oscillation is maintained but the amplitude is worst underestimated providing an average relative error of approximately 15%. Hence, the errors are substantially magnied showing the importance of having a good approximation of the T − g syn curve to be able to properly estimate the conductances.    Figure 7: Results of the estimation procedure applied to a membrane potential generated through the pyramidal cell model with additive noise, and using a simulated T − gsyn curve in the estimation procedure. Panel a shows the actual T − gsyn curve (blue curve), the T − gsyn curve obtained with the EIF method (black curve), a simulation of an experimental T − gsyn curve (red dots), and an exponential tting of the experimental curve (red curve). Panel b shows the scatter plot of the prescribed conductances versus the estimated ones; the red line represents the identity line. Panel c shows the prescribed conductances (black curve) and the estimated conductances (red curve). In panel a, the experimental curve (red dots) has been computed using 10 dierent trials where noise has been added to the period (a random normal distributed value with 0 mean and 0.5 standard deviation and dots has been interpolated using the exponential tting).

Estimation in the spiking regime with interspersed short subthreshold periods
In this section we aim to explore whether the estimation procedure can be also used when short periods of subthreshold activity interrupt the spiking regime. Note that, if those silent periods are suciently large (for instance, greater than 100 ms), one could use a combination of the spiking estimation procedure developed in Section 2.2 and a subthreshold estimation procedure (for instance, the so-called QIF method presented in [36], which takes into account the spike initialization and possible subthreshold activity). However, for shorter subthreshold regimes, this combination could be not possible (since some methods require a certain amount of data to be used) or could cause a lot of misestimations induced by recent spike activity. For these reasons, we consider a synaptic conductance trace (see Section 2.1.2), depicted in Figure 8a (black trace), which causes spiking regimes interleaved with subthreshold regimes (of 50 ms of duration, approximately) when this is used into the pyramidal cell model (see Figure 8b (black trace)). a b c d Figure 8: Results of the estimation procedure applied to a non-tonic membrane potential from a noisy pyramidal cell model. Panel a shows the prescribed conductances (black curve) and the estimated conductances using the same model (blue trace), the EIF model (green trace) and the experimentally tted curve (red trace). Panels b, c, d show the reconstructed membrane potential when each of the mentioned estimated conductances in panel a has been respectively used as input to the base model (the noisy pyramidal cell model).
Results presented in Figure 8a show that the estimation procedure provided in Algorithm 1 is able to reproduce the synaptic drive, either if we use as a base model to generate the T − g syn table the pyramidal cell model (blue trace) or the EIF model (green curve). In both cases, the reconstruction of the membrane potential using the estimated conductances maintains the spike frequency, thus keeping unchanged the ring rate (Figure 8b and 8c, respectively). On the other hand, when the simulated experimental curve is used, conductances around [0.02, 0.04], are not well estimated, as it happened when considering just the spiking regime, see Section 3.1.4. These misestimates alter the spike frequency in the membrane potential reconstructed using the estimated conductances as synaptic inputs (see Figure 8b), by adding extra spikes in the silent regime and displacing the spiking time-onset in the spiking regime.
Errors of the estimation in the dierent cases are presented in Table 5, where we can see also that both the pyramidal cell model and its EIF approximation lead to similar errors, while the experimental curve leads to higher errors.

Discussion
Estimations of the time course of synaptic conductances arriving at a single cell (or a population) reveal relevant information to fully understand the functional connectivity of the brain. Several strategies have been developed for the estimation during subthreshold activity regimes [6,1,39,4,40,34,30,29,5,20,22,14,10,23,25,37,36,8,26,21]. However, the estimation during oscillatory (spiking) regimes has received less attention [8,17,11,24]. The main drawback is the predominance of ionic currents over synaptic ones in this regime. Nevertheless, as we already showed in a caricature model, see [17], if one gets the main information about the contribution of ionic currents and is able to track their activity during the spiking regime, then it is possible to infer the synaptic activity. In this paper, we extend this idea to more realistic spiking models and provide a systematic procedure that can be applied to a variety of neural data, coming both from computational simulations and from experiments. The idea is to exploit the input-response properties of the cell; more precisely, the relation between constant synaptic conductances and the spiking frequency (the f − g syn curve).
In the spiking regime, the most natural measurable quantity is the oscillation frequency, which constitutes a valuable piece of information for our purpose. The frequency-input relationship encompasses the average response of the neuron which includes both the direct eect of the input stimulus and the resulting activation of ionic currents. The fact that f − I curves are usually monotone implies the one-to-one correspondence between frequency and input and favors the resolution of the inverse problem. However, while in the case of a steady input (I constant) that induces a constant-frequency repetitive ring it is straightforward to infer the input current from the output frequency, the presence of the neuron's voltage in the synaptic current (I = g syn (t) (v(t) − V syn ) in their minimal expression) makes the estimation of the total input current highly sensitive to the uctuations in g syn . Here, we propose to consider an alternative frequency-input relationship, the f − g syn curve, which is more appropriate for the estimation problem since it implicitly incorporates the voltage dependency of the synaptic current and mitigates the eects of the uctuations. In a mathematical model, we can obtain the f − g syn or T − g syn curve, T = 1/f , numerically (in some cases, even analytically or quasi-analytically, see [17] and Section 2.1.1). For experimental studies, retrieving the T − g syn curve requires the use of dynamic clamp techniques on the isolated cell for xed g syn values within a physiologically plausible interval. In Algorithm 1 we summarize the steps to perform the estimation of synaptic conductances from the interspike intervals of the voltage traces and the knowledge of the T − g syn curve.
The second ingredient of the method presented in this paper is the choice of a reference/base model which must be both exible to t to a variety of neural data (choosing a suitable parameter set in each case) and simple enough to be able to compute the f − g syn curve in a precise way. We found the exponential integrate-and-re model (EIF) to be a convenient base model: on the one hand, its oscillation period can be well approximated by means of the semi-analytical formula (6) Our procedure has been applied rst to a representative computational model driven by a prescribed input conductance trace. We verify that the use of the f −g syn curve brings excellent predictions of the g syn uctuations along the prescribed trace, see Figure 3. We also observe better estimations when the synaptic conductances change on a time-scale longer than the mean interspike interval. We also analyze the eects of the discretization and the measurement noise, see Tables 2 and 3. In particular, we observe that levels of noise below 1 mV do not compromise the estimations in a substantial way since the MSE of the estimations is very low (MSE< 10 −6 ) in all cases. In Section 3.1.3 we explore the goodness of the estimations under the assumption that the EIF (base) model provides a good t to the neural data. By using a model for a noisy pyramidal cell to generate the data, we observe again excellent estimations of the synaptic conductance (MSE< 10 −6 ), see Figure 5. We conclude that we could reasonably estimate the synaptic conductances from other neural data to which the EIF could be tted.
We next study potential diculties that can be found in more realistic scenarios. First, in Section 3.1.4, we cope with the problem of estimating synaptic conductances when there is no base model that decently ts the data. We propose a procedure to come up with an approximation of the T − g syn curve and apply it to data generated by a computational model. Again, we achieve good estimations (MSE< 10 −6 ), see Table 4. For real cells, the dynamic-clamp technique allows to obtain the period T from a g syn value and, therefore, it would be interesting to design a protocol so that it was feasible to obtain the T − g syn curve. Since our method is mainly addressed to the spiking regime, the second potential diculty is the presence of interspersed short subthreshold periods. We show, see Figure 8, that these short subthreshold periods do not contamine the estimation of synaptic conductances and provide results with small errors.
In order to check the robustness of the estimation procedure for other types of neurons, we have also considered the stellate cell model (see Section 3.1.2), to which the quadratic integrate and re model (QIF, see [15]) has been tted. Results (not shown) reveal that the estimation procedure can also be used in cases where the tting is not as perfect as those obtained for the pyramidal cell model.
Our method is related in spirit to the ring-clamp method presented in [8], where the authors also consider a point neuron and, by calibrating its activity, they experimentally build up an approximated function to describe the changes in both the membrane potential threshold and the amplitude of the spiking regime (instead of the spiking frequency). Their model is able to discern between excitatory and inhibitory conductances; however, since their calibration is purely experimental, errors could arise from the measurement noise because of repetitive trials.

Conclusion
Estimation of synaptic conductances from the time course of a neuron's membrane potential is a challenging inverse mathematical problem, with strong implications in the study of the brain's functional activity. Few methods are known that address the problem when the neuron is spiking. The aim of our research is building up a general method to estimate synaptic conductances based on the knowledge of the frequency versus synaptic conductance curve (f −g syn ). This strategy is complemented with the choice of a mathematical model for which the f − g syn curve can be accurately computed and, simultaneously, can t a variety of data sources. We have shown that the exponential integrate-and-re model (EIF) fullls these requirements and can be a good reference model for future applications.
We have successfully applied the method to data generated from mathematical neuron models driven by dierent synaptic input proles. We have tested dierent levels of knowledge of the f − g syn curve, namely, directly obtained from the neuron model, using the f −g syn curve corresponding to the base model and assuming no model behind. Of course, the latter gives the worst results, but still within high accuracy. The method has also proven to be robust under normal levels of measurement noise and alternation of spiking and subthreshold activity.
These results allow to be optimistic about the application of the proposed methodology to obtain trustable estimations of the time course of synaptic conductance from real membrane potential traces presenting heterogeneous activity (combining subthreshold activity with intervals of spiking). Therefore, the next challenge is to discern between the excitatory (g exc ) and the inhibitory (g inh ) conductances, which is relevant to understand how brain activity results from the interaction of dierent types of neurotransmission. For this purpose, one needs additional (uncorrelated) measurable quantities to be able to build up a one-to-one relationship between pairs of measurable elements and the two parameters g exc and g inh . Following the ideas in [8], we think that the amplitude of the spikes is a potential candidate since it is easily measurable and typically highly uncorrelated with the frequency.