Global activity shaping strategies for a retinal implant

Accurate quantifiable analysis of alternative global shaping strategies initially requires a forward model to calculate retinal activations resulting from a particular pattern of electrode settings. The model will use scaled units with S—spikes, T—time, I—electric current, L—length.

We first define a desired set of retinal activations, $\vec {r}_{{\rm d}}^{*}[{\rm S}{{{\rm T}}^{-1}}]$, and a modeled set of retinal activations ${{\vec {r}}_{{\rm d}}}~[{\rm S}{{{\rm T}}^{-1}}]$ that result from a given set of electrode current amplitudes, $\vec {s}~[{\rm I}]$. The retinal activation vectors $\vec {r}_{{\rm d}}^{*}$ and ${{\vec {r}}_{{\rm d}}}$ represent the set of values associated with ${{N}_{{\rm d}}}$ discrete retinal locations. The current amplitude vector, $\vec {s}$, represents the values associated with ${{N}_{{\rm e}}}$ electrode locations. These sets of retinal activations and current amplitudes could be spatially configured in an arbitrary way, for example the set of electrodes could be configured as hexagonal arrays. However, in the present investigation, the retinal and electrode configurations take the form of square arrays of dimensions ${{n}_{{\rm d}}}$, and ${{n}_{{\rm e}}}$ respectively, such that ${{N}_{{\rm d}}}=~n_{{\rm d}}^{2}$, and ${{N}_{{\rm e}}}=~n_{{\rm e}}^{2}$ and the spatial width of the array is ${{L}_{{\rm e}}}[L]$ is set to be 1. A ground electrode is assumed to be distal, as is the case for most retinal implants. We also define the generator signal, $\vec {G}~\left[ I \right]$ which is an intermediate variable used to calculate retinal activations, with the same length, ${{N}_{{\rm d}}}$, as the retinal activation vector. In the present investigation, the generator signal at any given location in the retinal tissue is the total local current at that location, which is the sum of the currents resulting from each electrode activation.

The forward model that converts electrode amplitudes, $\vec {s}$, to retinal activations, ${{\vec {r}}_{{\rm d}}}$, is defined in two steps: first, a linear step that maps the vector of electrode amplitudes, $\vec {s}$, to a vector of generator signal values, $\vec {G}$,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \vec {G}=g\left( {\vec {s}} \right),\nonumber \end{align} \tag{ 1 }$$

And, second, a nonlinear step that acts on each element of the generator vector, $\vec {G}$, to give the retinal activation vector, ${{\vec {r}}_{{\rm d}}}$,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{r}_{{\rm d},i}}=f\left( {{G}_{i}} \right),i=1,\ldots ,{{N}_{{\rm d}}}.\nonumber \end{align} \tag{ 2 }$$

In the retina, the biophysical basis for the utility of a linear-nonlinear model and electrical receptive field have been demonstrated experimentally [9, 26, 29–33] and loosely correspond to the subthreshold membrane current, spike generation threshold, and field of network inputs, respectively. The functions $f$ and $g$ are defined below.

The function $g$ in equation (1) is a transfer matrix:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \vec {G}=g\left( {\vec {s}} \right)={{{\bf W}}_{{\bf d}}}~\vec {s},\nonumber \end{align} \tag{ 3 }$$

where, ${{{\bf W}}_{{\bf d}}}$ is a matrix of dimensions ${{N}_{{\rm e}}}\times {{N}_{{\rm d}}}$, and gives the magnitude of the current at each retinal location resulting from the activation of each electrode. The rows of the matrix ${{{\bf W}}_{{\bf d}}}$ are the ERFs for each retinal patch, and the columns contain the retinal activation spread for each stimulating electrode. Gaussian electrode activation patterns were chosen to approximate the activation spread resulting from each electrode [12], with each element, ${{W}_{{\rm d},ij}}$, given by:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{W}_{{\rm d},ij}}=\left( 1+\rho {{\xi }_{ij}} \right){{e}^{\frac{-{{\left( {{x}_{i}}-{{x}_{j}} \right)}^{2}}-{{\left( {{y}_{i}}-{{y}_{j}} \right)}^{2}}}{2{{\sigma }^{2}}}}},\nonumber \end{align} \tag{ 4 }$$

where $({{x}_{i}},{{y}_{i}})~[L,L]$ is the spatial positions of the $i$th retinal patch, ${{\vec {r}}_{{\rm d},i}}$, and $({{x}_{j}},{{y}_{j}})~[L,L]$ the jth electrode, ${{\vec {s}}_{j}}$; $\sigma ~[L]$ is the spread of current, $\rho$ is the magnitude of variation in the activation spread, and ${{\xi }_{ij}}$ is a random number equally distributed between  −1 and 1. The additive variation, $\rho {{\xi }_{ij}}$, represents deviations in the Gaussian shape due to heterogeneity in the retinal tissue.

In its general form, the function $f$ (equation (2)) is the two-sided sigmoid function illustrated by the red line in figure 1(a) [9, 28]. With appropriate limits placed on the range of retinal activations and by appropriate linear scaling, the function $f$ can be approximated by a modulus $\left| {{G}_{i}} \right|$, as illustrated by the solid black line in figure 1(a). Using this approximation and equation (3), the final forward model described in equations (1) and (2) is now given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\vec {r}}_{{\rm d}}}=k\cdot |{{{\bf W}}_{{\bf d}}}\vec {s}|,\nonumber \end{align} \tag{ 5 }$$

where $k[{\rm S}{{{\rm T}}^{-1}}\,{{{\rm I}}^{-1}}]$ is set to 1. This forward model is illustrated in one dimension in figures 1(b)–(e). However, in general, these spatial interactions take place in two-dimensions, giving the retinal activation pattern.

