Inference on the effect of non homogeneous inputs in Ornstein-Uhlenbeck neuronal modeling

1.   Introduction and background

The first attempt to construct a stochastic model able to describe the spontaneous activity of a single neuron was due in 1964 to Gernstein and Mandelbrot in ^[1]. The authors assumed that the behavior of the neuronal membrane potential, subject to several simultaneous and independent acting potentials, is described by means of a stochastic diffusion process. They also proved that, by suitably choosing the involved parameters, the histograms of the experimentally recorded interspike intervals (ISI) can be approximated by the first passage time (FPT) probability density function (pdf) of a time-homogeneous Wiener process. Later, other models were proposed to include more characteristics of the behaviour of the membrane potential. Specifically, to describe the effect of inhibitory and excitatory inputs and to introduce the exponential decay of the membrane potential, models based on Ornstein-Uhlenbeck (OU) and Feller processes were considered (see, for example, ^[2] and reference therein).

Recently, studies also have been focused on the estimation of the parameters involved in neuronal models (see ^[3,4,5,6,7]) in order to understand information processing in neuronal mechanisms. An extensive review is also provided in ^[8]. In ^[3,4,6] estimators of the OU neuronal model input parameters are derived by looking at the firing regime of the process. Moreover, in ^[5], several estimation methods are compared in three different regimes: sub-threshold, threshold and supra-threshold. In the first regime the equilibrium point of the process is much smaller than the threshold and firing is caused only by random fluctuations and firings can be approximating by a Poissonian distribution. In the threshold regime the equilibrium point of the process is very closed to the firing threshold. In the supra-threshold regime the equilibrium point of the process is far above the firing threshold and the ISI are relatively regular (deterministic firing - which means that the neuron is active also in the absence of noise). In ^[5] it is shown that the moments method works best in moderate sub-threshold regime, the Laplace transform method only works in supra-threshold regime and the Integral equation method works in the entire parameter space. In ^[7] Bayesian estimators of the rates of the excitatory and inhibitory are provided to the aim to describe the stimulation effects on the synaptic input.

In ^[9,10,11] the OU process modeling the membrane potential is generalized assuming that in addition to the constant input there is a random component changing from two subsequent interspike intervals, generally caused by the naturally occurring variations of environment signaling or by other sources of noise not included in the model. In particular, in ^[9,10], the Gaussian assumption for the random effect is made and the the estimators of the parameters are explicitly obtained by MLE. In ^[11] non parametric estimators of the random effect are provided, by using a kernel strategy and the deconvolution method.

A wide literature on neuronal processing focuses on inputs arriving with generally non-constant frequencies and, in particular, on periodic inputs (see, for example, ^{[12,13,14,15,16,17]}). Generally, the periodic input signals are handled by means of a first passage time, viewed as the interval between two consecutive spikes.

