Adaptive stimulus design for dynamic recurrent neural network models

The continuous time recurrent neural networks have a similar structure to that of the discrete time counterparts that are often met in artificial intelligence studies. In Figure 1, one can see a general continuous time network that may have any number of neurons.

where $\beta_{e}$ and $\beta_{i}$ are the reciprocals of the time constants $\tau_{e}$ and $\tau_{i}$. They are taken to the right for easier manipulations of the equations. Note that this equation is written in matrix form to be conformed to the formal non-linear system forms. A descriptive illustration related to (5) is presented in Figure 1b. It should also be noted that, in (4) and (5) the weights are all assumed as positive coefficients and they have signs in the equation. So negative signs indicate that originating neuron is inhibitory (tend to hyper-polarize the other neurons in the network).

The theoretical response of the network in (4) will be the firing rate of the excitatory neuron as ${r_{e}=g_{e}\left(V_{e}\right)}$. In the actual environment, the neural spiking due to the firing rate ${r_{e}\left(t\right)}$ is available instead. While introducing this research, it is stated that this spiking events conform to an inhomogeneous Poisson process which is defined below:

where

is the mean number of spikes based on the firing rate $r_{e}(t)$ which varies with time, and $N(\tau)$ indicates the cumulative total number of spikes up to time $\tau$, so that $N\left(t+\Delta t\right)-N\left(t\right)$ is the number of spikes within the time interval ${\left[t,t+\Delta t\right)}$. In other words, the probability of having ${k}$ number of spikes in the interval ${\left(t,t+\Delta t\right)}$ is given by the Poisson distribution above.

Consider a spike train $(t_{1},t_{2},\ldots,t_{K})$ in the time interval $(0,T)$ (here $0\leq t_{1}\leq t_{2}\leq\ldots\leq t_{K}\leq T$ so $t$ and $\Delta t$ become $t=0$ and $\Delta t=T$). Here the spike train is described by a list of the time stamps for the $K$ spikes. The probability density function for a given spiking train $(t_{1},t_{2},\ldots,t_{K})$ can be derived from the inhomogeneous Poisson process [26, 27]. The result reads:

This probability density describes how likely a particular spike train $(t_{1},t_{2},\ldots,t_{K})$ is generated by the inhomogeneous Poisson process with the rate function $r_{e}\left(t,{\bf x},\theta\right)$. Of course, this rate function depends implicitly on the network parameters and the stimulus used.

The network parameters to be estimated are listed below as a vector:

which includes the time constants and all the connection weights in the E-I network. Our maximum-likelihood estimation of the network parameters is based on the likelihood function given by (8), which takes the individual spike timings into account. It is well known from estimation theory is that maximum likelihood estimation is asymptotically efficient, i.e., reaching the Cramér-Rao bound in the limit of large data size.

To extend the likelihood function in (8) to the situation where there are multiple spike trains elicited by multiple stimuli, consider a sequence of $M$ stimuli. Suppose the $m$-th stimulus ($m=1,\ldots,M$) elicits a spike trains with a total of $K_{m}$ spikes in the time window $[0,T]$, and the spike timings are given by ${S}_{m}=\left(t_{1}^{(m)},t_{2}^{(m)},\ldots,t_{K_{m}}^{(m)}\right)$. By (8), the likelihood function for the spike train $S_{m}$ is

where $r_{e}^{(m)}$ is the firing rate in response to the $m$-th stimulus. Note that the rate function $r_{e}^{(m)}$ depends implicitly on the network parameters $\theta$ and on the stimulus parameters. The left-hand side of (10) emphasizes the dependence on network parameters $\theta$, which is convenient for parameter estimation. The dependence on the stimulus parameters will be discussed in the next section.

We assume that the responses to different stimuli are independent, which is a reasonable assumption when the inter-stimulus intervals are sufficiently large. Under this assumption, the overall likelihood function for the collection of all $M$ spike trains can be written as

By taking natural logarithm, we obtain the log likelihood function:

Maximum-likelihood estimation of the parameter set is given formally by

Numerical issues related to this optimization problem will be discussed in Sections 2.5 and 2.6. In addition, some discussion on the local maxima problems is provided in Section 3.3.