Zoom In Zoom Out Reset image size

The units of the elements of each the vectors, ${{\vec {r}}_{{\rm d}}}$ and $\vec {s}$, and the matrix, ${{{\bf W}}_{{\bf d}}}$, were scaled such that an electrode stimulus amplitude value of $1\,[{\rm I}]$ results in a peak generator signal value of $1\,[{\rm I}]$, which in turn results in a retinal activation of $1\,[{\rm S}{{{\rm T}}^{-1}}]$ (figure 1(a)). In real implants, there is a limit on the maximum current $\vec {b}\,[{\rm I}]$ deliverable by each electrode for both power and safety reasons.

There are two regimes in which the model can operate: the 'narrow' activation spread regime, in which the influence of neighboring electrodes is isolated and does not significantly overlap, and the 'wide' activation spread regime, in which the retinal activation at any given point is the result of significant contributions from multiple electrodes. In the present investigation, the overlap parameter, $\gamma$, was defined by the ratio between one standard deviation in activation spread and the distance between neighboring electrodes in the square array, $\gamma =\frac{\sigma }{{{L}_{{\rm e}}}/{{n}_{{\rm e}}}}=\frac{\sigma {{n}_{{\rm e}}}}{{{L}_{{\rm e}}}}$, and since we choose units with ${{L}_{{\rm e}}}=1$, then $\gamma =~\sigma {{n}_{{\rm e}}}$. This parameter was used to quantify the overlap, where $\gamma =1$ approximately denotes the boundary between the narrow and wide activation spread regimes. The quantity, $\gamma$, and its dependence on $\sigma$, and ${{n}_{{\rm e}}}$ is illustrated in figure 1(f).

The next stage in the development of the model requires derivation of an inverse model as illustrated in figure 2(a). From the forward model's matrix ${{{\bf W}}_{{\bf d}}}$ equation (5), (figure 2(b)) of size ${{N}_{{\rm d}}}\times {{N}_{{\rm e}}}$ we derive 'measured' activation spreads ${{{\bf W}}_{{\bf t}}}$ (figure 2(c)) of dimension ${{N}_{{\rm t}}}\times {{N}_{{\rm e}}}$. This is a lower resolution and noisy version of ... The matrix ${{{\bf W}}_{{\bf t}}}$ is intended to represent an experimentally obtained ERF model and is the matrix used in the inverse model while the matrix ${{{\bf W}}_{{\bf d}}}$ is intended to model the actual retinal activation that results from a given pattern of electrodes. Separating these two models is a conservative check on the results obtained by the methods proposed in this investigation. ${{{\bf W}}_{{\bf t}}}$ was modeled by downsampling each of the 2D electrical receptive fields contained in ${{{\bf W}}_{{\bf d}}}$ from ${{N}_{{\rm d}}}$ elements to ${{N}_{{\rm t}}}$ elements using the MATLAB 2017a (Mathworks, MA, USA) function 'interp2', and then simulating additive noise by adding a quantity $\newcommand{\e}{{\rm e}} \eta {{\xi }_{ij}}$ to each element of the resulting activation spreads, where $\newcommand{\e}{{\rm e}} \eta$ is the magnitude of the 'estimation noise' and ${{\xi }_{ij}}$ is a random number between  −1 and 1. This quantity represents errors in measuring the true activation spreads. This estimation noise is critical in assessing the quality of the algorithms because, unlike the conventional method, shaping methods are susceptible to over-fitting under certain conditions, which can lead to artifacts in the resulting retinal activation pattern.

Zoom In Zoom Out Reset image size

All parameters are defined in table 1.

Table 1. List of parameters and their units of measurement. T—time, L—length, S—spikes, V—voltage, I—current.