In the present paper we consider a model for the activity of a single neuron based on a time non-homogeneous OU process and we propose a two-step procedure to infer the model. Specifically, we assume that, in the absence of input, the neuron membrane potential exponentially decays to a resting potential, with a time constant. Moreover, we assume that the neuron is stimulated by a sequence of inhibitory and excitatory postsynaptic potentials (PSP's) of constant amplitude occurring according to a Poisson's process with time-dependent rates. In this way we have two different process: one, referred as the control group, describes the membrane potential activity when the neuron is not stimulated, the other, denoted as a treated group, represents the neuronal membrane potential under the effect of the inhibitory and excitatory PSP's. Under the assumption of exponential decay, the control group is described by a time-homogeneous OU process, whereas the treated group is generally a non-homogeneous OU-type diffusion process.

The estimation procedure of the two processes consists of a first step in which the constant parameters of the control group are estimated, based on the Maximum Likelihood Estimation (MLE); in the second step of the procedure, by using the parameters estimated in the first step, the fit of the time-dependent functions involved in the treated group is performed by means of a generalized method of moments. Related to the second step, two methods are proposed: one based on the sample variance, while the second one is based on the sample covariance of the paths of the process.

We point out that the focus is on the time interval between two consecutive spikes in which the OU process is unlimited, in the sense that no absorbing or reflecting condition is considered. Nevertheless our procedures can be generalized in the case in which the control group as well as the treated group are described via a OU process restricted by suitable boundary conditions. This topic is beyond the aim of the present paper and it will be treated in future works.

The layout of the paper is the following. In Section 2 we describe the model and we derive some theoretical results for the resulting non-homogeneous diffusion process. Section 3 provides the two-step procedures to infer the neuronal model. In Section 4 some simulation experiments are performed to evaluate the effectiveness of the proposed procedures and to compare the performances of them.

2.   The model

We assume that, in the absence of stimuli, the behavior of the membrane potential of the neuron is described by a time-homogeneous OU diffusion process $\{X(t), \; t\geq t_0\}$ defined in $\mathbb{R}$ with drift and infinitesimal variance :

respectively. In Eq. (2.1), the parameter $\vartheta > 0$ is the time constant characterizing the exponential decay; $\mu$ is related to the resting potential $\varrho = \vartheta\, \mu$; finally, $\sigma^2 > 0$ represents the noise intensity. Hence, we assume that in the absence of inputs the membrane potential, described by the diffusion process $X(t)$, spontaneously decays to the resting potential $\varrho$ with a time constant $\vartheta$.

By stimulating the neuron with a sequence of constant inputs occurring with time-dependent rate, making use of standard procedure (see, for example, ^[23]), it can be proved that the evolution of the neuronal membrane potential is described via a generally time non-homogeneous OU diffusion process. Specifically, we assume that the neuron is reached by inhibitory and excitatory PSP's characterized by constant magnitude $\varepsilon$ occurring with rates:

where $a_i(t), \; a_e(t), \; u(t)$ are positive functions of time. Here, $\alpha_i(t)$ denotes the rate of the inhibitory inputs whereas $\alpha_e(t)$ represents the rate of the excitatory inputs. It can be proved that the dynamics of the neuronal membrane potential is described via a generally time non-homogeneous OU diffusion process $\{X_S(t), \; t\geq t_0\}$ defined in $\mathbb{R}$ whose infinitesimal moments are related to the rates given in (2.2). In particular, the drift and infinitesimal variance of $X_S(t)$ are

respectively, with

The sample-paths of $X_S(t)$ are solutions of the following stochastic differential equation:

where $t_0$ is the initial time; $x_0$ denotes the starting value and it is often assumed equal to the resting potential; $W(t)$ is a standard Brownian motion. The solution of (2.4) is

Clearly, the equation for the control group $X(t)$ can be obtained from (2.5) taking $m(t) = 0$ and $u(t) = \sigma^2$.

We note that (2.3) and (2.4) characterize a non-homogeneous OU process, known also as leaky-integrate-and-fire (LIF) model, describing the evolution of the membrane potential (see, for example, ^{[13,18,19,20]} and references therein). The theoretical study on the LIF model essentially concerns the analysis of the first passage time density viewed as the interspikes density, its asymptotic approximations (see, for example, ^[17,21]) and the analysis in the presence of a lower reflecting boundary (cf. ^[15,17,22]).

In the present paper we focus on the inference of model (2.3) and (2.4). Differently from ^[3,4,5,6], in which the focus is on the firing regime, in our assumptions, both $X(t)$ and $X_S(t)$ describe the neuronal membrane potential between two consecutive spikes, so no boundary conditions are assumed. This case was also considered in ^[9,10,11], in which the membrane potential, during a generic ISI, is a stationary OU process with a random equilibrium point having an assigned distribution; so the effect of the neuronal input remains constant within each ISI and it changes randomly from one ISI to another. Further, the input is able to influence only the drift of the process.

Moreover, we assume that the membrane potential, during the generic ISI, is a generally non-stationary process. Moreover, in our model, the neuronal inputs can influence both the drift and the infinitesimal variance by means of the two time-dependent functions $m(t)$ and $u(t)$.

2.1.   Probabilistic properties

The transition pdf $f_S(x, t|x_0, t_0)$ of $X_S(t)$ is solution of the Kolmogorov equation

and of Fokker-Plank equation

Moreover, the transition pdf satisfies the initial delta condition

By using the transformation:

the equations (2.6) and (2.7) can be reduced to the analogous equations of a standard Wiener process for which the end points of the diffusion interval, $\pm\infty$, are natural boundaries. So we obtain that $f_S(x, t|x_0, t_0)$ is a normal pdf:

with conditional mean and variance given by

respectively.

The transition pdf $f(x, t|x_0, t_0)$ of $X(t)$ can be obtained from (2.8) setting $m(t) = 0$ and $u(t) = \sigma^2$ so, for $t\geq t_0$, it is normal with conditional mean and variance

We note that

Further, the process $X(t)$ admits an asymptotic behaviour with steady-state density:

Under suitable hypothesis on the functions $m(t)$ and $u(t)$, also $X_S(t)$ admits an asymptotic behaviour. In particular, if the following limits:

exist finite, the steady state pdf $w_S(x)$ of $X_S(t)$ is normal with mean $\mu\, \vartheta+L_1$ and variance $L_2$.

Finally, we remark that $X_S(t)$ is a Gauss-Markov process with mean and covariance functions:

with $t > 0$ and $0 < \tau < t$, respectively. From (2.13) and (2.14) it is easy to obtain the mean and covariance for the control group $X(t)$:

3.   Fitting the model

In this section we consider a two-step procedure to infer the model. Firstly, we use the control group, modeled by the process $X(t)$ given in (2.1), to estimate the constant parameters $\vartheta, \mu$ and $\sigma^2$ via the MLE, and secondly, we use a treated groups described by $X_S(t)$ defined in (2.3) to fit the unknown functions $m(t)$ and $u(t)$. We suppose that the control group and the treated group are observed at the same discrete time instants. In the neuronal context this means to record the data of two groups of identical neurons, the first one without time-dependent stimuli and the second one with generally time-dependent inputs.

3.1.   Estimation of the constant parameters

In the following we consider a discrete sampling of the process $X(t)$ in (2.1) based on $d_1$ sample paths at the times $t_{i\, j}$ for $i = 1, \ldots, d_1$ and $j = 1, \ldots, n$ with $t_{i\, 1} = t_1$ for $i = 1, \ldots, d_1$. Let $x_{i\, j}$ be the observed values at times $t_{i\, j}$. We assume that the observations $x_{i\, j}$ are equally spaced at $\Delta = t_{i\, j}-t_{i\, j-1}\ \forall i, j$. Further, we assume $x_{i\, 1} = x_0$ for $i = 1, 2, \ldots, d_1$.

The likelihood function for the parameters $(\vartheta, \mu, \sigma^2)$ is

The MLE can be explicitly obtained (see, for example, ^[3,4]):

where

Since the process $X_S(t)$ is characterized by the same constant parameters with respect to the control group $X(t)$, the estimates $\widehat \vartheta, \widehat\mu, \widehat\sigma^2$ obtained by (3.1) will be used in the second step of the procedure to fit the unknown functions $m(t)$ and $u(t)$ in (2.4).

3.2.   Fitting the functions $m(t)$ and $u(t)$

In order to fit the functions $m(t)$ and $u(t)$ in (2.4) in Propositions 1, 2 and 3 we provide equations related to suitable characteristics of the process $X_S(t)$.

Proposition 1. The function $m(t)$ in (2.4) is given by

with

where $M_S(t|x_0, t_0)$ and $M(t|x_0, t_0)$ are the conditional means of the processes $X_S(t)$ and $X(t)$, respectively.

Proof. From (2.12), we obtain:

with $h(t)$ given in (3.3). From (3.4), deriving with respect to $t$, we derive (3.2).

Proposition 2. The function $u(t)$ in (2.4) is given by

where $V_S(t|t_0)$ is the conditional variance of the treated group $X_S(t)$.

Proof. Deriving both sides of (2.9) with respect to $t$, we obtain:

from which one has Eq. (3.5). Alternatively, the unknown function $u(t)$ in (2.4) can be obtained by looking at the covariance function $c_S(\cdot, \cdot)$ of the process $X_S(t)$ as shown in the following proposition.

Proposition 3. The function $u(t)$ in (2.4) is given by

Proof. Eq. (3.6) is obtained from (2.14) by deriving both sides with respect to $\tau$.

Proposition 1 can be used to fit the function $m(t)$ in (2.4), Propositions 2 and 3 are able to provide two different estimates of the function $u(t)$. In the following we specify the procedures.

3.3.   Procedure 1

Let us consider a discrete sampling of the process $X_S(t)$ in (2.4) based on $d_2$ sample paths at the times $t_{i\, j}$ for $i = 1, \ldots, d_2$ and $j = 1, \ldots, n$ with $t_{i\, 1} = t_1$ for $i = 1, \ldots, d_2$. Let $x^{(S)}_{i\, j}$ be the observed values at times $t_{i\, j}$, $i = 1, \ldots, d_2$ and $j = 1, \ldots, n$. We assume, as for the control group, that the step between two consecutive observations is $\Delta$ and $x^{(S)}_{i\, 1} = x_0$ for $i = 1, 2, \ldots, d_2$. Note that, in order to apply Eq. (3.2), data from the control and treated groups have to be observed at the same time instants. Considering Proposition 1 and Proposition 2, the following procedure can be formulated:

$1.$ From the data of the control group, estimate the parameters of process $X(t)$, obtaining the ML estimate $\widehat\vartheta, \; \widehat\mu, \; \widehat\sigma^2$.

$2.$ Let

be the sample means of $X(t)$ and $X_S(t)$, respectively, for $j = 1, 2, \ldots, n$.

$3.$ Obtain the points:

$4.$ Interpolate the points $h_j$ in (3.7) to obtain an estimate, $\widehat h(t)$, of the function $h(t)$ in (3.3).

$5.$ Consider

as approximation of the function $m(t)$.

$6.$ Obtain the points

$7.$ Interpolate the points $v_j$ in (3.8) to obtain an estimate, $\widehat v(t)$, of the function $V_S(t|t_0)$ in (2.9).

$8.$ A fitted function of $u(t)$ in (2.4) is:

3.4.   Procedure 2

An alternative procedure can be obtained by looking at Eq. (3.6) instead of Eq. (3.5). Precisely, by considering a discrete sampling of the process $X_S(t)$ as in Procedure 1, Procedure 2 works as before till point 5., whereas the function $u(t)$ is fitted by using Proposition 3. Therefore, the next steps change in

6$^\prime$. Obtain the points

7$^\prime$. Interpolate the points $c_j$ in (3.9) to obtain an estimate, $\widehat c(t)$, of the function $c_S(t-1, t)$ in (2.14).

8$^\prime$. A fitted function of $u(t)$ in (2.4) is:

We point out that the two suggested procedures differ between them only in the fit of the infinitesimal variance $u(t)$.

4.   Simulation results

In order to evaluate the proposed procedures we present some experiments in which $50$ sample-paths including 500 observations of the processes $X(t)$ and $X_S(t)$ with $t_0 = 0$, $\Delta = 0.1$, $x_0 = -70$ (representing the resting potential) are simulated.

We point out that the simulation experiments can be performed starting from the classical Euler's discretization of the corresponding stochastic differential equation. Since the process $X_S(t)$ is Gauss-Markov, an alternative way to simulate sample-paths of the process $X_S(t)$ in (2.4) is to consider the following discretization (see, for example, ^[24]):

with $t_k = t_0+k\Delta, \; k = 1, \ldots 499$. Here, $\xi_1, \xi_2, \ldots$ is a sequence of independent and identically distributed random variables with standard Gaussian distribution. Clearly, setting $m(t) = 0$ and $u(t) = \sigma^2$ for all $t$ in (4.1), sample-paths of $X(t)$ can be obtained.

We simulated the sample paths of the processes $X(t)$ and $X_S(t)$ by using both Eq. (2.4) and Eq. (4.1), the obtained estimates in two cases do not seem to show relevant differences, so in the following our simulation experiments are based on the simplest Euler's discretization.

We have chosen several functions $m(t)$ and $u(t)$, obtaining 4 non-homogeneous OU processes $X_S(t)$. Specifically, we have considered periodic functions in the drift and in the infinitesimal variance of the process, so to model periodic neuronal input widely admitted in the literature. Further, the infinitesimal variance of Case 3 was also chosen in ^[17]. Finally, some test cases are also considered to validate our procedures. In all the cases, for the process modeling the control group we have chosen the following values for the parameters: $\vartheta = 1, \; \mu = -70, \; \sigma = 0.05$, so the original homogeneous model has infinitesimal drift and variance

For each process we apply the procedures proposed in Sections 3.3 - 3.4 and we compare the true functions with their fitted versions. Moreover, to find out whether there are significant differences between the two fitting methods for the function $u(t)$, based on the sample variance and on the sample covariance, we compare the theoretical conditional mean and variance of the process $X_S(t)$ with the corresponding fitted functions. Finally, the simulation study includes the calculation of the Mean Absolute Error (MAE) for the conditional mean and the variance of $X_S(t)$.

First of all, by focusing on the control group, we remark that the use of MLE as described in Section 3.1 has provided the following results:

4.1.   Simulation results for the process $X_S^{(1)}(t)$

We consider the process $X_S^{(1)}(t)$ characterized by the following infinitesimal moments:

where $m_1(t) = c\, [b+\sin(\omega\, t)]$ and $u_1(t) = \sigma^2$. In this case we assume that the inputs rates influences only the drift of the treated process. For this process from (2.9), (2.12), (2.13) and (2.14) we have:

where $M(t|y, \tau)$ and $V(t|\tau)$, defined in (2.10) and (2.11), are the conditional mean and variance of the control group, whereas $c(\tau, t)$ is the covariance function of $X(t)$ defined in (2.15).

We assume $b = 0, \; c = 0.1$ and $\omega = 1$. In Figure 1 the function $\widehat m_1(t)$ (green line) along with the true function $m_1(t)$ (red line) are shown on the top. Smoothed version of $\widehat m_1(t)$ is also displayed. On the bottom the results obtained for $u_1(t)$ are plotted: on the left we show the results obtained making use of Procedure 1, whereas on the right the results obtained by Procedure 2 are plotted. We note that both the procedures give satisfactory results. Moreover, the approximation of $u_1(t)$ based on Procedure 1 shows an higher variability with respect to that one obtained by Procedure 2. Further, the approximation based on Procedure 2 seems to systematically underestimate the true function $u_1(t)$. To further investigate the comparison of the two suggested procedures, in Figure 2 the conditional mean (on the left) and variance (on the right) of the process (4.2) are plotted together with their fitted versions. We can see that he fitted conditional mean is able to capture the period structure of the conditional mean and also the amplitude of its periodicity. Further, the conditional variance of $X_S^{(1)}(t)$ in (4.2) seems to be well fitted by using Procedure 2, while Procedure 1 seems to overestimate the true conditional variance.

Finally, we simulate $N = 50$ replications of model (4.2) for several choices of the parameter $\sigma$ and of the observation step $\Delta$. For each Monte Carlo replication, indexed by $r = 1, 2, \ldots, N$, the absolute errors for both the procedures with respect to the mean and the variance of the process $X_S(t)$ have been calculated. In particular, denoting by $(\widehat M_S(t|x_0, t_0))_r$ the fitted conditional mean at the replication $r$ and by $(V^{(v)}_S(t|t_0))_r$ and $(V^{(c)}_S(t|t_0))_r$ the fitted conditional variance obtained by the Procedure 1 and 2, respectively, we consider the following absolute errors

Then, we derive the mean absolute errors as follows

We point out that the MAEs are dependent on time. In order to give a measure of the errors, we consider the mean over the time of the MAEs for each estimated functions, for example for the function $M(t)$

The obtained values of $mean(MAE_M(t)), mean(MAE_V^{(v)}(t))$ and $mean(MAE_V^{(c)}(t))$ are listed in Table 1 for several choices of the parameter $\sigma$ and of the observation step $\Delta$. Clearly, in almost all the cases the errors increases as $\sigma$ and the step $\Delta$ increase. This is more evident in the second column in which the errors on the conditional variance of $X_S^{(1)}(t)$ by Procedure 1 are shown. Further, comparing the errors on the conditional variance by the two procedures we can see that Procedure 2 provides always smaller mean errors, suggesting that it is to prefer when one is interested to the conditional mean and variance.

4.2.   Simulation results for the process $X_S^{(2)}(t)$

Let $X_S^{(2)}(t)$ be the stochastic diffusion process having infinitesimal moments:

with $m_2(t) = c\, t+d$ and $u_2(t) = \sigma^2$. Also in this case the inputs rates influences only the drift of the treated process. For the process $X^{(2)}(t)$ from (2.9), (2.12), (2.13) and (2.14) we have:

where $M(t|y, \tau)$ and $V(t|\tau)$, defined in (2.10) and (2.11), are the conditional mean and variance of the control group, whereas $c(\tau, t)$ is the covariance function of $X(t)$ defined in (2.15).

We assume $c = 0.01$ and $d = 0$. In Figure 3 and in Figure 4 the same analysis as in Figure 1 and 2 is performed. Also in this case the fitted functions seem to be very close to the true ones although Procedure 2 underestimates the true infinitesimal variance $u_2(t)$. Further, Figure 4 shows that the fitted conditional mean well recognizes the linear trend of the conditional mean $M_S^{(2)}(t|x_0, t_0)$ and only slightly overestimate its slope. By looking at the conditional variance $V_S^{(2)}(t|\tau)$, we observe similar results with respect the case in Section 4.1.

In Table 2 the errors on the conditional variance by the two procedures for the process $X_S^{(2)}(t)$ are listed for the same choices of Section 4.1, and they confirm the results of Table 1.

4.3.   Simulation results for the process $X_S^{(3)}(t)$

We consider the process $X_S^{(3)}(t)$ with infinitesimal moments:

where $m_3(t) = c\, (b+\sin(\omega\, t))$ and $u_3(t) = \sigma^2\, (1-e^{-2\, t/\vartheta})^2$. In this case the inputs rates influences both the drift and the infinitesimal variance of the treated process. For the process $X^{(3)}(t)$ from (2.9), (2.12), (2.13) and (2.14) we have:

with $M^{(1)}(t|y, \tau)$ given in (2.10). We set $c = 0.1$, $b = 1.2$ and $\omega = 1$. In Figure 5, on the top, the function $m_3(t)$ along with its approximation and its smoothed version are plotted. On the bottom the results obtained for $u_3(t)$ for the two procedures are shown. Also in this case, the approximation of $u_3(t)$ based on the sample variance of the process shows an higher variability with respect to that one obtained by Procedure 2 and this latter seems to underestimate the function $u_3(t)$. Figure 6, showing the conditional mean and variance, this latter approximated by the two procedures, shows the behaviour observed in Section 4.1 and 4.2. In Table 3 the mean errors of the conditional mean and variance are listed for various choices of $\sigma^2$ and of the step $\Delta$.

4.4.   Simulation results for the process $X_S^{(4)}(t)$

For the last simulation experiment we consider the process $X_S^{(4)}(t)$ with:

where $m_4(t) = 0$ and $u_4(t) = \sigma^2\, c\, (b+\sin(t))$.

In this case the inputs rates influences only the infinitesimal variance of the treated process. For the process $X^{(4)}(t)$ from (2.9), (2.12), (2.13) and (2.14) we have:

where $M(t|y, \tau)$ denotes the conditional mean of the control group modeled via $X(t)$. The numerical study is performed for $c = 0.1$, $b = 1.2$ and $\omega = 1$. The fitted functions along with the true one are plotted in Figure 7. As shown in Figure 8, both the suggested procedures are able to detect the periodicity in the infinitesimal variance, although Procedure 2 provides a fitted conditional variance slightly delayed with respect the true one. In Table 4 the related mean errors are listed for various choices of $\sigma^2$ and of the step $\Delta$.

5.   Conclusions

In this paper we have considered a time non-homogeneous OU process. It can be viewed as a model for the membrane potential activity of a single neuron, stimulated by a sequence of inhibitory and excitatory post-synaptic potentials with generally time-dependent rates. We have proposed a two-step procedure to infer the model. It uses a control group modelled by means of a homogeneous OU process, and a treated group. In the first step the constant parameters of the model are estimated making use of a MLE applied to the control group. In the second step we are able to fit the non-homogeneous terms of the process. This second step uses the estimates obtained in the first step and, by a generalized method of moments, it provides the fit of the time-dependent functions in the process describing the treated group. To estimate the infinitesimal variance of the non-homogeneous OU process, we propose two procedures: the first one is based on the conditional variance and the second one is based on the covariance function of the treated group. Several numerical simulations show that the first procedure seems to be able to well capture the time-dependent infinitesimal variance, whereas the second one better fits the conditional variance of the process.

Acknowledgements

This research is partially supported by MIUR – PRIN 2017, Project "Stochastic Models for Complex Systems" and by INDAM-GNCS. The authors acknowledge the constructive criticism of anonymous reviewers on an earlier version of this paper.

Conflict of interest

The authors declares no conflict of interest.