Mathematical model & model estimation

Our objective was to find a mathematical description able to accurately capture the spiking probability of RGCs to subretinal stimulation using a MEA. We characterized neural responses using an N-dimensional linear subspace of the stimulus space, combined with a nonlinearity describing the intrinsic nonlinear firing properties of neurons. Using STC analysis, we derived the lower dimension stimulation subspace that led to a short-latency response in the neuron. By projecting the raw and spike-triggered stimuli onto the lower dimension subspace, we estimated the intrinsic nonlinearity. The probability of generating a spike, given stimulus St→, was estimated as

where N represents the static nonlinear function operating on arguments in μA and v→i (i = 1,2,…,N) represent the N significant principal components. To find v→i (i = 1,2,…,N), the stimulus data were first separated into a matrix containing only stimuli generating a short-latency response, S_D, and a matrix containing all stimuli, S_T (Fig 3A). We found the low-dimensional linear subspace that best captured the difference between the spike-triggered stimuli and the raw ensemble by performing principal component analysis (PCA) on the covariance matrix of the spike-triggered ensemble,

and comparing it to the variance of the raw ensemble which was approximately σ² in all stimulus directions due to the Gaussian nature of S_T. PCA on C_s produce a set of eigenvectors giving a rotated set of axes in stimulus space and a corresponding set of eigenvalues giving the variance of the spike-triggered ensemble along each of the axes. Eigenvalues that are greater than the variance of the input stimuli reveal the directions where the spike-triggered stimuli have experienced an increase in variance from the raw ensemble. Similarly, eigenvectors that are smaller than the variance of the input stimuli reveal directions where the spike-triggered stimuli have experienced a decrease in variance from the raw ensemble. The eigenvalues that are sufficiently different from the raw ensemble correspond to eigenvectors (v→i,i=1,2,…,N) pointing in directions in the stimulus space that carry information about the spiking probability of the neuron.

(a) Discretized sequence of the neural response and stimulus. Each stimulus consists of a combination of biphasic pulses applied to all 20 electrodes. Stimuli that evoked a spike in the neuron are recorded in the stimulus matrix S_D. (b) STC was conducted on the stimuli generating a response, S_D, to separate the stimulus space into a positive and negative region (+ and ×). The x-axis corresponds to the first eigenvector (v→1); the y-axis corresponds to the second eigenvector (v→2). Not all stimuli generated a response in the neuron. Shown in black are the total applied stimuli, S_T, which are overlaid by stimuli S_D (white crosses). Also shown are the projections of the electrical receptive fields, w→+ (large diamond) and w→− (large circle).

To test if the neural response could be accurately characterized by a one-dimensional model, we examined how many eigenvalues resulting from PCA were significantly different to chance [20]. We compared the eigenvalues obtained by PCA on spike-triggering stimuli to a distribution of eigenvalues for a randomly chosen ensemble of stimuli. This was done by randomly time-shifting the spike train and performing PCA on the corresponding randomized spike-triggered stimuli to recover a new set of eigenvalues. By repeating these 1000 times, we construct a distribution of eigenvalues and set a confidence criterion outside of which we presumed the magnitude of the true eigenvalues to be significant. The confidence criterion used was two standard deviations, or a 95% confidence interval. If the greatest or least eigenvalue fell outside these bounds, we rejected the null hypothesis that the spike-triggered ensemble was not significantly different to the full ensemble. We then projected out the axis corresponding to this eigenvalue from the raw data. We repeated the test until all remaining eigenvalues fell within the bounds of the null distribution, suggesting that the remaining eigenvalues were insignificant in affecting the variance of the neuron. Components having an eigenvalue significantly greater than the variance of the randomly time-shifted ensemble were considered to be components that generate an excitatory response on the cell. Conversely, components that are significantly smaller than the variance of the raw ensemble were considered to be components that suppressed the cell’s response.

For the majority of cells, we found that a simplification to one-dimension (v→1) accurately captured the spike-triggering information, thereby reducing the equation to one dimension. Using this simplification, Eq (2) becomes