Parameters
Electrode amplitudes (${{N}_{{\rm e}}}$)	$\vec {s}\,[{\rm I}]$
Maximum electrode amplitude	${{s}_{\max }}\,[{\rm I}]$
Electrode array width, 1	${{L}_{{\rm e}}}=1\,[{\rm L}]$
'True' transfer matrix (${{N}_{{\rm d}}}\times {{N}_{{\rm e}}}$)	${{{\bf W}}_{{\bf d}}}\,[{\rm S}{{{\rm T}}^{-1}}\,{{{\rm I}}^{-1}}]$
Desired retinal array dimensions	${{n}_{{\rm d}}}$
Desired retinal vector length, $n_{{\rm d}}^{2}$	${{N}_{{\rm d}}}$
Desired retinal activations (${{N}_{{\rm d}}}$)	$\vec {r}_{{\rm d}}^{*}\,[{\rm S}{{{\rm T}}^{-1}}]$
Modelled desired retinal activations (${{N}_{{\rm d}}}$)	${{\vec {r}}_{{\rm d}}}\,[{\rm S}{{{\rm T}}^{-1}}]$
Mean squared error, $\left\Vert \vec {r}_{{\rm d}}^{*}-\vec {r}_{{\rm d}}^{2} \right\Vert$	${{\varepsilon }_{{\rm d}}}\,[{{{\rm S}}^{2}}\,{{{\rm T}}^{-2}}]$
Current spread	$\sigma \,[{\rm L}]$
Variation magnitude in ${{{\bf W}}_{{\bf d}}}$	$\rho$
Target transfer matrix (${{N}_{{\rm t}}}\times {{N}_{{\rm e}}}$)	${{{\bf W}}_{{\bf t}}}$
Target retinal array dimensions	${{n}_{{\rm t}}}$
Target retinal vector length, $n_{{\rm t}}^{2}$	${{N}_{{\rm t}}}$
Target retinal activations (${{N}_{{\rm t}}}$)	$\vec {r}_{{\rm t}}^{*}\,[{\rm S}{{{\rm T}}^{-1}}]$
Modelled target retinal activations (${{N}_{{\rm t}}}$)	${{\vec {r}}_{{\rm t}}}\,[{\rm S}{{{\rm T}}^{-1}}]$
Mean squared error, $\left\Vert {{{\vec {r}}}_{{\rm t}}}^{*}-{{{\vec {r}}}_{{\rm t}}}^{2} \right\Vert$	${{\varepsilon }_{{\rm t}}}\,[{{{\rm S}}^{2}}\,{{{\rm T}}^{-2}}]$
Estimation noise magnitude in ${{{\bf W}}_{{\bf t}}}$	$\newcommand{\e}{{\rm e}} \eta$
Binary partitioning vector (${{N}_{{\rm t}}}$)	$\vec {\varphi }$
Diagonal matrix form of $\vec {r}_{{\rm t}}^{*}$ (${{N}_{{\rm t}}}\times {{N}_{{\rm t}}}$)	${{{\bf R}}^{*}}\,[{\rm S}{{{\rm T}}^{-1}}]$
Maximum allowed electrode amplitudes	$\vec {b}\,[{\rm I}]$
Electrode array dimensions	${{n}_{{\rm e}}}$
Electrode vector length, $n_{{\rm e}}^{2}$	${{N}_{{\rm e}}}$
Neural activation constant	$k=1\,[{\rm S}{{{\rm T}}^{-1}}\,{{{\rm I}}^{-1}}]$
Electrode locations	$\left( {{x}_{j}},~{{y}_{j}} \right)\,[{\rm L},{\rm L}]$
Retinal locations	$\left( {{x}_{i}},~{{y}_{i}} \right)[{\rm L},{\rm L}]$
Linear model	$g$
Nonlinear model	$f$
Random value	$\newcommand{\e}{{\rm e}} \xi \epsilon (-1,1)$
Generator function values (${{N}_{{\rm t}}}$)	$\vec {G}\,[{\rm I}]$
Current overlap parameter, $\sigma {{n}_{{\rm e}}}/{{L}_{{\rm e}}}$	$\gamma$
Linear inversion objective function	${{\Theta }_{{\rm lin}}}\,[{{{\rm S}}^{2}}\,{{{\rm T}}^{-2}}]$
Quadratic program objective function	${{\Theta }_{{\rm quad}}}~[{{{\rm S}}^{2}}{{{\rm T}}^{-2}}]$
Binary partition objective function	${{\Theta }_{{\rm bin}}}\,[{{{\rm S}}^{2}}\,{{{\rm T}}^{-2}}]$
Eigenvalue fraction	${{F}_{{\rm eig}}}$
Left eigenvectors, retina (${{F}_{{\rm eig}}}{{N}_{{\rm e}}}\times {{N}_{{\rm t}}}$)	${\bf P}$
Right eigenvectors, electrode (${{N}_{{\rm e}}}\times {{F}_{{\rm eig}}}{{N}_{{\rm e}}}$)	${\bf Q}$
Eigenvalues (${{F}_{{\rm eig}}}{{N}_{{\rm e}}}\times {{F}_{{\rm eig}}}{{N}_{{\rm e}}}$)	${\bf D}$
Pseudoinverse of ${{{\bf W}}_{{\bf t}}}$ (${{N}_{{\rm e}}}\times {{N}_{{\rm t}}}$)	${\bf W}_{{\bf t}}^{+}$
Mean error over 32 images	${{\varepsilon }_{{\rm d}32}}\,[{{{\rm S}}^{2}}\,{{{\rm T}}^{-2}}]$
Partitioned retinal activity, ${{{\bf R}}^{*}}\vec {\varphi }$ (${{N}_{{\rm t}}}$)	$\vec {p}_{{\rm t}}^{*}\,[{\rm S}{{{\rm T}}^{-1}}]$

The mean squared error, ${{\varepsilon }_{{\rm t}}}\,[{{{\rm S}}^{2}}\,{{{\rm T}}^{-2}}]$, (figure 2(b)) quantifies the magnitude of the error between a target retinal activation, $\vec {r}_{{\rm t}}^{*}$, and the retinal activation resulting from electrical stimulation, ${{\vec {r}}_{{\rm t}}}$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\varepsilon }_{{\rm t}}}=\left\Vert \vec {r}_{{\rm t}}^{*}-\vec {r}_{{\rm t}}^{2} \right\Vert.\nonumber \end{align} \tag{ 6 }$$

This is used as an objective function, and each inverse model is developed by attempting to minimize this objective function under different assumptions.

As a test of the results of each inverse model, the electrode patterns are applied using equation (5), and the resulting pattern compared to the desired pattern to give ${{\varepsilon }_{{\rm d}}}\,[{{{\rm S}}^{2}}\,{{{\rm T}}^{-2}}]$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\varepsilon }_{{\rm d}}}=\left\Vert \vec {r}_{{\rm d}}^{*}-\vec {r}_{{\rm d}}^{2} \right\Vert.\nonumber \end{align} \tag{ 7 }$$

This process is illustrated in figure 2(c). The procedure is intended to replicate the real conditions by confirming that elements of the retina that are not included in ${{{\bf W}}_{{\bf t}}}$ remain constrained by the inverse models that we develop.

In this investigation, 32 images were used as the basis for the desired retinal activation patterns. These images are shown in figure 3. The mean value of ${{\varepsilon }_{{\rm d}}}$ over over these 32 images was defined as ${{\varepsilon }_{{\rm d}32}}$. The images were downsampled to give the desired retinal activation points ${{N}_{{\rm d}}}$, and the target retinal activation points ${{N}_{{\rm t}}}$. This was done using the MATLAB 2017a (MathWorks, MA, USA) function 'imresize' with the 'antialiasing' option, which uses a low-pass filter adapted to the resolution of the desired array dimensions.

Zoom In Zoom Out Reset image size

The conventional approach to calculating electrode activations in retinal implants is to simply downsample the desired activation pattern to give an array of retinal activation 'pixels' with the same dimensions as the electrode array. This is completed using the MATLAB function 'imresize' with the 'antialiasing' option. These electrode values, which are between 0 and 1 may lead to oversaturation of the retina, especially in the case of wide activation spread ($\gamma >1$). To normalize the retinal activation, the initial electrode settings are provided to the forward model to calculate an initial set of retinal activations. The ratio of the resulting total mean retinal activation and the target total mean retinal activation is then used to scale the initial set of electrodes.

In calculating the inverse model we make an approximation in which the modulus in equation (5) can be removed, and choose $k=1$. The inverse model will then operate using the assumption that ${{r}_{{\rm i}}}={{G}_{{\rm i}}}$, the dashed line in figure 1(a). The consequences of this approximation are discussed below. Replacing the true set of activation spreads, ${{{\bf W}}_{{\bf d}}}$, with the measured set of activation spreads, ${{{\bf W}}_{{\bf t}}}$, this approximation means that equation (5) becomes ${{\vec {r}}_{{\rm t}}}={{{\bf W}}_{{\bf t}}}\vec {s}$. The error, ${{\varepsilon }_{{\rm t}}}$ equation (6), can then be re-written as an objective function, ${{\Theta }_{{\rm lin}}}[{{{\rm S}}^{2}}\,{{{\rm T}}^{-2}}]$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\Theta }_{{\rm lin}}}=\underset{{{{\vec {s}}}_{i}}\epsilon \mathbb{R}}{\mathop{\min }}\,~{{\left\Vert {{{\vec {r}}}_{t}}^{*}-{{{\bf W}}_{{\bf t}}}\vec {s} \right\Vert}^{2}}.\nonumber \end{align} \tag{ 8 }$$

The solution of this equation, $\vec {r}_{{\rm t}}^{*}={{{\bf W}}_{{\bf t}}}\vec {s}$, can be arranged to give $\vec {s}={\bf W}_{{\bf t}}^{-1}\vec {r}_{{\rm t}}^{*}$. However a true inverse of the matrix ${{{\bf W}}_{{\bf t}}}$ is not possible because in general it is a non-square matrix, and we must use a pseudo-inverse, ${\bf W}_{{\bf t}}^{+}$.

The pseudo-inverse of ${{{\bf W}}_{{\bf t}}}$ is calculated using singular value decomposition (SVD). Under SVD, ${{{\bf W}}_{{\bf t}}}={\bf PD}{{{\bf Q}}^{T}}$, where the matrices P and Q are unitary and contain the left and right eigenvectors of W_t; these are associated with the generator signal pattern and the electrode pattern, respectively. D is a rectangular diagonal matrix and contains the eigenvalues associated with those eigenvectors. The pseudo-inverse of W_t can then be calculated as ${\bf W}_{{\bf t}}^{+}={\bf Q}{{{\bf D}}^{-1}}{{{\bf P}}^{T}}$, and the electrode activation currents are given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \vec {s}={\bf W}_{{\bf t}}^{+}\vec {r}_{{\rm t}}^{*}.\nonumber \end{align} \tag{ 9 }$$

The number of eigenvectors used in the inversion can be varied, as explained in the results, and the fraction of eigenvectors used is defined as ${{F}_{{\rm eig}}}$.

The linear inversion algorithm does not place constraints on electrode amplitudes. High electrode amplitudes can be avoided by using a quadratic program with linear constraints approach, which has the benefit of allowing the imposition of an explicit limit on the electrode current amplitudes.

A second inverse model can be developed by substituting ${{\vec {r}}_{{\rm t}}}={{{\bf W}}_{{\bf t}}}\vec {s}$ in equation (6), and ignoring the constant term $\vec {r}_{{\rm t}}^{*T}\vec {r}_{{\rm t}}^{*}$ in the subsequent quadratic expansion. The error, ${{\varepsilon }_{{\rm t}}}$ equation (6), can then be re-written as an objective function, ${{\Theta }_{{\rm quad}}}[{{{\rm S}}^{2}}\,{{{\rm T}}^{-2}}]$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\Theta }_{{\rm quad}}}~=\underset{{{{\vec {s}}}_{i}}\epsilon (-2,2)}{\mathop{\min }}\,\left( -2\vec {r}_{{\rm t}}^{*T}{{{\bf W}}_{{\bf t}}}\vec {s}+{{{\vec {s}}}^{T}}{\bf W}_{{\bf t}}^{\boldsymbol{T}}{{{\bf W}}_{{\bf t}}}\vec {s} \right).\nonumber \end{align} \tag{ 10 }$$

This equation fits the form of the standard quadratic program minimization, $\left( \min \left( {{c}^{T}}\vec {v}+\frac{1}{2}{{{\vec {v}}}^{T}}{\bf M}\vec {v} \right);{\bf A}\vec {v}<\vec {b} \right)$, with ${\bf M}=2{\bf W}_{{\bf t}}^{\boldsymbol{T}}{{{\bf W}}_{\boldsymbol{t}}}$ and ${{c}^{T}}=-2\vec {r}_{{\rm t}}^{*T}{{{\bf W}}_{{\bf t}}}$. The problem is convex because $2{\bf W}_{{\bf t}}^{\boldsymbol{T}}{{{\bf W}}_{{\bf t}}}$ is positive semi-definite. Each element of the vector of electrode strengths, $\vec {s}$, must be limited to conform to the maximum allowed electrode amplitudes. To achieve this two-sided boundary, A is set to be a pseudo identity matrix with values of 1 in the top left quadrant and values of  −1 in the bottom right quadrant, and $\vec {b}$ is the vector of the limits on the electrode strength magnitudes. With these substitutions, the MATLAB 2017a (Mathworks, MA, USA) function 'quadprog' was used to solve equation (10). This function uses an interior point convex method [34, 35].

It is important to note that removing the modulus from the forward model (equation (5)) in order to derive the two inverse models (the linear inversion model (equation (9)) and the quadratic programming model (equation (10))) prevents these inverse models from finding solutions where negative generator signal could have been utilized to induce retinal activation. A binary quadratic program approach restores the modulus to the inverse model allowing it to discover these negative generator signal solutions.

The modulus in equation (5), which was removed in calculating the linear inversion and quadratic program inverse models, can be retained in the calculation of the inverse model. Using equation (5) with $k=1$, and W_d replaced with W_t, we wish to solve $\vec {r}_{{\rm t}}^{*}=\left| {{\boldsymbol{W}}_{{\bf t}}}\vec {s} \right|$ for $\vec {s}$. This is equivalent to:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{{\bf R}}^{*}}\vec {\varphi }={{{\bf W}}_{{\bf t}}}\vec {s},\nonumber \end{align} \tag{ 11 }$$

where ${{{\bf R}}^{*}}$ is a diagonal matrix containing the vector, $\vec {r}_{{\rm t}}^{*}$, and $\vec {\varphi }={\rm sgn} ({{{\bf W}}_{{\bf t}}}\vec {s})$, which is here referred to as the binary partitioning vector. We use ${{\vec {p}}_{{\rm t}}}$ [${{{\rm S}}^{1}}\,{{{\rm T}}^{-1}}$] as the notation for the partitioned retinal activation resulting from a set of electrode amplitudes, and its target value $\vec {p}_{{\rm t}}^{*}$ [${{{\rm S}}^{1}}\,{{{\rm T}}^{-1}}$] is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \vec {p}_{{\rm t}}^{*}={{{\bf R}}^{*}}\vec {\varphi }={{\vec {\varphi }}^{T}}{{{\bf R}}^{*}}.\nonumber \end{align} \tag{ 12 }$$

The solution to equation (11) is $\vec {s}={\bf W}_{{\bf t}}^{+}{{{\bf R}}^{*}}\vec {\varphi }$, and so the resulting partitioned retinal activation will be

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\vec {p}}_{{\rm t}}}={{{\bf W}}_{{\bf t}}}{\bf W}_{{\bf t}}^{+}{{{\bf R}}^{*}}\vec {\varphi }={\bf PD}{{{\bf Q}}^{\boldsymbol{T}}}{\bf Q}{{{\bf D}}^{-{\bf 1}}}{{{\bf P}}^{\boldsymbol{T}}}{{{\bf R}}^{*}}\vec {\varphi }={\bf P}{{{\bf P}}^{T}}{{{\bf R}}^{*}}\vec {\varphi }.\nonumber \end{align} \tag{ 13 }$$

Using the partitioned form of retinal activation the error can be defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\tilde{\varepsilon }}_{{\rm t}}}=\left\Vert \vec {p}_{{\rm t}}^{*}-\vec {p}_{{\rm t}}^{2} \right\Vert.\nonumber \end{align} \tag{ 14 }$$

Substituting equations (12) and (13) into (14) gives

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\tilde{\varepsilon }}_{t}}={{\left\Vert {{{\vec {\varphi }}}^{T}}{{{\bf R}}^{*}}-{\bf P}{{{\bf P}}^{{\bf T}}}{{{\bf R}}^{*}}\vec {\varphi } \right\Vert}^{2}}=~{{\vec {\varphi }}^{T}}[{{{\bf R}}^{*}}\left( {\bf I}-{\bf P}{{{\bf P}}^{\boldsymbol{T}}} \right){{{\bf R}}^{*}}]\vec {\varphi }.\nonumber \end{align} \tag{ 15 }$$

The minimization of this quantity is then performed over the partitioning vector $\vec {\varphi }$, by introducing the objective function, ${{\Theta }_{\rm bin}}\left[ {{S}^{2}}{{T}^{-2}} \right]:$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\Theta }_{\rm bin}}=~\underset{{{{\vec {\varphi }}}_{i}}\epsilon \{-1,1\}}{\mathop{\min }}\,~{{\vec {\varphi }}^{T}}[{{R}^{*}}\left( I-P{{P}^{T}} \right){{R}^{*}}]~\vec {\varphi }.\nonumber \end{align} \tag{ 16 }$$

Constant elements in this equation that do not vary with sign changes in elements of $\vec {\varphi }$ can be removed from the equation without affecting the minimization. Since $\newcommand{\e}{{\rm e}} {{\vec {\varphi }}_{i}}\epsilon \{-1,1\}$, the diagonal elements of ${{\Theta }_{{\rm bin}}}$ will contain $\vec {\varphi }_{i}^{2}=1$, a constant that will not influence the minimization over $\vec {\varphi }$. This allows us to remove all diagonal elements of ${{\Theta }_{{\rm bin}}}$; i.e. the identity matrix I is removed, and PP^T is replaced with with ${\bf P}{{{\bf P}}^{{\bf T}}}-{\bf diag}\left( {\bf P}{{{\bf P}}^{\boldsymbol{T}}} \right)$ giving

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\Theta }_{{\rm bin}}}=\underset{{{{\vec {\varphi }}}_{i}}\epsilon \{-1,1\}}{\mathop{\min }}\,{{\vec {\varphi }}^{T}}[{{{\bf R}}^{*}}({\bf diag}\left( {\bf P}{{{\bf P}}^{\boldsymbol{T}}} \right)-{\bf PP}{{{\bf }~ }^{{\bf T}}}){{{\bf R}}^{*}}]\vec {\varphi }.\nonumber \end{align} \tag{ 17 }$$

Since $\newcommand{\e}{{\rm e}} {{\vec {\varphi }}_{i}}\epsilon \{-1,1\}$ this equation fits the form of the binary optimization problem, ${{\vec {v}}^{T}}\mathbf{M}\vec {v}$, with ${\bf M}={{{\bf R}}^{*}}({\bf diag}\left( {\bf P}{{{\bf P}}^{\boldsymbol{T}}} \right)-{\bf PP}{{{\bf }~ }^{{\bf T}}}){{{\bf R}}^{*}}$ and $\vec {v}$  =  $\vec {\varphi }$. There are ${{2}^{{{N}_{{\rm t}}}}}$ possible solutions to this equation. The binary optimization problem is NP-hard even when M is positive definite, with solutions often falling into local minima and unable to 'escape' to a globally optimum solution.

Our approach to solving the local minima problem is to use a Hopfield network in which the weights between the neurons are given by the values contained in the matrix M. The initial values of the vector $\vec {\varphi }$ are chosen to be random and then updated according to the rule, ${\rm sgn} ({\bf M}{{\vec {\varphi }}_{k}})\to {{\vec {\varphi }}_{k+1}}$. This process is iterated until the vector $\vec {\varphi }$ is static.

In initial investigations using the Hopfield network, it was observed that solutions frequently fell into local minima, with many unnecessary and detrimental partitions between positive and negative current regions resulting in high value of ${{\varepsilon }_{{\rm d}}}$. To mitigate this problem, simulated annealing was applied. Elements of the vector $\vec {\varphi }$ were randomly 'flipped' on each cycle before that vector was used as the initial conditions of a new Hopfield network optimization. The number of annealing cycles was set to 12 000, with an initial probability that any given element of $\vec {\varphi }$ would be flipped set at 0.5. This probability was linearly reduced so that it reaches $1/N_{{\rm t}}^{2}$ as the iterations reach 12 000. For each of the images this process was repeated 10 times and the solution with lowest value of ${{\tilde{\varepsilon }}_{{\rm t}}}$ selected as the output.

The conventional and linear inversion algorithms were applied to calculate electrode strengths for a particular desired retinal activation pattern (figure 4). Then the forward model equation (5) was used to convert those electrode strengths into their retinal activation patterns (figure 4).

Zoom In Zoom Out Reset image size

The conventional algorithm stimulation pattern consisted of entirely positive values. This is expected given that it is simply a downsampled and intensity-scaled version of the original image. The linear inversion algorithm produced a stimulation pattern with relatively high electrode amplitudes of both polarities. These higher amplitude electrode strengths summate to produce a similar overall level of retinal activation as the conventional strategy, but with sharper contrast. This qualitative difference is the benefit of a global shaping strategy.

In the results shown in figure 4 the linear inversion algorithm utilized all eigenvectors associated with the pseudoinverse of ${{{\bf W}}_{{\bf t}}}$ in equation (9). The use of different numbers of eigenvectors leads to different solutions (figure 5).

Zoom In Zoom Out Reset image size

It can be seen (figure 5(b)) that the use of fewer eigenvectors resulted in a retinal activation pattern of higher error (${{\varepsilon }_{{\rm d}}}=0.031$) and lower spatial frequency, while the use of more eigenvectors created a retinal activation pattern of lower error (${{\varepsilon }_{{\rm d}}}=0.015$) and with features of higher spatial frequency.

Figure 6 explores the effects on the solution of different numbers of eigenvectors. It shows the form of eigenvalues and eigenvectors under four different parameter combinations as a function of ${{F}_{{\rm eig}}}$ (figures 6(a)–(c)). The normalized value of ${{\varepsilon }_{{\rm d}32}}$ under 648 parameter combinations as a function of ${{F}_{{\rm eig}}}$ is also shown (figure 6(d)). These parameter combinations are intended to represent different device conditions as explained below.

Zoom In Zoom Out Reset image size

In principle, the greater the number of eigenvectors used in the calculation of the pseudoinverse the more accurate the result will be. However, in practice there are two limiting factors to the number of eigenvectors used. First, if the matrix W_t is affected by noise, then some eigenvectors will be corrupted, and their inclusion in the calculation of the pseudo-inverse will lead to overfitting. Second, some eigenvectors require very large changes in the amplitude of electrodes to make a significant change in retinal activation. Both these limitations apply to eigenvectors associated with small eigenvalues. Under SVD, the eigenvectors are ordered by the magnitude of their associated eigenvalues. When using a reduced number of eigenvectors we chose to remove those associated with the smallest eigenvalues first.

The left eigenvectors (figure 6(a)) are the retinal activation pattern that results from the electrode amplitude pattern contained in the associated right eigenvectors (figure 6(c)). The eigenvalue magnitudes (figure 6(b)) can be interpreted as the ratio of the amplitude of the associated pair of left and right eigenvectors. In other words, if an eigenvalue is relatively small then the associated electrode pattern (right eigenvector) will have to be of relatively higher amplitude if it is to have as much impact on amplitude of the generator signal pattern. It is also apparent that generator signal patterns with higher spatial frequency are associated with small eigenvalues (figure 6(a)), and as mentioned above, it is the small eigenvalues that are removed first. Given these observations, it is possible to predict that reducing the fraction of eigenvectors used, ${{F}_{{\rm eig}}}$, will lead to solutions of lower spatial frequency and low electrode amplitudes, and this is exactly what was observed (figure 5). In a device with wide activation spreads, $\sigma =0.2$, no variation, $\rho =0$, and no noise, $\newcommand{\e}{{\rm e}} \eta =0$ (parameter combination 4 in figure 6(a)), the lowest magnitude eigenvalues have smaller values than with the other parameter combinations shown. Therefore, using the fact that lower eigenvalues predict higher electrode amplitudes, it can be predicted that linear inversion will create higher amplitude electrode patterns under these conditions. With narrower activation spreads, the smallest eigenvalues increase in magnitude (parameter combination 3 in figure 6(a)) which would reduce the resulting electrode pattern's amplitude. With more variation or measurement noise (parameter combinations 3, and 4 in figure 6(a)) the higher eigenvalues also increase in magnitude which again predicts lower electrode amplitudes.

In the presence of ERF variation or measurement noise, the eigenvectors associated with the smallest eigenvalues become spatially unstructured and noisy so that they can no longer contribute to improving the resolution of the evoked retinal activity pattern (figure 6(a), parameter combination 1). This observation implies that over-fitting in the presence of noise could be avoided by excluding eigenvectors associated with eigenvalues of small magnitude in calculating the pseudoinverse. However, as seen in figure 5, the use of fewer eigenvectors can reduce the accuracy of the solution by placing an upper limit on the spatial frequencies of eigenvectors used. It was therefore necessary to test the effect of using a different fraction of eigenvalues under the full range of parameters (figure 6(d)).

Each of the parameters $\rho$, and $\newcommand{\e}{{\rm e}} \eta$ were given six values covering the ranges, $0.05<\rho <0.3$, and $\newcommand{\e}{{\rm e}} 0.05<\eta <0.3$, and $\sigma$ was given values such that $1<\gamma <16$, giving 216 parameter combinations, and these were tested with three different electrode numbers ${{N}_{{\rm e}}}$  =  2, 8, and 32 making a total of 648 different combinations of the four dimensions. These 648 combinations can be thought of as representing potential devices.

For each of these 648 devices, the mean value of ${{\varepsilon }_{{\rm d}}}$ across 32 images, ${{\varepsilon }_{{\rm d}32}}$, was calculated for each value of ${{F}_{{\rm eig}}}$ and then divided by the smallest ${{\varepsilon }_{{\rm d}32}}$ for all ${{F}_{{\rm eig}}}$ for that combination of parameters. This is referred to as the 'Normalized ${{\varepsilon }_{{\rm d}32}}$'. Each parameter combination can be thought of as creating the conditions for a specific device. The normalization was undertaken because the aim was to find the best choice of ${{F}_{{\rm eig}}}$ for a given device, independent of the absolute value of ${{\varepsilon }_{{\rm d}32}}$ for that device for different values of ${{F}_{{\rm eig}}}$.

It was found that the choice of ${{F}_{{\rm eig}}}$  =  0.6 led to an acceptable result under all devices explored (figure 6(d)). In the worst case device, it led to a 27% increase in the ${{\varepsilon }_{d32}}$ compared to the best choice of ${{F}_{{\rm eig}}}$ for that device (figure 6(d)). This value of ${{F}_{{\rm eig}}}$ was used for all other calculations in the present investigation.

Figure 7 compares the linear inversion and conventional algorithms under different activation spread conditions by varying the values of $\sigma$ and ${{n}_{{\rm e}}}$. The conventional and linear inversion algorithms gave different results for different electrode overlaps. For narrow activation spreads (first row in figure 7), an increase in electrode density for both the conventional and linear inversion algorithms reduced error, ${{\varepsilon }_{{\rm d}}}$, in the resulting retinal activation patterns. However, the same does not apply under wide activation spreads (last row, figure 7).

Zoom In Zoom Out Reset image size

The conventional algorithm's error value does not change with the number of electrodes $({{\varepsilon }_{{\rm d}}}=0.043)$. In contrast, the linear inversion algorithm can utilize additional electrodes to shape the total current (${{\varepsilon }_{{\rm d}}}=0.043\to {{\varepsilon }_{{\rm d}}}=0.030$). This apparent difference between the conventional and linear inversion algorithms was also explored quantitatively; the parameters $\sigma$, $\rho$, and $\newcommand{\e}{{\rm e}} \eta$ were given three values covering the previous ranges, giving 27 parameter combinations. The value of ${{\varepsilon }_{{\rm d}32}}$ for the conventional and linear inversion algorithms was calculated for electrode numbers (${{n}_{{\rm e}}}$) between 2 and 32 (figures 8(a) and (b)).

Zoom In Zoom Out Reset image size

The accuracy of the conventional algorithm appears to improve with increasing number of electrodes until the stimulation overlap reaches a certain level and then the improvement ceases (figure 8(a)). In contrast, the linear inversion algorithm shows reducing improvement in ${{\varepsilon }_{{\rm d}32}}$ with increasing number of electrodes (figure 8(b)).

The ratio of the ${{\varepsilon }_{{\rm d}32}}$ for the conventional and linear inversion algorithms (figure 8(c)) demonstrates that, for γ  >  1, the linear inversion algorithm provides lower ${{\varepsilon }_{{\rm d}32}}$ in most circumstances. The exception is under the highest activation spread, $\sigma =1$, and highest measurement noise, $\newcommand{\e}{{\rm e}} \eta =0.3$, parameter combination.

Figure 9 shows the results of linear inversion under a range of noise and current-spread conditions. In the case with no noise, $\newcommand{\e}{{\rm e}} \eta =0$, the linear inversion algorithm produces a solution in which the electrodes take high amplitudes. For illustration purposes, we select a value of 2 as the maximum value, ${{s}_{\max }}$, twice the value of 1, which produces the maximum retinal activation of 1. So each element of $\vec {b}$ is set to be 2. With this selection it can be seen that under certain parameter combinations the electrodes exceed this maximum (box in figure 9). This is because wide activation spread and low noise conditions lead to very small eigenvalues being used in the pseudoinverse (figure 5(a), parameter combination 4). Use of a quadratic program with linear constraints on the electrode amplitudes can eliminate the high electrode amplitudes that arise using linear inversion seen in figure 9.

Zoom In Zoom Out Reset image size

Figure 10 compares the linear inversion and quadratic program solutions at both high noise, $\newcommand{\e}{{\rm e}} \eta =1$, and low noise, $\newcommand{\e}{{\rm e}} \eta ={{10}^{-5}}$, conditions.

Zoom In Zoom Out Reset image size

In the case of the linear inversion algorithm, low estimation noise and wide activation spread led to the selection of high-amplitude electrodes (figure 10(a)). However under the same conditions, the quadratic program avoided high electrode amplitudes while the value of $E$ increased only marginally, (${{\varepsilon }_{{\rm d}}}=0.012\to {{\varepsilon }_{{\rm d}}}=0.014$), (figure 10(b)). These results were examined quantitatively across the same range using the value ${{\varepsilon }_{{\rm d32}}}$ (figure 10(b)). It can be seen (dashed red line in figure 10(b)) that the quadratic program avoids the increase in electrode amplitude with reduced noise, $\newcommand{\e}{{\rm e}} \eta <0.05$. In both the qualitative (figure 10(a)) and quantitative (figure 10(c)) analyses it can be seen that increase in ${{\varepsilon }_{{\rm d32}}}$ that results from limiting the electrode amplitude is small. It can also be seen that outside these low-noise circumstances the quadratic program gives an identical outcome to the linear inversion algorithm by using the same restricted eigenvalue version of ${{{\bf W}}_{{\bf t}}}$, with ${{F}_{{\rm eig}}}=0.6$. The quadratic program avoids overfitting in the case of high noise wide activation spread conditions and also places an explicit limit on the amplitudes of the electrodes, avoiding overstimulation in low noise, wide activation spread conditions.

The limitation of the quadratic program approach is computational cost. In MATLAB, the process of calculating the electrode pattern for a single image using the quadratic program inverse model was ~1 s, which compares to the linear inversion time of ~ $1$ ms, and the conventional algorithm time of ~10 ms.

In calculating the linear inversion and the quadratic program solutions, the modulus operation in equation (5) was removed to allow for simple linear inversion equation (9). This decision limited those two algorithms to use only positive generator signal values to induce retinal activation (positive values of ${{G}_{i}}$ in figure 1(A)). Partitioning the target retinal activation pattern into positive and negative regions restores the modulus present in equation (5), allowing for retinal activation to be induced by either a positive or negative generator signal (figure 1(a)).

Figures 11(a)–(d) shows examples of target retinal activity patterns, $\vec {r}_{{\rm t}}^{*}$, potential beneficial partitions of these patterns, $\vec {p}_{{\rm t}}^{*}$, and the resulting retinal activation pattern using linear inversion, ${{\vec {r}}_{{\rm t}}}$. Borders between positive and negative regions of the generator signal pattern appear in the retinal activation pattern as closed loops of zero retinal activity. This is because the generator signal values must pass through zero between regions of positive and negative values, for example see figure 11(c).

Zoom In Zoom Out Reset image size

The fractures in the solutions derived from a partitioned target can lead to lower ${{\varepsilon }_{{\rm d}32}}$ in cases where the original target possesses such features (figures 11(a) and (b)). However, there are two ways in which partitioning could create ambiguities and artefactual features in the solution. In the first case, a partition assists to more accurately create an unclosed line of zero retinal activity in the target activation pattern, but the partition border must be a closed loop, and this leads to an artifact by creating a line of zero retinal activity where none exists (figure 11(c)). In the second case, the configuration of a closed loop in the target image may allow for a beneficial partition, but fails to do the same in an equivalent line (figure 11(d)). In target retinal activation patterns based on natural images, these two problems may not be encountered frequently. For example, in natural images, the algorithm may take advantage of ambiguities in the target patterns (figure 11(b)).

The binary quadratic program, solved using a Hopfield network with simulated annealing, was applied to the same 32 images throughout this investigation. Figures 12(a)–(e) compares a sample of the lowest and highest ${{\varepsilon }_{{\rm d}}}$ outcomes of using partitions to the results without a partition. Figure 12(e) shows the histogram of the difference in ${{\varepsilon }_{{\rm d}}}$ for each of the 32 images used.

Zoom In Zoom Out Reset image size

In 6 out of the 32 images the partition led to a detrimental outcome as measured by ${{\varepsilon }_{{\rm d}}}$. In 15 out of the 32 images the best outcome resulted from an unpartitioned target pattern. However improvements in ${{\varepsilon }_{{\rm d}}}$ were found 11 out of 32 images (figure 12). Examination of the most detrimental partitions, as measured by the value of ${{\varepsilon }_{{\rm d}}}$ (figures 12(a)–(e) bottom two rows), reveals that in qualitative terms, the result of partitioning is not necessarily worse than without.

The process of simulated annealing cycles of the Hopfield network for each image requires a large computational investment of ~100 s in MATLAB for each image.