The optimal design method generates the stimuli by maximizing a utility function, or an objective function. The basic idea is that these stimuli are designed so as to elicit responses that are most informative about the network parameters. In optimal design method, the utility function $U({\bf x},\theta)$ depends on the stimulus parameters $\bf x$, but typically also on the model parameters $\theta$. An intuitive explanation of the dependence on the model parameter is best illustrated with an example. Suppose we want to estimate a Gaussian tuning curve model with unknown parameters although we may have some idea about the sensible ranges of these parameters. To estimate the height of the tuning curve accurately, we should place a probing stimulus around the likely location of the peak. To estimate the width, the probing stimulus should go to where the tuning curve is likely to have the steepest slope. For the baseline, we should go for the lowest response. This simple example illustrates two facts: first, optimal design depends on our knowledge of possible parameter values; second, the elicited responses in an optimally design experiment are expected to vary over a wide dynamic range as different parameters are estimated.

Once the utility function $U({\bf x},\theta)$ is chosen, the optimally designed stimulus may be written formally as:

where the network parameters $\theta$ can be obtained by maximum-likelihood estimation from the existing spike data as described in the preceding section. Here the stimulus is specified by vector $\bf x$, which is a set of parameters rather than the actual stimulus itself. Direct computation of the actual time-varying stimulus is not easy because no closed analytical form of the objective function is available and furthermore the computation of the optimal control input generally requires a backward integration or recursion. Instead of struggling with this difficulty, one can restrict the stimulus $I$ to a well known natural form such as sum of phased cosines as shown below:

where $A_{n}$ is the amplitude, $\omega_{n}$ is the frequency of the $n$-th Fourier component, and $\phi_{n}$ is the phase of the component. We choose a base frequency ${\omega_{1}}$ and set the frequencies of all other components as the harmonics: $\omega_{n}=n\omega_{1}$ for $n=1,\ldots,N$. Now the stimulus parameters can be summarized by the stimulus parameter vector:

We sometimes refer to $\bf x$ as the stimulus, with it understood that it really means a set of parameters that uniquely specify the actual stimulus $I$.

Some popular choices of the objective function are based on the Fisher information matrix, which is generally defined as:

where $p({\bf r}\mid{\bf x},\theta)$ is the probability distribution of the response $\bf r$ to a given stimulus $\bf x$. One should here note that, the term $p({\bf r}\mid{\bf x},\theta)$ is different from the terms $p\left(t_{1},t_{2},\ldots,t_{K}\right)$ in (8) and $p\left(S_{m}\mid\theta\right)$ in (10). The latter two are derived from the spiking likelihood function presented in [26, 27] whereas the definition $p({\bf r}\mid{\bf x},\theta)$ is derived from (6). However only the related results are provided in this text in order to save space. In (17), $\left<\;\right>$ represents the weighted average over all possible responses $\bf r$ to all available sets of stimuli $\bf x$. This average is calculated based on the response probability distribution $p({\bf r}\mid{\bf x},\theta)$. Stimulus variable $\bf x$ can be represented by a group of parameters such as the Fourier series parameters in (15) which are $A_{n},\ \omega_{n}$ and $\phi_{n}$. The frequency parameter may or may not be fixed depending on the requirements or computational burden.

The Fisher information matrix reflects the amount of information contained in the noisy response ${\bf r}$ about the model parameters ${\theta}$, assuming a generative model given by the conditional probability $p({\bf r}\mid{\bf x},\theta)$. So the stimulus designed by maximizing a certain measure of the Fisher information matrix (17) is expected to decrease the error of the estimation of the parameters ${\theta}$.

The utility function $U$ can be chosen as a scalar function of the Fisher information matrix $\bf F$. One theoretically well founded popular choice is the D-optimal design:

although the determinant of the Fisher information matrix is not always easy to optimize. The A-optimal design is based on the trace of the Fisher information matrix and is easier to optimize:

Another alternative is the E-optimal design where the objective function is the smallest eigenvalue of the Fisher information matrix. In this paper the A-optimality measure of the information matrix is preferred. There is an obvious reason for this preference. As the computational complexity of the optimization algorithms are expected to be high, the necessity of numerical derivative computation should be avoided as much as possible. Since it is not easy to evaluate the derivatives of the eigenvalues and determinants by any means other than numerical approximations it will be convenient to apply a criterion like A-optimality which has a direct relationship like the sums of the diagonal elements.

In the beginning of this section we mentioned that an optimally design stimulus is expected to depend on which parameter is supposed to be estimated. Since a scalar utility function in (18) or (19) depends on all the parameters $\theta=[\theta_{1},\theta_{2},\ldots]$, optimizing a single scalar function is sufficient to recover all the parameters. When a sequence of stimuli are generated by optimal design, the stimuli may sometimes alternate spontaneously as if the optimization is performed with respect to each of the parameters one by one [12].

In our network model, the recorded spike train has an inhomogeneous Poisson distribution with the rate function ${r_{e}\left(t\right)}$. We write this rate as $r_{e}\left(t,{\bf x},\theta\right)$ to emphasize that it is a time-varying function that depends on both the stimulus $\bf x$ and the network parameters $\bf x$. For a small time window of duration $\Delta t$ and centered at time $t$, the Fisher information matrix entry in (17) is reduced to:

Since the Poisson rate function $r_{e}$ varies with time, the A-optimal utility function in (19) should be modified by including integration over time:

Here the time window $\Delta t$ is ignored because it is a constant coefficient that does not affect the result of the optimization.

For convenience, we can also define the objective function with respect to a single parameter $\theta_{k}$ as follows:

The objective function in (21) is identical to $\sum_{k=1}^{8}U_{k}$.

The optimization of the D-optimal criterion in (18) is not affected by parameter rescaling, or changing the units of parameters. For example, changing the unit of parameter 1 (say, from msec^-1 to sec^-1) is equivalent to rescaling the parameter by a constant coefficient: $\theta_{1}\to c\theta_{1}$. The effect of this transformation is equivalent to a rescaling of the determinant of the Fisher information matrix by a constant, namely, $\det{\bf F}\to(\det{\bf F})/c^{16}$, which does not affect the location of the maximum of (18). By contrast, the criterion function in (19) or (21) are affected by parameter rescaling. A parameter with a smaller unit would tend to have larger derivative value and therefore contribute more to (21) than a parameter with a large unit. To alleviate this problem, we use $U_{k}$ one by one to generate the stimuli. That is, stimulus 1 is generated by maximizing $U_{1}$, and stimulus 2 is generated by maximizing $U_{2}$, and so on. Once the 8th stimulus is generated by maximizing $U_{8}$, we go back and use $U_{1}$ to generate the next stimulus, and so on. Finally, an alternative way to get rid of scale dependence is to introduce logarithm and use $U=\sum_{k}\ln U_{k}$ as the criterion, which, however, may become degenerate when $U_{k}$ approaches 0.

As just explained in the previous section about the computational issues in this research, the gradient computation decreases the computation durations considerably. The main issue with this fact is the lack of closed form expressions like in the case of static non-linear mappings as the model. In researches such as [33], [2] and [34] the gradients are computed as a self contained differential equation which is formed by taking the derivatives of the model equations (3) and (4) from both sides.

Compiling all the information in this section one can write the gradient of the Fisher Information Measure (i.e. the Fisher Information Matrix with a certain optimality criterion such as A-Optimality). In the beginning of this section, it is stated that the sensitivity levels of the firing rate w.r.to different network parameters are different and thus it would be convenient to maximize the fisher information for a single parameter at a time.

The maximization of the Fisher Information Measure in trace form requires its gradients and they involve second order cross derivatives of the membrane potentials with respect to parameters ${\left(\theta_{k}\right)}$ and stimulus parameters ${\left(I_{l}\right)}$.

To compute the gradients needed for optimizing the objective function based on Fisher Information or for performing maximum-likelihood estimation, we need to evaluate the derivatives of the mean firing rate ${r_{e}=g_{e}\left(V_{e}\right)}$ with respect to the network parameters $\theta_{k}$ in (9) or the stimulus parameters $x_{j}$ in (23). The first and the second order derivatives are:

These formulas are expressed in terms of the derivatives $\frac{\partial{V_{e}}}{\partial{\theta_{k}}}$, $\frac{\partial{V_{e}}}{\partial{x_{l}}}$ and $\frac{\partial^{2}V_{e}}{\partial x_{l}\partial\theta_{k}}$, which are regarded as dynamical variables that can be solved from the three differential equations (27)–(29). Here the initial conditions were always assumed to be the equilibrium state. The initial values of the derivatives can be set to zero as recommended in [2].

Lastly, for maximum likelihood estimation, one need the gradient of the log likelihood function of spike trains in (12):

which is written in terms of derivatives that can be evaluated by (31).

Up to this point, the theoretical grounds of this research are presented. In order to achieve the results one needs two separate maximization algorithms targeting (13) and (14). Optimization algorithms have varieties in certain aspects. One of these classifications is the gradient requirements. The global optimization algorithms such as genetic algorithms (MATLAB^® ga), simulated annealing (MATLAB^® simulannealbnd) and pattern search (MATLAB^® patternsearch) do not require gradient computations. However due to their stochastic nature, they do a lot of computations which increase the computation duration drastically. In addition, that stochastic nature leads to different results at the end of different runs. Disabling this randomness, might have adverse effects on the optimization performance and their ability of finding a global optimum solution. In order not to be trapped in these issues we will stick to gradient based optimization routines such as constrained interior-point gradient descent. MATLAB^® provides these algorithms through its optimization toolbox. Interior point and similar methods are available in fmincon function. We can write the following facts about the optimization algorithm:

• 1.

  MATLAB^®’s fmincon allows the user to set constraints on the solution.In the optimal design the stimulus amplitudes ${A_{n}}$ will have an upper bound (numerical details will be discussed Section 3.2. The stimulus amplitude is expected to tend to the upper bound. The parameter estimation process will also be benefited from similar bounds as they are assigned to be positive in (5). In addition assignment of bounds on parameters will aid in prevention of overstimulation of the model.

• 2.

  fmincon and similar algorithms are local optimizers. In order to find a good optimum, the algorithms are often repeated with multiple initial guesses and the best one is chosen according to the value of the objective and gradient value at the termination point. This is especially important in the optimal design part as Fisher Information measure in (21) may have lots of local minima. The same factors do exist for likelihood function (12) however as the number of samples ${M}$ increases the likelihood function tend to converge to the same optimum for different initial guesses. See Section 3.3 for a detailed discussion of this issue.

• 3.

  The local optimization algorithms need the gradient of the objective functions. In MATLAB^®’s fmincon and other similar packages have the ability to compute the gradients numerically. However this capability may not be an advantage in this sort of complicated problems for two main reasons. First of all, numerical gradient computation increases the risk of singularities in the solution. Secondly, the computational duration increases dramatically. It is often noted that, the numerical differentiation in fmincon at least doubles the computational duration. This is unacceptable in the context of this research.

• 4.

  Longer simulation times due to the computational complexity of the overall problem can be eliminated to a certain extent (about the 5 times shorter duration) by utilizing the MATLAB^®’s parallel computation facilities (such as replacing standard for loops by parfor loops).

• 5.

  The MATLAB^® version used in this research is R2013b and the fmincon version included in this package has the interior-point method of constrained non-linear optimization by default. This is not the only method available in fmincon and it can always be modified by changing the algorithm option. However, the interior-point method appears to be the most efficient among all available options.

In this section, we will summarize the overall procedure to give an insight about how the optimal design and parameter estimation algorithm works together in a automatized loop. Before giving a start, it is worth to discuss the initial parameter problem. In the beginning of an experiment (or simulation in our case), one usually has no idea about the true or approximate values of the network parameter vector $\theta$. However, the optimal design always needs a current estimate. Thus one should assign a randomly chosen initial parameter vector. A good rule to get that value is to draw one sample from a parametric subspace of which members are uniformly distributed between the lower and upper bounds of the parameters. One can find the details about the parametric bounds in Table 3 Setting a parametric bound on the stimulus amplitude coefficients ($A_{n}$ in (15)) might be critical. Too large amplitudes may lead to instabilities or very large responses which may break the optimal design procedures. Because of this fact, an upper bound on $A_{n}$’s will most often required. A few number of initial simulations will be required to determine a nominal value for those bounds. If no instabilities are detected one will be fine with the decided value of the bounds on $A_{n}$ (which is termed as $A_{\max}$) and can continue the optimization with these bounds. Mostly, the optimization procedure will fail when the assigned value of the bounds are too large.

The algorithmic details of the overall process are summarized as shown below:

1. 1.

   Set $i=1$

2. 2.

   Generate a random parameter vector $\theta^{0}$ from a uniform distribution between $\theta_{min}$ and $\theta_{max}$ (from Table 3).

3. 3.

   Set $\hat{\theta}=\theta^{0}$

4. 4.

   Set $k=1$

5. 5.

   Optimize $U_{k}$ based on $\hat{\theta}$ to generate stimulus $\bf x$

6. 6.

   Using the stimulus $\bf x$ perform a maximum likelihood estimation to find a new estimate of $\theta$ and replace $\hat{\theta}$ by the new estimate.

7. 7.

   Set $k\rightarrow{k+1}$ and go to Step 5 . If $k>8$ set $i\rightarrow{i+1}$

8. 8.

   If $i>N_{itr}$ stop and report the result as $\hat{\theta}^{final}=\hat{\theta}$ otherwise go to Step 4

This section is devoted to the detailed presentation of the simulation set-up. An numerical example will be presented which will demonstrate our optimal design approach. In the example application, the algorithms presented in Sections 2.3 and 2.4 are applied to probe an EI network.

Our model in (5) has two more important components which are the gain functions $g_{e}\left(V_{e}\right)$ and $g_{i}\left(V_{i}\right)$. These are obtained by setting $j$ in (2) by either ’$e$’ or ’$i$’. So one has $6$ additional parameters $\left[\Gamma_{e},a_{e},h_{e},\Gamma_{i},a_{i},h_{i}\right]$ which have direct effect on the neural model behaviour. This research targeted the estimation of the network parameters only. Because of that, the parameters of the gain functions are kept as fixed and they have the values $\Gamma_{e}=100,a_{e}=0.04,h_{e}=70,\Gamma_{i}=50,a_{i}=0.04,h_{i}=35$.

It is well known that optimization algorithms such as fmincon requires an initial guess of the optimum solution. A suitable choice for the initial guesses can be their assignment from a set of initial conditions generated randomly between the optimization bounds. In the optimization of stimuli, the initial amplitudes can be uniformly distributed between $\left[0,A_{max}\right]$ and phases can be uniformly distributed between $\left[-\pi,\pi\right]$. Although we do not have any constraints on the phase parameter, we limit the initial phase values to a safe assumed range. In the analysis of the optimal phase results we will wrap the resultant phase values to a range $\left[-\pi,\pi\right]$ using modulo function. Thus assuming an initial range in the same interval will be meaningful.

In Section 2.3, one recall from (11) that the likelihood estimation should produce better results when the number of samples (i.e. $M$ in (11)) increases. Because of this fact, the likelihood function is based on data having all spikes generated since the beginning of the simulation. The number of repeats determines $M$. If simulation is repeated $N_{itr}$ times ($N_{itr}$ iterations), one will have an $M$ value of $M=8N_{itr}$ due to the fact that each iteration has $8$ optimal designs sub-steps (with respect to each parameter $\theta_{k}$. Read Section 2.4). So, if one has 15 iterations $N_{itr}=15$, $M=120$ which means likelihood has $120$ samples. This also means that optimal design and subsequent parameter estimation will also be repeated $120$ times.

Having all necessary information from Section 3.1, one can perform an optimal design and obtain a sample optimal stimulus and associated neural responses Figure 3. It is noted that the optimal stimulus in Figure 3 top panel has a periodic variation as it is modelled as a phased cosines form (as it is equivalent to real valued Fourier series).

In addition to the fundamental responses, the second half of Figure 3 displays the variation of the parametric sensitivity derivatives $\frac{\partial{r_{e}}}{\partial{\theta_{i}}}$ which are generated by integrating (28) and then substituting to (31). The variation of the sensitivity derivatives support the idea of optimization of the Fisher Information Metric with respect to a single parameter (see (22)) as the sensitivity (or gradient) of the Fisher Information Metric in (22) with respect to a single parameter $\theta_{k}$ varies widely from parameter to parameter. . It appears that, the weight parameters have a very high sensitivity which are more than $10$ times that of the $\beta_{*}$ parameters. Also some of the parameters affecting the behaviour of inhibitory (I) unit $\beta_{i},w_{i}$ and also the E-I interconnection weights $w_{ei},w_{ie}$ have a reverse behaviour. This fact appears as a negative values variation of the sensitivity derivatives. Self inhibition coefficient $w_{ii}$ does not show this behaviour as it represent an inhibitory effect on the inhibitory neuron potential (which favours excitation).

One of the most distinguishing result related to optimal design is the statistics of the optimal stimuli (i.e. the amplitudes $A_{n}$ and phases $\phi_{n}$). This analysis can be performed by generating the histograms of $A_{n}$ and $\phi_{n}$ from available data as shown in Figure 4. The Figure 4A shows the flat uniform distribution of the random stimuli amplitudes and phases. This is expected as the stimuli is generated directly from a uniform distribution. Concerning the optimal stimuli, the histograms shown in Figure 4B reveals that optimal design has a tendency to maximize the amplitude of the stimuli towards the upper bound. This reassures that optimal design tends to maximize the stimulus power which is expected increase the efficiency of the parameter estimation process and it distinguishes optimal stimuli from their randomly generated counterparts.

The results mentioned above coincide the findings of [35]. Here, what happens is related to the topological boundary property [36] associated with the optimal design results. It is found in [35] that, maximum firing rate response is always lies on the topological boundary of the collection of all allowable stimuli provided that the neurons have increasing gain functions and convergent synaptic connections between layers. Another interpretation of this result can be expressed as follows: The maximum and minimum responses of each individual neuron should arise from the stimulus boundary (from all available sets of stimuli) and the entire boundary of the pattern of responses elicited in a particular layer should also a result of the stimulation in the boundary level. This result is valid regardless of the neuron being feed-forward or recurrent. For recurrent neural networks we assume that the network is stable. Some mathematical proofs related to the stated results are also presented in [35]. In Figure 4B, the amplitudes of the Fourier stimulus components are collected at the boundary level which is $A_{\max}=120$. This is the topological boundary of the set of stimuli defined by (15).

It would also be convenient to stress the fact that the greater the likelihood, the better the fitting of the data to the model. This fact is demonstrated by two regression lines imposed on the box diagram Figure 5. One of those lines correspond to the optimal design and the other correspond to the random stimuli based tests. Both regression lines path through the origin point (0, 0) approximately. This means that the ratio of the log likelihood values in the two cases is approximately a constant, regardless of the number of stimuli. The regression lines are represented by equations $l=257.0896M+76.1500$ for optimal and $174.5375M+31.1000$ for the random stimuli. So the two lines have a slope ratio of approximately $1.4730$. The significance of this number can be explained by a simple example. If one desires to attain the same level of likelihood with $M=120$ optimally designed stimuli, the required number of random stimuli to be generated is equal to a value about $180$.

Another statistical comparison of the maximized likelihoods can be performed by Wilcoxon rank-sum tests [37, 38, 39]. When applied, one will be able to see the difference which is highly significant. Regardless of the number of samples the p-values remained at least $10^{-30}$ times smaller than the widely accepted probability significance threshold of $p=0.05$ or $5\%$.

The likelihood function provides an overall measure of how well a model fits the data. We have also tested the mean errors of individual parameters relative to their true values. The main finding is that, for each individual parameter, the error is typically smaller for the optimally designed stimuli than the error for the random stimuli. This result can be observed from the bar charts presented in Figure 6. The heights of the bars show the mean error levels for randomly and optimally generated stimuli respectively. One can get the benefits of the rank-sum test on the statistical properties of the parameter estimates. Computation of rank-sum p-values for each individual parameter corresponding to the case of $120$ samples $\left(M=120\right)$ yields: