The effect of stimulus level on excitation patterns of individual electrode contacts in cochlear implants

Objective: Spread of excitation (SOE) in cochlear implants (CI) is a measure linked to the speciﬁcity of the electrode-neuron interface. The SOE can be estimated objectively by electrically evoked compound action potential (eCAP) measurements, recorded with the forward-masking paradigm in CI recipients. The eCAP amplitude can be plotted as a function of the roving masker, resulting in a spatial forward masking (SFM) curve. The eCAP amplitudes presented in the SFM curves, however, reﬂect an interaction between a masker and probe stimulus, making the SFM curves less reliable for examining SOE effects at the level of individual electrode contacts. To counter this, our previously published deconvolution method estimates the SOE at the electrode level by deconvolving the SFM curves (Biesheuvel et al., 2016). The aim of this study was to investigate the effect of stimulus level on the SOE of individual electrode contacts by using SFM curves analyzed with our deconvolution method. Design: Following the deconvolution method, theoretical SFM curves were calculated by the convolution of parameterized excitation density proﬁles (EDP) attributable to masker and probe stimuli. These SFM curves were subsequently ﬁtted to SFM curves from CI recipients by iteratively adjusting the EDPs. We ﬁrst improved the EDP parameterization to account for stimulus-level effects and validated this updated parameterization by comparing the EDPs to simulated excitation density proﬁles (sEDP) from our com- putational model of the human cochlea. Secondly, we analyzed SFM curves recorded with varying probe stimulus level in 24 patients, all implanted with a HiFocus Mid-Scala electrode array. With the decon- volution method extended to account for stimulus level effects, the SFM curves measured with varying probe stimulus levels were converted into EDPs to elucidate the effects of stimulus level on the SOE. Results: The updated EDP parameterization was in good agreement with the sEDPs from the computational model. Using the extended deconvolution method, we found that higher stimulus levels caused signiﬁcant widening of EDPs ( p < 0.001). The stimulus level also affected the EDP amplitude ( p < 0.001) and the center of excitation ( p < 0.05). Concerning the raw SFM curves, an increase in current level led to higher SFM curve amplitudes ( p < 0.001), while the width of the SFM curves did not change signiﬁcantly ( p = 0.62). Conclusion: The extended deconvolution method enabled us to study the effect of stimulus level on excitation areas in an objective way, as the EDP parameterization was in good agreement with sEDPs from our computational model. The analysis of SFM curves provided new insights into the effect of the stimulus level on SOE. We found that the EDPs, and therefore the SOE, mainly became wider when the stimulus level increased. Lastly, the comparison of the EDP parameterization with simulations in our computation model provided new insights about the validity of the deconvolution method. (


Introduction
Modern cochlear implants (CI) are multi-channel devices with multiple electrode contacts on the electrode array. Taking advantage of the tonotopic organization of the cochlea, these different electrode contacts can provide different pitch percepts to the CI recipient ( Bonham and Litvak, 2008 ). In practice, however, it can https://doi.org/10.1016/j.heares.2022.108490 0378-5955/© 2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ) be difficult or even impossible to discriminate sounds produced by two adjacent electrode contacts ( Biesheuvel et al., 2019 ). A likely explanation for this phenomenon is the relatively large overlap between excitation areas of adjacent electrode contacts ( Bierer and Litvak, 2016 ). Previous research has shown that the overlap in excitation areas can be estimated objectively using electrically evoked compound action potentials (eCAP) generated by the auditory neurons ( Abbas et al., 1999 ). When measured via the forward-masking paradigm with fixed probe and roving masker stimulus contacts, the eCAP amplitude can be plotted as a function of the position of the masker contact, leading to the well-known spatial forwardmasking (SFM) curve ( Abbas et al., 1999 ;Cohen et al., 2003 ;Hughes and Abbas, 2006 ;van der Beek et al., 2012 ). However, the eCAP amplitudes presented in the SFM curve reflect an interaction between a masker and probe stimulus, making the SFM curve less reliable for examining excitation patterns or spread of excitation (SOE) attributable to individual electrode contacts. To overcome this problem, we developed a deconvolution method that translated SFM curves into excitation patterns attributable to the individual masker and probe stimuli ( Biesheuvel et al., 2016 ). The aim of this study was to investigate the effect of stimulus level on the SOE of individual electrode contacts.
Several studies have investigated the effects of current level on SFM curves. Cohen et al. (2004) reported that the widths of the SFM curves did not vary significantly with probe stimulus current, while the amplitudes of the curves varied systematically with stimulus level. Van der Beek et al. (2012) also did not find a statistically significant effect of stimulus level on the width of the SFM curves. Although some subjects showed a change in curve width with changing stimulus level, the majority of their subjects had SFM curves of similar widths at different stimulus levels. Looking at individual subjects, Hughes & Stille (2010) concluded that stimulus level has a significant effect on the width of the SFM curve. However, across all subjects, significant stimulus level effects were observed in only 34% of the cases. Lastly, Abbas et al. (2004) presented SFM curves with, in most subjects, a clear growth in eCAP amplitude and a tendency to become wider at higher stimulus levels. They suggested that the higher amplitudes were caused by more overlap in neurons excited by the masker and probe and that the increased width of the SFM curve could be partially attributed to an increased SOE. However, it is still unknown whether this increase in SOE is reflected by a higher excitation density (i.e., more excited nerve fibers locally) or by a wider spread along the electrode array, nor it is known whether it comes from the masker or the probe.
Most studies estimated SOE by determining the width of the SFM curve at a specific amplitude level ( Abbas et al., 2004 ;Cohen et al., 2003 ;Hughes and Abbas, 2006 ;Snel-Bongers et al., 2012 ;van der Beek et al., 2012 ). However, the width of a SFM curve is difficult to analyze if the curve is not complete, especially at the apical and basal ends of the electrode array ( Biesheuvel et al., 2016 ). Furthermore, a SFM curve reflects the interaction between a probe and masker stimulus, but the width of such a curve cannot reveal the individual contribution of each of the two stimuli. We developed a method to process SFM curves into so-called excitation density profiles (EDP) for individual electrode contacts using mathematical convolution ( Biesheuvel et al., 2016 ). An EDP, depicted in Fig. 1 A, reflects the proportion of excited neurons as a function of location along the electrode array. The deconvolution is performed using an iterative process. In the basis, parameterized EDPs (reflecting masker and probe excitation areas) are convolved to generate predicted SFM curves. Next, the predicted SFM curves are fitted to measured SFM curves (in CI recipients) by iteratively adjusting the EDPs. After many iterations, the method results in EDPs and SFM curves that predict the measured SFM curves quite well. A similar method has been presented by Cosentino et al. (2015Cosentino et al. ( , 2016 , who analyzed SFM curves using matrix algebra and Gaussian functions, also called the panoramic eCAP method. In our previous study, the EDPs were kept as simple as possible to clearly demonstrate the principle of (de)convolution ( Biesheuvel et al., 2016 ). When using the deconvolution method to study stimulus-level effects on SOE, we have to verify whether our EDP parameterization is suitable for that purpose. As in most parameterizations, defining the optimal number of EDP parameters is a challenging task; too few parameters could lead to a suboptimal representation of the excitation area, while too many parameters likely result in a good fit between the predicted and measured SFM curves but produce non-physiological EDPs. We have considered which EDP parameterizations were relevant for studying stimulus-level effects. First of all, we think that the EDPs must be able to get smaller and narrower with lower stimulus levels. This property has already been realized in the original parameterization, where the EDP width can be adjusted at the top of the Gaussian and the EDP amplitude is linked to the width parameter ( Biesheuvel et al., 2016 ). Based on the literature, we think there are two other useful parameterizations: (1) allowing variation in the center of excitation, and (2) allowing a gradual change in the slopes of the excitation pattern. The rationale for allowing some variation in the center of excitation is that studies have found an effect of stimulation level on perceived pitch ( Arnoldner et al., 2006 ;Carlyon et al., 2010 ;Shannon, 1983 ;Townshend et al., 1987 ), and that pitch might be determined by the centroid of the excitation area ( Laneau et al., 2004 ;McDermott and McKay, 2005 ;McKay et al., 1999 ). Kalkman et al. (2014) studied place pitch versus electrode location in a computational model of the implanted human cochlea; they found that the elicited pitch of electrode contacts after the first turn decreases as stimulus level increases. So, the above-mentioned studies indicate that an extra parameter that allows a shift in the center of the EDP might be useful. In a subsequent study, Kalkman et al. (2015) added spatial variability of the auditory nerve fibers' cell bodies to the model and examined neural excitation using excitation density curves. They found that the excitation density was lower at lower stimulus levels and spanned a shorter length along the cochlea. Their simulations further showed that the sides or edges of the excitation areas changed slightly with changing stimulus level: the density decay may be steeper or shallower. To include these effects in the EDP parameterization, we wanted to allow a gradual change in the sides of the EDPs.
The main goal of this study was to examine the effect of stimulus level on SOE by using the deconvolution method ( Biesheuvel et al., 2016 ). Our hypothesis was that the SOE is narrower at lower stimulus levels. The rationale for this hypothesis was that excitation density simulations in our computational model showed narrower SOE at lower current level ( Kalkman et al., 2015 ). Furthermore, the cell bodies are approximately positioned in a row in Rosenthal's canal. In that case, we expect using the physics of current spread that less stimulation would lead to narrower excitation areas, followed by less excitation density. We needed two stages to study our hypothesis. In the first stage, we had to improve the EDP parameterization to account for stimulus level effects. In the absence of exact information on real neural excitation patterns, we used the output of our computational model to validate the EDP parameterization ( Kalkman et al., 2015( Kalkman et al., , 2014. In this model, we can study the effect of stimulus level with the simulated excitation density, which strongly resembles the EDP from the deconvolution method (conceptually both describe the proportion of excited neurons along the electrode array). Henceforth the simulated excitation densities from the computational model were indicated with sEDP. In the second stage of this study, we measured SFM curves in human subjects and applied the deconvolution method to these SFM curves. For each sub- ject, two sets of SFM curves were collected. The first set consisted of SFM curves measured at electrode contact 3 through 16 with both masker and probe stimulus levels at 1200 μA, henceforth indicated as SFM 1200 . The second set consists of SFM curves measured at electrode contact 3 (apex), 9 (middle) and 15 (base), with a masker level of 1200 μA and probe levels ranging from 600 to 1100 μA, indicated as SFM 60 0-110 0 . Subsequently, the SFM curves of each subject were all together deconvolved into EDPs, revealing the SOE at individual electrode contacts and at different current levels of that subject.

Parameterization of EDPs
Essentially, an EDP is a Gaussian function extended with other parameters to reflect excitation density. The original deconvolution method included 19 variables to estimate the EDPs across the electrode array of a subject; 16 variables for the EDP plateau width and amplitude at each electrode contact, one variable for the EDP slope at all contacts (the slope depicts the density decay represented by the σ of the Gaussian function), one scaling variable to scale all predicted SFM curves to the measured SFM curves, and one variable as offset for the predicted SFM curves representing the noise floor (see also Fig. 1 ) ( Biesheuvel et al., 2016 ).
For studying stimulus level effects, the parameterization of the EDP width, the scaling factor and the offset parameter remained unchanged. The EDP slope was implemented slightly differently and a shift in the center of the EDPs (represented by the μ of the Gaussian function) was added. The shift and slope parametrization accounting for stimulus level effects was only applied on contact 3, 9 and 15, for which SFM 60 0-110 0 curves were collected as well. The slope and shift parameters were estimated at probe stimulus levels of 600 μA and 1200 μA and the intermediate stimulus levels were estimated using linear interpolation. Interpolation between the most extreme values 600 μA and 1200 μA is most robust. The interpolation limited the number of additional EDP parameters, and we did not expect large differences in slope and shift of EDPs obtained with succeeding stimulus levels at fixed electrode sites.
In more detail, we defined two EDP slope parameters: (1) slope 1200 , which estimated the slope at 1200 μA for all EDPs across the electrode array, and (2) slope 600 , which estimated the slope of the EDPs at 600 μA. The fitting range of both slope parameters was [0.1-6] and it had no unit. Subsequently, the EDP slope at the current levels between 600 μA and 1200 μA was estimated using linear interpolation. Regarding the relative shift in the center of excitation caused by the stimulus level, we assumed that the shift occurs at the higher stimulus levels because of deeper excitation in the spiral ganglion, while there is no shift at the lower current levels. The shift was represented by two parameters: (1) threshold shift , which estimated the current level where the shift was zero, and (2) shift 1200 , which estimated the shift of the EDP at the 1200 μA level. The fitting range of threshold shift was [0-600] μA and the range of shift 1200 was [-1-+ 1] electrode contacts. Subsequently, the shift as a function of stimulus level was defined as a linear relationship between the points (threshold shift μA; 0 electrode contacts) and (1200 μA; shift 1200 electrode contacts). To avoid unrealistic shift in the EDPs across the electrode array, for example, a + 1 shift at contact 9 and no shift at contact 10, the shift in the EDPs across the electrode array was estimated using a polynomial function. Given that the degree of the polynomial is limited by the number of data points, we could fit first-and second-order polynomials using the three shift 1200 values at contacts 3, 9, and 15. An overview of the final EDP parameterizations is shown in Table 1 .

Computational model of the human cochlea
The updated EDP parameterization was validated using excitation density simulations in our computational model of the implanted human cochlea, developed at the Leiden University Medical centre ( Briaire and Frijns, 20 0 0a,b, Frijns et al., 1995, 20 0 0, Kalkman et al., 2014, 2015, 2022. The model consists of two parts: a volume conduction model and a deterministic active nerve fiber model. The volume conduction model uses the boundary element method to simulate electrical potential distributions in a realistic three-dimensional geometry of a human cochlea implanted with a CI electrode array. Next, electrical potentials are determined along nerve fiber trajectories that have been defined in the cochlear geometry with the help of histological data. The nerve fiber model then simulates neural responses resulting from the electrical po-  (3,9,15), a total of 22 + 3 × 8 = 46 parameters had to be optimized. tential fields generated by the volume conduction model, by modeling each auditory neuron as an active double-cable electrical network using the human-based Schwarz-Reid-Bostock neural kinetics scheme ( Kalkman et al., 2022 ;Schwarz et al., 1995 ).
In humans, there is a high variability in neural responses. The actual eCAP-data to be measured depend on both the electric potential distribution and the neural activity. The electric potential distribution in the scala tympany is quite broad but uniform across subjects, while the neural activity is highly variable across subjects, depending on the neural survival and neural interactions ( Tang et al., 2011 ). To create a realistic validation model, we simulated excitation density profiles, denoted as sEDPs, at different probe current levels using five different cochlear geometries each modeled with three different states of neural health. Each cochlear geometry contained a representation of a HiFocus Mid-Scala electrode array with sixteen electrode contacts in ideal midscalar position. Three sets of 3200 auditory neurons were modeled for each geometry, each representing a different state of neural health: one set of intact nerve fibers, one set of fibers where the lengths of the unmyelinated peripheral terminal nodes were shortened from 10 μm to 1 μm, and one set of fibers from which the peripheral processes were completely removed. The peripheral processes of the intact neurons were spread out evenly along the basilar membrane (BM) and their cell bodies were spatially distributed along Rosenthal's canal in the manner described by Kalkman et al. (2015) . The stimuli used in the neural simulations were anodic-first biphasic pulses with 32 μs phase duration. At electrode contact 3, 9 and 15, the stimulus amplitude ranged from 600 μA to 1200 μA in steps of 100 μA, and at the other contacts across the electrode array the stimulus amplitude was 1200 h μA. The excitation density results from the neural model were processed into sEDPs, which show the excitation density as a function of cochlear position along the BM. Excitation density was defined as the percentage of excited neurons along 1 mm segments of the BM, centered at the positions of the tips of the peripheral processes of each modeled neuron. For the degenerated nerve fibers, excitation densities were determined as if the neurons still retained their peripheral processes. Further details of the computational model are discussed by Kalkman et al. (2015Kalkman et al. ( , 2014.

Analysis
In order to validate the EDP parameterization, we analyzed how well the EDPs could fit the sEDPs using the updated parameterization. To increase the sEDP variability and to strengthen the validation, the sEDPs were simulated in five different cochleae, each modeled with three neural states. The sEDPs at electrode contact 1, 2 and 16 were excluded from the analysis, because these sEDPs at the distal electrode contacts were often irregular (contact 1 and 2) or incomplete (contact 16) and fitting the EDPs to these sEDPs would lead to large fitting errors. The irregularities in sEDPs at contact 1 and 2 are likely caused by cross-turn stimulation in the apex and the sEDPs at contact 16 were often incomplete, as their simulation was hindered by the boundaries of the cochlea model, which basally ends at the round window. Moreover, the sEDPs at contact 1, 2 and 16 were strictly not necessary for validating the EDP parameterization to account for stimulus level effects on contact 3, 9 and 15. On the remaining electrode contacts (3 to 15), still 43 sEDPs (9%) were incomplete, e.g., at contact 14 or 15, and 43 sEDPs (9%) at lower stimulus levels did not show any excitation. These sEDPs did not lead to a unique solution and they were also excluded in the validation. We fitted EDPs to the remaining 379 sEDPs and calculated the similarity between each EDP and sEDP using the Jaccard index, which was defined as: The symbols ∩ and ∪ represent, respectively, the intersection and union of the sEDP and EDP. We compared the sEDPs and EDPs with respect to their slope (density decay at the sides of the EDP), full width at half maximum (FWHM), and the shift in excitation center using Bland-Altman plots and mountain plots ( Giavarina, 2015 ). In a Bland-Altman plot, the difference between two variables is plotted as function of their average. A mountain plot is a complementary representation of the difference plot. It shows the distribution of the differences by computing a percentile (p) for each ranked difference. The difference is plotted as function of p while p < 50, and otherwise as 100-p. Based on the outcomes from the Jaccard index, the Bland-Altman plots, and the mountain plots, we improved the EDP parameterization step by step, resulting in the parameterization as described in Section 2.1.1 .

Results
Fig. 2 shows an example of sEDPs from the computational model using different current levels (gray lines). The sEDPs originated from electrode contact 9 in cochlea model 4 with completely degenerated dendrites. When we fitted the EDP as parametrized conform our 2016 paper to these sEDPs, it turned out that there was a clear mismatch between the sEDPs and EDPs (panel A). This showed that an EDP parameterization update was necessary. After optimizing the EDP parametrization, we calculated the similarity between the sEDPs and fitted EDPs using Eq. (1) . Across all cochlea models, the similarity between the sEDPs and fitted EDPs increased from an 84.0% median (10th percentile: 69.1%, 90th percentile: 95.0%) using the 2016 parameterization, to a 92.2% median (10th percentile: 82.3%, 90th percentile: 96.5%) when using the new parameterization. Fig. 2 B shows that with the new parameterization the EDPs fitted much better to the sEDPs than in Fig. 2 A.
In Fig. 3 , we evaluated the updated EDP parametrization by comparing the EDPs with the sEDPs. We fitted EDPs to the sEDPs and explored the differences with respect to the slope (i.e., density decay at the sides of the EDP), full width at half maximum (FWHM), and the shift in excitation center. Each data point represents a single EDP comparison modelled in a specific cochlea and neural condition. Note that the differences on the y-axes were calculated as sEDP minus EDP and that data presented in this figure are the final results after optimizing the EDP parameterization. None of the differences plotted in panel A, C and E came from a normal distribution, as revealed by the Anderson-Darling test for normal distribution. Therefore, percentiles rather than the regular means and 95% confidence bounds are shown in the Bland-Altman plots. To gain more insight, the EDP 1200 data (blue) and EDP 60 0-110 0 data (red) were plotted separately. Panel A shows that the slope of the EDP 60 0-110 0 patterns differed from EDP 1200 ones, which ar- Fig. 2. Example of sEDPs (gray) and EDPs (black), both showing the excitation density along the basilar membrane. Panel A shows that the original EDPs, which were parameterized conform our 2016 paper, did not optimally fit to the sEDPs obtained at different current levels. Panel B shows that the updated EDPs, which account for stimulus level effects, fitted the sEDPs much better. EDP indicates excitation density profile; sEDP, simulated excitation density profile.
gued for allowing the slope to change towards lower current levels. Panel C confirms the importance of a width parameter, as it shows that the width of both the EDP 60 0-110 0 and EDP 1200 patterns varies highly across the different cochleae and neural conditions. Panel E shows that the center of excitation of the sEDPs varied between -0.47 and 0.23 electrode contacts, indicating that the parameterization of the shift was useful as well. The mountain plots at the right side of Fig. 3 show that for most EDPs the difference from the sEDPs was small; the peak of the mountain plot is steep around 0. Note that we rotated the mountain plot 90 °clockwise to align the y-axes of the Bland-Altman and mountain plots.
Based on the aforementioned outcomes, we assessed that the new EDP parameterization was suitable for the purpose of answering the aim of the study.

Patients and data collection
The patient data consist of intra-operative SOE recordings from 24 CI recipients, collected in the period from April 2015 to March 2017. The demographics of these subjects are shown in Table 2 . All subjects were implanted with a HiRes90K device with a Hi-Focus Mid-Scala electrode array from Advanced Bionics (Valencia, CA, USA). The electrode array contained 16 electrode contacts, numbered from apex (1) to base (16). The SFM curves were collected using the Bionic Ear Data Collection System (BEDCS) research software from Advanced Bionics, controlled by a custommade MATLAB (Mathworks, Inc., Natick, MA) interface. The neural responses were measured using monopolar biphasic pulses (anodic first) with a duration of 32 μs per phase. The sampling rate of the eCAP recording system was 56 kHz and its gain 100. The SFM curves were recorded using the forward-masking paradigm with a fixed probe contact and roving masker ( Biesheuvel et al., 2016 ;Cohen et al., 2003 ;Hughes, 2013 ). The masker-probe stimulus interval was 500 μs and the recording contact was two electrodes away from the probe contact in the apical direction. To save time, the probe and signature frames were recorded once for each SFM curve, and they were reused in all eCAP measurements of that SFM curve. This was possible because the probe and signature frames are independent from the masker location. To reduce the amount of noise in the recordings, the reused frames were measured with 64 averages and the masker and masker-probe frames with 16 averages ( Klop et al., 2009 ). The eCAPs were processed as described by Biesheuvel et al., (2016) , followed by a visual check to avoid erroneous processing of abnormal eCAPs (double peaks, shallow P 2 , etc.). In case of wrong peak detections, the peaks were corrected manually. For each patient, two sets of SFM curves were collected. The first set consisted of SFM curves measured at probe contacts 3 through 16 with both masker and probe stimulus levels at 1200 μA. Note that no SFM 1200 curves were measured for probe contacts 1 and 2 due to the requirement for an apically located record- Fig. 3. Quantitative analysis showing how well the updated EDP parameterization could fit sEDPs simulated at different current levels in five different cochleae and three different neural states. The panels at the left side represents Bland-Altman plots, in which the difference between an sEDPs and EDP parameter (y-axis) was plotted as a function of its average (x-axis). Each dot represents the result of fitting an EDP to a sEDP, whereby the EDP 1200 data (blue) and the EDP 60 0-110 0 data (red) were plotted separately. The slope (panel A) is calculated using the density decay between 40 and 60% of the maximum EDP amplitude. The black solid line shows the median and the dashed lines show the 2.5 and 97.5 percentiles. At the right side, mountain plots are plotted, wherein each ranked difference (y-axis) is plotted as function of percentile, or as 100-percentile if percentile > 50 (x-axis). EDP indicates excitation density profile; sEDP, simulated excitation density profile; e, electrode contact. ing contact. The second set consisted of SFM curves measured with a masker level of 1200 μA and different probe current levels. For most subjects, the probe current levels ranged from 600 to 1100 μA with a step size of 100 μA. For subjects S0125, S0126, and S0129, the probe current levels ranged from 600 to 10 0 0 μA with a step size of 200 μA because of time limitation during surgery. The SFM 60 0-110 0 curves were recorded at three probe contacts along the array: 3 (apex), 9 (middle), and 15 (base). In more than half of the patients, we could not measure the SFM 60 0-110 0 curves on all three probe contacts (see Table 2 ). The reasons were that sometimes there were time limitations during surgery, and sometimes the amplitudes of the SFM 1200 curves were already so small ( < 0.2 mV) that reliable measurement of SFM 60 0-110 0 curves was not possible.

Deconvolution of SFM curves
The SFM curves were analyzed according to the deconvolution method of Biesheuvel et al., (2016) , using a custom MATLAB application. Note that, if the maximum amplitude from the SFM curves decreased below 0.15 mV, these SFM curves were not included for further analysis, since the low eCAP amplitudes and bad signal-tonoise ratio of these curves caused unreliable EDP estimations. Due to this, 44 SFM curves (15%) were not included for further analysis. For each subject, the SFM 1200 and SFM 60 0-110 0 curves were all deconvolved at the same time in one minimization routine. The SFM 1200 curves were deconvolved into EDP 1200 patterns, reflecting the excitation areas of a 1200 μA stimulus at each electrode contact. Next, using the EDP 1200 estimations as masker EDPs, the SFM 60 0-110 0 curves were deconvolved into the probe EDP 60 0-110 0 patterns.

Analysis
Linear mixed models (LMMs) were used to test the effects of current level and electrode contact on the SFM 60 0-110 0 curves and EDP 60 0-110 0 estimates. The variables of interest were the maximum amplitude of the SFM curves, the width of the SFM curves (i.e., FWHM), the amplitude of the EDPs, the width of the EDPs (FWHM) and the shift in the center of the EDPs. Note that, if a side of the SFM curve did not drop to the FWHM-level, the width of the curve was set to the limit of the array in the apical direction (contact 3), or in the basal direction (contact 15) ( Abbas et al., 2004 ). The effect of current level and electrode contact on each of these variables was tested using separate models, including both current level and electrode contact as fixed effects, while subject was included as random effect. In addition, we repeated the EDP analyses on sEDPs as well in order to compare the human-based results with the output of the computational model. The output from the linear mixed models is displayed in Table 3 .

Results
An example of the deconvolution of SFM 60 0-110 0 curves is shown in Fig. 4 . The data were recorded at electrodes 3, 9, and 15 in subject S0119. The SFM curves in the upper half of the figure illustrate that an increase in stimulus level led to an increase in Table 3 Output from the linear mixed model, analyzing the effect of stimulus level on the amplitude, width and shift of the excitation patterns. SFM indicates spatial forward masking; EDP, excitation density profile; sEDP, simulated excitation density profile.

Electrode Current
Amplitude SFM IO curve p = 6.47 × 10 −13 * * p = 3.97 × 10 −50 * * Amplitude EDP IO p = 0.20 p = 2.54 × 10 −11 * * Amplitude sEDP IO p = 1.98 × 10 −7 * * p = 9.37 × 10 −26 * * Width SFM IO curve p = 4.17 × 10 −20 * * p = 0.62 Width EDP IO p = 5.33 × 10 −13 * * p = 2.58 × 10 −24 * * Width sEDP IO p = 5.91 × 10 −5 * * p = 2.74 × 10 −70 * * Shift EDP IO p = 0.14 p = 3.80 × 10 −2 * Shift sEDP IO p = 0.66 p = 0.11 * p < 0.05;. * * p < 0.001. amplitude as well. No major change in the width of the SFM curves is visible as a result of changing the stimulus level. Deconvolution of the SFM curves resulted in the EDPs plotted directly underneath them, in the lower row of the figure. The area under all EDP curves became greater when the stimulus level grew, which implies that lower stimulus levels excited fewer neurons. However, the EDPs at electrode 3 mainly increased in amplitude, while the EDPs at electrode 9 and 15 became wider. Fig. 5 shows the effect of current level on all EDP 60 0-110 0 series in more detail for all subjects. We plotted the area under the EDP, the amplitude, the FWHM, and the shift in center of excitation of each EDP as a function of stimulus level. In all cases, the area under the EDP became bigger with increasing current level ( Fig. 5 A-C), which is consistent with more neurons being excited at higher current levels. The smaller EDPs at lower current levels are mainly caused by narrower EDPs ( Fig. 5 G-I, increasing line towards higher current level) in combination with a still high EDP Fig. 4. Example of SFM curves recorded at different probe stimulus levels on electrode contacts 3, 9 and 15 (respectively panel A, B and C) in subject S0119. The error bar represents the eCAP amplitude with measurement error and the color shade represents the descending probe current level. The deconvolution of these SFM curves resulted in the EDPs directly shown below the curves (panel D-F). SFM indicates spatial forward masking; eCAP, electrically evoked compound action potential; EDP, excitation density profile. . Furthermore, the effect of current level on the shift in center of excitation was variable across subjects and electrode contacts ( Fig. 5 J-L). The LMM analyses revealed that, in general, the current level had a significant effect on the EDP amplitude ( p < 0.001), on the EDP width ( p < 0.001), and on the shift in EDP center ( p < 0.05). More details from the LMMs can be found in Table 3 . Effect of probe stimulus level on the SFM curves (blue), the EDPs (red) and the sEDPs (green). The parameters of interest (y-axis) were plotted as function of probe current level (x-axis). The line shows the average and the vertical bars indicate the standard deviation. The columns show the results for electrode contact 3, 9 and 15 separately. SFM indicates spatial forward masking; EDP, excitation density profile; sEDP, simulated excitation density profile; e, electrode contact.
We also compared the EDP 60 0-110 0 outcomes with the effect of current level on the SFM curves and the sEDPs from the computational model. Fig. 6 shows how the amplitude, FWHM and shift (if applicable) of the EDP 60 0-110 0 , sEDP 60 0-110 0 and SFM 60 0-110 0 curves changed as function of stimulus level. Observe that the sEDP amplitudes in Fig. 6 panel A and C show an irregular pattern; they increased from 600 to 700 μA, decreased until 800 μA and then increased again. This effect is caused by some missing data at 600 and 700 μA, as some sEDPs did not show any excitation at the lowest stimulus levels. So, the remaining data at 600 and 700 μA came from model simulations with a relatively high excitation density at these levels. The LMM analyses showed that current level had a significant effect on the amplitude of the SFM curve ( p < 0.001), while the FWHM of the SFM curve did not change significantly ( p = 0.62). Concerning the sEDPs, the analysis shows that current level had a significant effect on the amplitude ( p < 0.001) and the width ( p < 0.001), but not on the shift ( p = 0.11). See Table 3 for more details from the linear mixed-model analyses, including the effect of the electrode contact on the SOE.

Discussion
In this study, we investigated how excitation density profiles of cochlear implants change as a result of varying the stimulus level. Our hypothesis was that lower stimulus levels cause narrower excitation areas or at least excite fewer neurons. We studied this effect using our newly extended method for deconvolving SFM curves into so-called EDPs, which reflect the neural excitation by individual electrode contacts. The results of the present study support our hypothesis: lower stimulus levels excite fewer neurons, mainly represented by narrower EDPs.
In our previous study ( Biesheuvel et al., 2016 ), the parameterization of the EDPs was kept simple to clearly demonstrate the principle of (de)convolution. As we have explained in the introduction, the EDP parameterization had to be extended to account for stimulus level effects. Initially we implemented an EDP parameter for 'slope' and 'shift' per stimulus level. However, it turned out that this parameterization caused the EDPs to become physiologically unrealistic. For example, it was possible that the EDP at 1200 μA was centered around the stimulating electrode contact, while the EDPs at lower stimulus levels (60 0-110 0 μA) were all shifted by a full electrode contact; such a large and sudden shift cannot be explained physiologically. Therefore, we implemented a shift parameter that was linearly interpolated and must be zero somewhere at a lower stimulus level ( < 600 μA). In that way, we were able to study relative shifts caused by the stimulus level. Consequently, we could not account for another physiological phenomena that there might be an absolute shift at the lower stimulus levels depending on the position of the electrode contact in the cochlea. However, we think that an absolute shift depending on the position of the electrode array in the cochlea might be negligible, as all subjects had a HiFocus Mid-Scala electrode array which can be seen as a fixed factor.
Regarding the slope, it could happen that the EDP was shaped like a Gaussian curve at 1200 μA, while the EDPs at all other lower stimulus levels were almost rectangular. From a mathematical point of view, these rectangular EDPs are a valid solution in the deconvolution method. However, rectangular EDPs are physically and physiologically implausible ( Biesheuvel et al., 2016 ). Defining the EDP parameterization turned out to be a challenging task and, in order to verify whether the EDP parameterization was realistic, we compared the EDP parameterization with sEDPs from our computational model of the human cochlea ( Kalkman et al., 2022( Kalkman et al., , 2015( Kalkman et al., , 2014. Using the Bland-Altman analysis ( Fig. 3 ) and Jaccard index, we have improved the EDP parameterization stepwise and, finally, we came to the parameterization as described in Section 2.1.1 . Based on the results shown in Fig. 3 and the similarity median of more than 92%, we think that the final EDP parameterization presented in this study is good for the purpose of answering the central question of this study. Note that the EDPs are robust enough to predict excitation areas simulated in five anatomically different cochleae, including three states of neural health (i.e., intact nerve fibers, neurons with a shortened peripheral terminal node, and neurons that suffered a complete loss of their peripheral processes). Therefore, we conclude that the extended deconvolution method can now be used to study the effect of stimulus level on excitation areas.
Before we go into detail about the effect of stimulus level on SOE, we want to discuss another aspect regarding the validity of the deconvolution method. Thus far, the deconvolution method was based on eCAP data collected with anodicleading stimuli, while both the literature and current CI systems often use cathodic-leading stimuli. The deconvolution method, however, was developed using anodic-first stimuli and we continued on this in the current study ( Biesheuvel et al., 2016 ). We think that the used anodic-leading stimuli did not have major consequences for the outcomes in the current study; we performed within-subject-analyses with anodic-leading data only. Further, research from Hughes et al., (2016) showed no difference in eCAP amplitude measured with either anodicfirst or cathodic-first forward masking in AB recipients and our method is based on these eCAP amplitudes. Nevertheless, the excitation patterns underlying the anodic-first and cathodic-first stimuli may differ. A preliminary analysis comparing the sEDPs with sEDPs simulated under equal conditions but with cathodic-first stimuli, showed a similarity of 81.6% (according to Eq. (1) ). Visual inspection revealed that the differences between cathodic-first and anodic-first sEDPs were mainly in the form of changes in the widths and amplitudes of the excitation patterns. Both the width and amplitude are parameterized in our EDPs and the deconvolution method. So, we think that our method can deal with the different stimulus polarities correctly. In a follow-up study, the stimulus polarity can be an area of interest, especially if we want to compare the EDPs with behavioral data collected with cathodicleading stimuli.
Regarding the stimulus level effects on SOE, our hypothesis was that a lower stimulus level will excite fewer neurons, due especially to a narrowing of the excitation areas. The results presented in this study confirm this hypothesis ( Fig. 5 ). Both the EDPs and sEDPs became significantly wider when the stimulus level increased ( p < 0.001). A few cases showed that the EDPs did not necessarily become wider with increasing stimulus level. Higher stimulus level could initially also lead to higher EDP amplitudes, as shown in Fig. 4 D, reflecting deeper excitation in the spiral ganglion. In general, the stimulus level was always related to the area under the EDP, which reflects the number of excited neurons. This is consistent with previous findings that the area under the EDP is highly correlated to the eCAP amplitude, given that eCAP amplitudes depend on current level ( Biesheuvel et al., 2016 ). Fig. 6 D-F shows that the sEDPs are narrower, or more selective, than the EDPs found in patients. Unfortunately, we do not yet have a definitive explanation for this outcome. It is possible that neural excitation patterns in the computational model are too selective compared with SOE data collected in humans. In that case, the results of the present study invite us to critically review the computational model, especially the parameters that are involved in simulating SOE. Meanwhile, it is important to realize that, although EDPs and sEDPs are conceptually very similar, strictly speaking they are different. The sEDPs are derived from simulated excitation patterns; the modeled nerve fibers are characterized by their positions along the BM and each excited neuron contributes to the sEDPs equally. Additionally, the widths of the sEDPs are affected by the length of the segment along the BM over which the excitation density is averaged. This length was set to 1 mm for this study, but this is essentially an arbitrary value. In contrast, EDPs are obtained from mathematical deconvolution of SFM curves recorded in humans and they are a function of electrode contact spacing, rather than distance along the BM. In short, although the EDP and the sEDPs are both intended to represent neural excitation densities, their differences should be kept in mind when comparing the two to each other directly. Regardless, for the purpose of this study, which is validation of the EDP parameterization, the fact that the sEDPs appears to be more selective does not matter, since the width of the EDP is a variable in the deconvolution method as well.
The updated EDP parameterization allows the center of excitation to shift as a function of stimulus level. The bottom row of Fig. 5 shows how the center of the EDPs shifted with the current level in our patient data. Note that for all three electrode contacts a number of EDPs shifted to the maximum of the fitting range, plus or minus a full contact at 1200 μA. Although it is undesirable for the minimization routine to reach the minimum or maximum of the fitting range, extending the fitting range would be unrealistic. It turned out that when the shift parameter was not constrained, some EDPs tended to shift 4 or 5 electrode contacts. It is not likely that an EDP would shift so far (even more than one electrode contact) when the current level is increased from 600 μA to 1200 μA. We think we encountered a limitation of the deconvolution method that is likely caused by the incomplete SFM curves at distal electrode contacts. At electrode contacts 3 and 15, the SFM curves are asymmetric and most of the masker contacts involved in the deconvolution are unilateral from the probe contact. When the probe stimulus level was lowered the eCAP amplitudes became smaller and we would expect smaller probe EDPs. However, in the deconvolution method, smaller eCAP amplitudes are coded by decreased overlap between the masker and probe EDPs. We have seen that this can be achieved by narrowing the EDPs, but it can also be represented in a mathematically correct way by shifting the probe EDPs away from the masker EDPs, especially for distal electrode contacts (see also Fig. 5 J-L). However, the shift parameter did not reach the boundaries of its fitting range in all sub-jects. So, we can also speculate that this problem of reaching the fitting range is caused by a within-subject factor, rather than by the edges of the SFM curves. For example, in the case of a dead neural region we would expect a smaller EDP, but a dead neural region can also be represented mathematically by shifting the adjacent EDPs away. Unfortunately, in the context of this study, we cannot rule out such dead-region effects. To estimate dead regions, we would need a different EDP parametrization followed by a validation step using an independent measure for neural status.
Dead neural regions, however, are interesting as they are directly related to the clinical applicability of the deconvolution method. If the deconvolution method can detect poor neural survival, such information can be used to optimize the cochlear implant fitting, e.g., by de-activating less effective electrode contacts. Thus far, our studies regarding the deconvolution method are quite fundamental. Nevertheless, the present results about stimulus level effects on excitation patterns can help audiologist with CI fitting, especially making them aware that higher stimulus levels may limit the spatial resolution in most subjects. In a next step, it would be valuable to work on the detection of poor neural regions. Garcia et al., (2021) recently published a study in which they investigated the assessment of neural health using the panoramic eCAP method. They concluded that the panoramic eCAP method can detect neural survival patterns with high accuracy (at least 90%). However, a disadvantage of their method is that they used computer simulations with a backward approach, so that both the generated and solved dead regions are based on the same model assumptions. Instead, it would be better to validate the detection of dead regions by an independent model, in order to increase the clinical feasibility. For instance, it would be possible to model poor neural regions with our computational model. Subsequently, based on these poor neural regions, eCAP-based SFM curves can be simulated, and these curves can be translated back to EDPs using the deconvolution method. Next, the obtained EDPs can be compared with the simulated dead regions to validate them.
The observed shift in the center of the EDPs might indicate a change in pitch percept, assuming that the center of excitation is involved in pitch perception by cochlear implant users ( Laneau et al., 2004 ;McDermott and McKay, 2005 ;McKay et al., 1999 ). Such a pitch shift with intensity is a well-known phenomenon in normal hearing ( Stevens, 1935 ;Terhardt, 1979 ), but also in line with several studies reporting an effect of stimulation level on perceived pitch in electrical hearing ( Arnoldner et al., 2006 ;Carlyon et al., 2010 ;Shannon, 1983 ;Townshend et al., 1987 ). Shannon (1983) asked a CI user to judge pitch as function of loudness level and this subject perceived an increase in pitch with increased loudness. Townshend et al., (1987) studied pitch perception in three implant subjects and they found that pitch percepts were primarily affected by the place and rate of stimulation, while the stimulus level appeared to be a secondary effector. Among these three patients, one perceived an increased pitch at higher stimulus levels, while the other two perceived a decreased pitch at higher levels. About 20 years later, other studies still show varying results. Arnoldner et al., (2006) showed an increase in pitch at higher current levels in 10 patients, whereas only one patient showed the opposite effect. Carlyon et al., (2010) studied the effect of stimulus level and place of stimulation on temporal pitch perception by cochlear implant users. They found that in 16 of 21 cases the pitch increased with signal level, while the opposite effect was observed in the other five cases. If we assume that pitch perception is linked to the excitation area, the results of the present study show similar outcomes, especially at electrode contacts 3 and 15 where the shift in the excitation center varied between -1 and + 1 ( Fig. 5 G-I). Carlyon et al., (2010) concluded that the stimulus-level effects observed in their study cannot entirely be assigned to changes in current spread with increasing level and they also speculated that temporal coding might play a role in pitch perception. So, to further elucidate the relationship between stimulus level and pitch perception, we suggest combining the eCAP-based EDPs with psychophysical measures of pitch perception.

Conclusions
The deconvolution method enables us to study the effect of stimulus level on excitation areas in an objective and physiological way. The validation step showed that the updated parameterization of the EDPs describes the excitation areas simulated in our computational model well. Using this updated EDP parameterization, eCAP-based SFM curves measured with different probe current levels can be deconvolved in EDPs, thereby unraveling the effect of stimulus level on the excitation areas. It turned out that excitation areas are smaller, mainly narrower, at lower stimulus levels, which could not be derived from the SFM curves included in this study. Lastly, the comparison of the EDP parameterization with simulations from our computational model provided new insights about the validity of the deconvolution method.

Funding
This study was supported by the Dutch Technology Foundation STW (project number 11693) and Advanced Bionics.