Results in the literature indicate that response thresholds to electrical stimulation for some cell types might differ depending on pulse polarity [37]. To explore difference in response to pulse polarity, we allowed the probability to be described by two different nonlinear functions and we found the electrical receptive fields (ERFs) for stimuli having a net effect that was either cathodic-first or anodic-first. Eq (4) then becomes

where N+ and N− represent static nonlinear functions and w→+ and w→− represent the ERFs for stimuli with positive projections (v→1∙St→>0, net anodic-first) and negative (v→1∙St→<0, net cathodic-first), respectively. R_s represents the spontaneous firing rate in Hz and c represents a scaling factor of units Hz^-1. To find the nonlinearities and the ERFs, the first principal component (v→1) was used to divide the stimulus space into positive and negative regions by projecting all stimuli of S_D and S_T onto the first principal component (Fig 3B). Positive and negative regions were defined by the stimuli having either a positive or negative projection onto the first principal component. This produced two spike-triggered stimulus matrices, SD+ and SD−. The means of the matrices are analogous to the spike-triggered average for net anodic-first and net cathodic-first stimuli [16], and provide an estimate of the ERFs, w→+ and w→−, respectively. Fig 3B shows an example of the stimuli projected onto the first two principal components and the ERFs, w→+ and w→−. After the stimuli were separated into two regions, the nonlinear functions, N+ and N−, were recovered by projecting all stimuli onto the normalized vectors w→+ and w→− and segmenting the projected stimuli into 30 bins (15 for the net anodic-first and 15 for the net cathodic-first regions) such that each bin contained an equal number of spikes. By comparing the number of spike-eliciting stimuli to the total number of stimuli in each bin, an estimate of the spike probability was generated. The nonlinear function from Eq (5) was then fit to the data, with the following equations for the sigmoidal curves:

where x+=w→+∙St→ and x−=w→−∙St→. Coefficients a₊ and a₋ represent scaling factors that determine the saturation amplitudes, b₊ and b₋ represent the gain of the sigmoidal curves, and c₊ and c₋ represent the thresholds (50% of the saturation level) for the net anodic-first and net cathodic-first stimulation, respectively. Note that the vectors w→+ and w→− might not necessarily be parallel to each other, nor parallel to v→1. This may result in electrodes that differentially influence the neuron’s response to anodic-first or cathodic-first stimulation.

To test the similarity between the positive and negative ERFs, we calculated the correlation coefficient between them. A correlation coefficient close to -1 indicated that the two ERFs are approximately equal in magnitude but opposite in sign, and therefore the cell was equally influenced by the same combination of electrodes. A value closer to 0 indicates that the two ERFs have no correlation, and therefore the cell is not influenced by the same electrodes. Positive correlation coefficients were not expected and did not occur. The spatial extent over which a cell was influenced by electrical stimulation was estimated by computing a weighted mean of the distance between the cell and the electrodes. The distance between the cell and each electrode center was weighted by the electrode’s influence on the cell as given by the ERFs. The weighted mean for both ERFs was given by,

where wi+ and wi− are the weights given by w→+ and w→− respectively, and d_i is the distance between the cell and electrode i.

To test which electrodes significantly affected the cell's response, w→+ and w→− were recalculated 1000 times by projecting the data onto the first eigenvector of the randomly time-shifted distribution of eigenvectors from the significance test. A distribution for w→+ and w→− was constructed from which the true w→+ and w→− could be compared. Electrodes from the true w→+ and w→− were compared to the root mean square (RMS) of the distribution and electrodes that were larger than the RMS bounds were considered significant.

For cells where more than one significant principal component was obtained from the significance test, we compared the variance explained by the first principal component to that of the next most significant component. This was done by comparing the separation of the first eigenvalue e₁ from the mean of the randomized distribution of eigenvalues (e¯rnd) with the separation between the next most significant eigenvalue (e₂) and the same mean. The strength was defined as

and gives a relative measure of how much larger e₁ is compared to the next most significant eigenvalue. e¯rnd was calculated from the first iteration of the hypothesis test.