The relation between polarity sensitivity and neural degeneration in a computational model of cochlear implant stimulation

The main aim of this computational modelling study was to test the validity of the hypothesis that sensitivity to the polarity of cochlear implant stimulation can be interpreted as a measure of neural health. For this purpose, the effects of stimulus polarity on neural excitation patterns were investigated in a volume conduction model of the implanted human cochlea, which was coupled with a deterministic active nerve fibre model based on characteristics of human auditory neurons. The nerve fibres were modelled in three stages of neural degeneration: intact, with shortened peripheral terminal nodes and with complete loss of the peripheral processes. The model simulated neural responses to monophasic, biphasic, triphasic and pseudomonophasic pulses of both polarities. Polarity sensitivity was quantified as the so-called polarity effect (PE), which is defined as the dB difference between cathodic and anodic thresholds. Results showed that anodic pulses mostly excited the auditory neurons in their central axons, while cathodic stimuli generally excited neurons in their peripheral processes or near their cell bodies. As a consequence, cathodic thresholds were more affected by neural degeneration than anodic thresholds. Neural degeneration did not have a consistent effect on the modelled PE values, though there were notable effects of electrode contact insertion angle and distance from the modiolus. Furthermore, determining PE values using charge-balanced multiphasic pulses as approximations of monophasic stimuli produced different results than those obtained with true monophasic pulses, at a degree that depended on the specific pulse shape; in general, pulses with lower secondary phase amplitudes showed polarity sensitivities closer to those obtained with true monophasic pulses. The main conclusion of this study is that polarity sensitivity is not a reliable indicator of neural health; neural degeneration affects simulated polarity sensitivity, but its effect is not consistently related to the degree of degeneration. Polarity sensitivity is not simply a product of the state of the neurons, but also depends on spatial factors.


Introduction
In most clinical stimulation strategies, cochlear implants (CIs) use symmetric biphasic pulses to electrically excite the auditory neurons located in the cochlea ( Loizou, 2006 ). These symmetric biphasic pulses consist of a rectangular cathodic and anodic phase of equal duration and amplitude, which is the simplest way to construct a charge balanced electrical pulse that can be used safely, without causing harmful electrochemical reactions in the cochlea ( Brummer et al., 1983 ;Donaldson and Donaldson, 1986 ). However, nerve fibres in general are not equally responsive to anodic and cathodic stimuli, since the ions involved in the transmission of neural signals all have positive charge and are therefore not indifferent to the polarity of an external electric field ( Frankenhaeuser and Huxley, 1964 ;Colombo and Parkins, 1987 ). Indeed, psychophysical and electrophysiological studies in humans and animals have confirmed that there are differences in thresholds and dynamic range between the two polarities ( Miller et al., 1999 ;Macherey et al., 2006 ;van Wieringen et al., 2008Undurraga et al., 2010 ;Carlyon et al., 2013Carlyon et al., , 2013 Cazals, 2016 ;Carlyon et al., 2017 ); however, this sensitivity to stimulus polarity is not yet fully understood. Recently, several studies have hypothesized that polarity sensitivity can be interpreted as an indicator of neural health ( Mesnildrey, 2017 ;Carlyon et al., 2018 ;Hughes et al., 2018 ;Goehring et al., 2019 ;Jahn and Arenberg, 2019a ;Mesnildrey et al., 2020 ). If this is the case, then polarity sensitivity experiments would be a valuable tool for assessing the state of the auditory neurons in individual CI subjects. The aim of the present study is therefore to shed more light on this subject by investigating the effects of CI stimulus polarity and neural degeneration in a computational model of the human cochlea.
Objective measures and psychophysical experiments have either suggested or demonstrated that at high stimulus amplitudes, human subjects are more sensitive to anodic pulses than to cathodic ones ( Macherey et al., 2006van Wieringen et al., 2008 ;Macherey et al., 2010 ;Undurraga et al., 2010 ;Carlyon et al., 2013Carlyon et al., , 2013Carlyon et al., 2017 ;Hughes et al., 2018 ). Conversely, animal experiments have shown that, at high amplitude, animals such as cats and guinea pigs are more sensitive to cathodic pulses rather than anodic ones ( Miller et al., 1999 ;Macherey and Cazals, 2016 ), implying that polarity sensitivity is affected by species specific cochlear anatomy, neurophysiology, aetiology/duration of deafness or some combination of these. Additionally, computational modelling studies have suggested that different parts of the auditory neurons respond to different polarities; specifically, that cathodic pulses tend to stimulate the peripheral processes of the neurons and that anodic pulses are more likely to stimulate the central axons ( Rattay, 1999 ;Rattay et al., 2001aRattay et al., , 2001bPotrusil et al., 2020 ). Furthermore, computational modelling studies have also shown that the degree of peripheral process stimulation can vary depending on both the location of the stimulating contact along the cochlear duct as well as its distance from the modiolus ( Frijns et al., 1995 ;1996 ;2001 ;Hanekom, 2001 ;Briaire and Frijns, 2006 ;Smit et al., 2010 ;Kalkman et al., 2015 ;Potrusil et al., 2020 ), raising the question of how much the positions of the electrode contacts affect their polarity sensitivity.
Recent studies have investigated polarity sensitivity in human CI subjects at perceptual threshold and have found inconsistent results. Although anodic thresholds were usually lower than cathodic thresholds, there was a great amount of variability, both between subjects and within each subject ( Mesnildrey, 2017 ;Carlyon et al., 2018 ;Goehring et al., 2019 ;Jahn and Arenberg, 2019a ;Mesnildrey et al., 2020 ). The insights provided by computational modelling has prompted the researchers of these experiments to hypothesise that their inconsistent results could partially be explained by differences in neural health (i.e. the state of the surviving auditory neurons), specifically that greater sensitivity to cathodic stimuli indicates the presence of healthy peripheral processes. The underlying reasoning is that degeneration of the peripheral processes will affect the cathodic thresholds considerably, while their absence should not make much difference for the anodic thresholds, since those are more dependant on the central axons. This kind of retrograde degeneration, the gradual disappearance of peripheral processes, has indeed been observed in human auditory neurons after hearing loss ( Kujawa and Liberman, 2009 ;Lin et al., 2011 ;Wu et al., 2019 ), indirectly adding to the plausibility of the hypothesis.
Despite this, the evidence supporting the neural health hypothesis is inconclusive, as there are still unknown factors that could play a role. A potentially complicating factor is the fact that all experiments in humans must be performed with charged balanced multiphasic pulses; many of the cited studies use either pseudomonophasic pulses or triphasic pulses as charge-balanced approximations of monophasic pulses, since true monophasic stimuli are not safe for clinical use. This raises the additional question whether the use of these pulses has distorted the results or not, particularly for triphasic pulses, since their 'undesired' phases still have considerable amplitudes. Computational modelling can provide valuable insight in these matters, as they allow researchers to manipulate factors and variables that are unchangeable or unknown in human CI subjects and make it possible to simulate clinically unsafe stimuli such as monophasic pulses.
The main goal of the present study was to test the validity of the hypothesis that polarity sensitivity can be interpreted as an indicator of neural health in CI stimulation. For this purpose, the study used a computational model of the implanted human cochlea, which was developed over the years at the Leiden University Medical Centre ( Frijns et al., 1995( Frijns et al., , 1996Briaire and Frijns, 20 0 0a , 20 0 0b ;Frijns et al., 20 0 0 , 20 01 ;Briaire and Frijns, 20 05 ;20 06 ;Frijns et al., 20 09a , 20 09b , 2011 ;Westen et al., 2011 ;Snel-Bongers et al., 2013 ;Kalkman et al., 2014Kalkman et al., , 2015, updating its neural model with a human-based kinetics scheme and an accompanying auditory nerve fibre morphology, represented as an electrical double cable . This model was used to simulate neural responses to monophasic, biphasic, triphasic and pseudomonophasic pulses in five cochlear geometries implanted with lateral, mid-scalar and medial electrode arrays, under three different states of neural health (i.e. healthy intact nerve fibres, neurons with minor degeneration at the peripheral tip and fibres that were completely missing their peripheral processes).

Methods
The present study used an updated version of a previously published model of the implanted human cochlea, which consisted of a volume-conduction model for simulating electrical field distributions and a deterministic active nerve fibre model for simulating the resulting neural responses ( Kalkman et al., 2014 ;2015 ).

Model geometry
The volume-conduction model is based on the Boundary Element Method (BEM) and is used to calculate quasi-static electrical potential fields in a three-dimensional geometrical representation of a human cochlea implanted with an electrode array. The model contains five human cochlear geometries, labelled CM1 through CM5, the first two of which were derived from histological crosssections while CM3, CM4 and CM5 were based on μCT reconstructions. CM1 through CM4 are the same cochlear geometries used previously in Kalkman et al. (2015) ; CM5 is a not previously published geometry and is depicted in Fig. 1 .
Electrode arrays were modelled in a lateral, mid-scalar and medial (peri-modiolar) position in each of the five cochlear geometries, leading to a total of 15 model geometries. The modelled array dimensions were based on the HiFocus1J electrode (Advanced Bionics, Valencia, CA, USA) and have rectangular plate contacts with a longitudinal length of 0.5 mm, a transversal length of 0.4 mm, and an electrode spacing of 1.1 mm (centre to centre). The arrays were inserted as deeply as possible without touching the walls of the scala tympani and the number of modelled electrode contacts for each individual array was chosen so that the electrodes span the maximum possible cochlear range in all modelled cochleae. Table 1 shows the active angular ranges and number of contacts for each geometry used in this study.
The auditory nerve fibre trajectories in each cochlear geometry are modelled in the manner described in our previous studies ( Kalkman et al., 2014 ;2015 ). The model geometries all contain 3200 neurons, the peripheral tips of which are evenly distributed along the basilar membrane (BM) while the cell bodies are spatially distributed in Rosenthal's canal in pseudo-random fashion without breaking their tonotopic organisation (i.e. from base to apex the fibres retain the same order in Rosenthal's canal Fig. 1. Illustration of one of the cochlear geometries used in this study (CM5). Panel a shows a mid-modiolar plane through a human cochlea, reconstructed from μCT imaging of a human temporal bone, obtained from Advanced Bionics and the University of Antwerp. Black lines in panel a indicate boundaries between cochlear structures used to define the compartments of the model geometry; the green lines indicate the outlines of the neural trajectories. Panel b shows a ray-traced image of a cut-through of the final three-dimensional geometry, implanted with an electrode array in lateral position. Modelled intact neural trajectories are shown in yellow.

Table 1
Insertion angles of the most basal and apical electrode contacts of all model geometries used in the study (angles are measured in degrees from the round window), as well as the number of electrode contacts of each array geometry ('#').

Lateral
Mid - as they have along the BM). The modelled neurons are grouped into 80 bundles of 40 nerve fibres; these nerve bundles condense when they enter the Habenula Perforata, expand around the spiral ganglion, and condense again when they exit Rosenthal's canal. The trajectory of the central line of each nerve bundle is determined by applying data from a histological study by Stakhovskaya et al. (2007) , which found a relationship between the positions of the auditory neurons' peripheral tips along the organ of Corti and the positions of their cell bodies along the spiral ganglion.

Neural model
The auditory nerve fibre model from our previous studies, which was largely based on animal data, was found to insufficiently behave like human (auditory) neurons in experimental studies, particularly in terms of predicting action potential shape, latencies and conduction velocity ( Briaire and Frijns, 2005 ;Bachmaier et al., 2019 ). To remedy some of these shortcomings, a new auditory nerve fibre model was used which incorporated some human anatomical features and properties that were absent in the previous model. First, the new neural model described nerve fibres as an electrical double cable, rather than a single cable, in order to include peri-axonal current between the axon and the first layer of myelin ( Berthold, 1978 ;Berthold and Rydmark, 1983 ). In addition, the previously used generalised Schwarz-Eikhof-Frijns (gSEF) kinetics, which were derived from rat and cat experiments, were replaced with the so-called Schwarz-Reid-Bostock (SRB) neural kinetics scheme ( Schwarz et al., 1995 ), which was based on human data. This neural model was previously published using a homogenous nerve fibre morphology ; for this study the morphology was adapted to match the human auditory neuron (illustrated in Fig. 2 ), incorporating anatomical details from histological observations ( Spoendlin and Schrott, 1989 ).
As shown in Fig. 2 a, the modelled nerve fibres consisted of a thinly myelinated cell body flanked by two 1 μm long Ranvier nodes, connecting it to a peripheral process and a central axon. The central axon consisted of 22 myelinated segments separated by Ranvier nodes; the lengths of these segments increased from 150 μm to 350 μm in steps of 50 μm, in central direction starting from the cell body. The peripheral process of the base fibre illustrated in Fig. 2 consisted of 6 myelinated segments of equal length, likewise separated by Ranvier nodes, with a 10 μm long unmyelinated terminal at the tip. Since the peripheral processes vary in length throughout the cochlea, the lengths of the myelinated peripheral segments were adjusted accordingly. However, an upper limit was placed on the peripheral segment lengths as it was found that longer segments would block the propagation of action potentials. Therefore, in extremely long peripheral processes (particularly those in the apex) segments were added to ensure that their lengths would not exceed 280 μm, equivalent to the length limitation that was necessary for the previous version of the neural model ( Kalkman et al., 2014 ). As a result, the lengths of the peripheral segments in the cochlear geometries ranged from 131 μm to 280 μm, with the number of segments ranging from 6 to 11.
In order to simulate the effects of neural degeneration, the neurons were modelled in three conditions: completely intact (as described above and depicted as the top fibre in Fig. 2 a), with shortened peripheral terminals and with complete loss of the peripheral processes. For the short terminal condition the 10 μm unmyelinated nodes at the tips of the peripheral processes of all modelled neurons were shortened to 1 μm (middle fibre in Fig. 2 a), in the same manner as in a previous study ( Snel-Bongers et al., 2013 ), where it was shown that this seemingly minor change in neural morphology can have a noteworthy impact on neural activation. This condition was based on animal studies that showed that the initial stage of auditory neural degeneration is loss of the peripheral terminal ( Kujawa and Liberman, 2009 ;Lin et al., 2011 ). For the third neural condition, the entire peripheral process was removed from all modelled neurons, leaving only a 1 μm long Ranvier node on the peripheral side of the cell body (bottom fibre in Fig. 2 a). This condition represents a state of severe neural degeneration where it is nevertheless still possible to excite the auditory nerve fibres electrically.

Fig. 2.
Schematic representation of the model's intact auditory nerve fibre morphology (not to scale). The top fibre depicted in subfigure a illustrates how the base intact model fibre is built out of myelinated segments; the peripheral process consists of six myelinated segments (left side) and has an axon diameter of 2 μm. The peripheral segments are all the same length (L p ), which varies throughout the cochlea, ranging from 131 μm to 280 μm. It is followed by a cell body of 30 μm length and 10 μm diameter, flanked by two nodes of Ranvier. The central axon (right side) consists of 23 segments and has an axon diameter of 3 μm. The first myelinated segment after the cell body is 150 μm long and each segment after that if 50 μm longer than the previous, up to a maximum of 350 μm. All nodes of Ranvier are 1 μm long, except the terminal node of the peripheral process (left end of the figure), which is 10 μm long for an intact neuron. The middle fibre in subfigure a represents a partially degenerated neuron, in which the peripheral terminal node has been reduced to a length of 1 μm; the bottom fibre shows a more severe state of degeneration where the entire peripheral process has been removed. Subfigure b gives a more detailed illustration of a myelinated segment flanked by nodes of Ranvier (NR), showing how the segment is split into three subsegments of equal length ( λ). The myelin sheath on each segment consists of a number of layers: 20 for the peripheral process, 1 for the cell body and 70 for the central axon. Furthermore, there is a peri-axonal layer located between the axon and the myelin sheath, which is 4 nm wide. Subfigure c shows the modelled electrical diagram of a node of Ranvier (left part) and the first myelinated subsegment. The following symbols are used: C m for membrane capacitance; P Na for voltage dependant sodium permeability; g Kf and g Ks for voltage dependant conductance for the fast and slow potassium currents; G L for membrane leak conductance; V L for leakage equilibrium potential; G a , G p , G e1 and G e,2 for axonal, peri-axonal and extracellular conductance; C my for myelin capacitance; G my for myelin conductance.

Stimulus configuration
With the implanted cochlear geometries and the model neurons defined, the volume-conduction model was used to calculate the electrical potentials at each Ranvier node and myelinated internode (see Fig. 2 b) resulting from injecting a monopolar 1 μA current at one of the modelled electrode contacts. Since the volume-conduction model was purely resistive, the potentials could simply be multiplied by an electrode stimulus pulse to obtain time-dependant nerve fibre potentials, which could then be used as external field for the neural simulations. To study the effects of stimulus polarity on neural excitation, simulations were run with four types of pulse shape: monophasic, biphasic, triphasic and pseudomonophasic pulses, as illustrated in Fig. 3 . The monophasic pulses were simulated to isolate polarity effects without the complication of opposite polarity phases interacting with each other; these were compared to biphasic pulses (the standard pulse shape in clinical settings), as well as triphasic and pseudomonophasic pulses (clinically safe pulse shapes that are frequently used as approximations of monophasic pulses in clinical experiments ( van Wierin-gen et al., 2005 ;Macherey et al., 20 06 , 20 08 ;van Wieringen et al., 2008 ;Macherey et al., 2010 ;Carlyon et al., 2013 ;Cazals, 2016 , 2017 ;Mesnildrey, 2017 ;Carlyon et al., 2018 ;Goehring et al., 2019 ;Arenberg, 2019a , 2019b ;Mesnildrey et al., 2020 )).
The simulated monophasic pulses consisted of a single 40 μs rectangular pulse with either anodic polarity (MA; Fig. 3 a1) or cathodic polarity (MC; Fig. 3 a2). The biphasic pulses used in this study consisted of two phases of 40 μs duration with no interphase gap; the two phases were rectangular pulses with equal amplitude but opposite polarity. Biphasic pulses lead with either the anodic phase (BA; Fig. 3 b1) or with the cathodic phase (BC; Fig. 3

b2).
Triphasic pulses consisted of a central 40 μs phase flanked by two opposite polarity phases with the same duration but half amplitude and no interphase gaps. These triphasic pulses were named by the polarity of their central phase, which was expected to be its dominant phase due to its higher amplitude; Fig. 3 c1 shows the triphasic anodic pulse (TA) used in this study, while Fig. 3 c2 shows the triphasic cathodic pulse (TC). Pseudomonophasic pulses were essentially modified biphasic pulses; the amplitude of the second phase of a biphasic pulse was reduced with an arbitrary Fig. 3. Illustration of the different pulse shape definitions: monophasic (a1&a2), symmetric biphasic (b1&b2), triphasic (c1&c2) and pseudomonophasic pulses (d1&d2). The top row of panels (a1-d1) represent the pulses where the anodic phase is expected to be dominant (or simply the first phase, in the case of biphasic pulses), while the bottom row (a2-d2) represents their cathodic counterparts. The main pulse duration for all pulses in the study is 40 μs, for pseudomonophasic pulses this is the duration of the first phase (the duration of the second phase is variable, but it is 200 μs long in the example in panels d1&d2). See the methods section for further details. ratio R ps , while the duration of the second phase was increased with the same ratio to maintain charge balance (i.e. with the first phase set at 40 μs, the second phase was R ps •40 μs long). For a sufficiently large value of R ps , this pulse was expected to behave like a monophasic pulse, with the benefit of being charge balanced and therefore usable in clinical practice. An example of a pseudomonophasic anodic pulse (PSA) with R ps = 5 is shown in Fig. 3 d1, and the corresponding pseudomonophasic cathodic pulse (PSC) is shown in Fig. 3 d2. Three different values of R ps were used for the simulations: R ps = 2 (with the resulting pseudomonophasic pulses referred to as PSA2 and PSC2), R ps = 4 (PSA4 and PSC4), and R ps = 8 (PSA8 and PSC8).

Model output
Using the pulse definitions above as monopolar stimuli on each electrode contact in all 15 geometries, simulated neural responses were calculated for all modelled nerve fibres under the three described conditions of neural degeneration. For each simulation it was recorded whether an action potential propagated to the central end of the given nerve fibre, and if so, at which node the action potential originated. Stimulus amplitudes (i.e. the amplitude of the main phase of each pulse type) ranged from 0.05 mA to 5 mA in steps of 0.404 dB; once the calculations for this amplitude range were completed, additional iterative simulations were done to determine the threshold of excitation of each individual fibre more precisely, to within 0.1% of its precise value. Of these additional iterations only the lowest above-threshold value for the stimulus amplitude and the corresponding node of excitation were recorded.
For each individual electrode contact and stimulus pulse, simulation results were compiled into so-called excitation profiles, which were represented as two-dimensional colour maps that indicate which model neurons were excited at a given stimulus amplitude, with different colours indicating in which part of the fibre the excitation took place: peripheral process, cell body or central axon (excitation in one of the two Ranvier nodes on either side of the cell body was considered to be excitation at the cell body ( Frijns et al., 1995( Frijns et al., , 1996 ; Briaire and Frijns, 20 0 0a )).
As in our previous studies, simulated loudness levels were quantified by the amount of space the excited neurons occupied along the BM. Each modelled nerve fibre was considered to occupy a certain amount of length along the BM at their peripheral tips (for this purpose, degenerated nerve fibres were treated as if they were still intact) and thus the simulated loudness level of a given stimulus was determined by summing the lengths along the BM for all excited neurons. Based on comparisons of our model to psychophysical data, perceptual threshold was considered to be equivalent to 1 mm excitation along the BM in the model ( Snel-Bongers et al., 2013 ), while maximum comfortable loudness (MCL) was considered to be comparable to an excitation width of 4 mm; the stimulus amplitude needed to reach perceptual threshold and MCL in the model will therefore be referred to as I 1mm and I 4mm , respectively.
In order to quantify polarity sensitivity, several studies have defined the so-called polarity effect (PE) as the difference between the perceptual thresholds of cathodic and anodic stimuli ( Mesnildrey, 2017 ;Carlyon et al., 2018 ;Goehring et al., 2019 ;Jahn and Arenberg, 2019a ;Mesnildrey et al., 2020 ), using the following equation: Here, PE is expressed in dB, I cathodic is the current amplitude at threshold for the cathodic stimulus, and I anodic is the current amplitude at threshold for the anodic stimulus.  Fig. 4 are typical examples that illustrate how stimulus polarity and neural degeneration determined the location of excitation along the auditory nerve fibres. In general, MA pulses excited neurons almost exclusively along the central axon, which can be seen in Fig. 4 a, where all excitation is coded blue. This was the case regardless of array type, electrode contact or degree of neural degeneration (data not shown), though these factors did affect the I 1mm and I 4mm levels. It should be noted, however, that neural degeneration had only a minor impact on simulated perceptual loudness levels for MA stimulation.

Excitation patterns
By contrast, in simulations with MC pulses, neurons near the stimulating electrode contact were rarely excited along the central axon. For intact and short terminal neurons, excitation occurred mostly along the peripheral processes (green area in Fig. 4 b and c), while a complete loss of peripheral processes moved the place of excitation to the cell bodies (red areas in Fig. 4 d). When MC pulses excited neurons at the central axon, they generally did so at high stimulus levels (above I 4mm ) at sites deep in the modiolus, in neurons associated with higher/lower cochlear turns (blue areas in Figs. 4 b-d). This kind of unintended excitation was reported in previous studies, where it was referred to as cross-turn stimulation ( Frijns et al., 1995 ;Briaire and Frijns, 20 0 0a , 20 01 ;Kalkman et al., 2014 ). Notably, the threshold for cross-turn stimulation is considerably lower for MC stimulation than it is for MA stimuli (compare blue areas in the top right corners of Figs. 4 b-d to the same area in Fig. 4 a).
Additionally, the dependency on the peripheral processes also meant that I 1mm and I 4mm levels of MC stimuli were much more affected by neural degeneration than those of MA pulses. One notable observation from the excitation profiles is that, for lateral contacts, the I 1mm threshold was generally lower for short terminal neurons than it was for intact neurons, due to the presence of a spatially restricted area of excitation at the peripheral tips of the neurons closest to the stimulating contact (visible as a sharp peak left of the dashed lines at 18 mm along the BM in Fig. 4 c). This region of excitation was much more present in short terminal neurons and lead to a lower I 1mm than was found for intact neurons (compare Fig. 4 b to 4 c).

Simulated loudness levels and pe values in different cochlear geometries
The top row of panels in Fig. 5 shows I 1mm levels and resulting PE values of MA and MC pulses, stimulated on all lateral array contacts in each cochlear geometry, modelled with intact neurons.  ( Verbist et al., 2010 ), which is measured in degrees from the round window. The data points connected by dot-   Fig. 5 (panels a2-c2) repeats the same plots, but using I 4mm levels, rather than I 1mm levels (the PE values in panel c2 are therefore calculated at simulated MCL, rather than threshold).
Simulated thresholds, MCLs and PE values of individual geometries generally followed the same trend as the model average, with some individual variations (compare dotted curves to the solid black curves in Fig. 5 ). The model average PE in panel c1 showed a nearly monotonically decreasing curve, going from 2.3 dB at the basal end of the cochlea (which means that MA thresholds were 2.3 dB lower than MC thresholds), to −1.5 dB at 495 °( which means MC thresholds were 1.5 dB lower than MA thresholds), while the curves for the individual geometries showed erratic variations to the average trend ( Fig. 5 c). Most notably, although the average PE curve only crosses the 0 dB line once, individual PE curves can cross it multiple times (blue circles, green squares and purple triangles for CM1, CM3 and CM5 respectively in Fig. 5 c), which indicates that the favoured polarity of these geometries changed several times along the lengths of their arrays. PE values calculated using I 4mm levels were similar to those calculated at I 1mm , though PE values were shifted downwards by about 1 dB between roughly 90 °and 360 °insertion angle.

Comparing degenerated neurons to healthy neurons
The effects of neural degeneration on polarity sensitivity are demonstrated in Fig. 6 . PE values were determined using monophasic stimuli on all electrode contacts of each electrode array and cochlear geometry, for all three neural conditions and were plotted against each other. In Fig. 6 a, PE values of individual contacts stimulating short terminal neurons are plotted along the ordinate, while the abscissa indicates the corresponding PE values for the intact neurons. Fig. 6 b shows the same plot, but with PE values of neurons with completely degenerated peripheral processes along the ordinate. In both plots the data points are split up by array type (blue circles for lateral contacts, red triangles for midscalar contacts and green squares for medial contacts) and the grey diagonal lines indicate the line where the PE values for the degenerated neurons were equal to those of intact neurons. The error bars along the axes indicate the means and standard deviations of the PE values for the corresponding neural condition, split up by array type (the black error bar represents all three array types combined). Asterisks next to the error bars denote statistically significant differences between data sets ( p < 0.05).
In both plots of Fig. 6 it is apparent that the modelled changes in neural condition did not have a consistent effect on PE values, as the data is spread out on both sides of the diagonal, which means that neural degeneration increased the PE for some contacts (data points above the diagonals) but decreased it for others (data points below the diagonals). However, there was a clear effect of array type, reflected by the way the data points for each type of array are clustered together in Fig. 6 . This clustering was mainly due to the fact that PE values for intact neurons were strongly affected by electrode array type (compare the error bars along the abscissae), while mean PE values for degenerated neurons were less affected (compare the error bars along the ordinates). For lateral contacts, neural degeneration tended to shift the PE downward (blue circles in Figs. 6 a and b), as opposed to PE values of medial electrodes, which were more likely to be shifted upwards (green squares). For mid-scalar electrode contacts, shortening the terminal node of the modelled neurons had a relatively minor im- . Each point in the plots represents an individual electrode contact from one of the model geometries, the data has been split up by array type (blue circles for lateral arrays, red triangles for mid-scalar arrays and green squares for medial arrays). Error bars along the axes indicate the means and standard deviations of the PE values, with the black error bars representing the total data set and the other bars corresponding to the data sets indicates by their colours and symbols. Black asterisks next to the error bars indicate statistical significance between the data sets.
pact on the PE, as the corresponding data points are mostly located close to the diagonal, but on average slightly shifted upwards (red triangles in Fig. 6 a), whereas completely removing the peripheral processes mostly shifted PE values downwards (red triangles in Fig. 6 b). Paired t-tests revealed that these changes in the PE values between neural conditions were all statistically significant for each individual array type; when the data for all three array types were merged, paired t-tests still showed significant differences between the neurons without peripheral processes and the other two neural conditions, but not between the intact and short terminal neurons.

Comparing multiphasic pulses to monophasic pulses
Next, the effect of pulse shape on the calculated PE values at I 1mm was examined; the results are plotted in Fig. 7 . The figure shows PE values obtained with simulated perceptual thresholds of multiphasic pulses on the ordinate, plotted against the corresponding PE values of monophasic pulses on the abscissa. Data from five different cathodic/anodic multiphasic pulse pairs are plotted: three implementations of pseudomonophasic pulses, namely PSC8/PSA8 (Pseudomonophasic 8; blue circles), PSC4/PSA4 (Pseudomonophasic 4; red upward pointing triangles), PSC2/PSA2 (Pseudomonophasic 2; green squares), biphasic pulses BC/BA (yellow diamonds) and triphasic pulses TC/TA (purple downward pointing triangles). Each data point represents an individual electrode contact in one of the 15 geometries, under one of the three neural conditions, and the grey diagonal represents the line where the PE values of the multiphasic pulses were equal to the PE values obtained with monophasic pulses. PE values obtained with pseudomonophasic 8 pulses were very similar to those found with true monophasic pulses, as indicated by the fact that their data points (blue circles) follow the grey diagonal closely in Fig. 7 . Inspection of simulated threshold levels show that, on average, I 1mm of PSC8 and PSA8 pulses were slightly shifted up from their monophasic counterparts by about 0.4 dB and 0.3 dB respectively, resulting in PE values that were on av- Fig. 7. Effect of pulse shape on simulated polarity sensitivity of each electrode contact in all model geometries (cochleae and array types) and all three neural conditions. The abscissa indicates PE values obtained with monophasic pulses, while the ordinate indicates the corresponding PE values derived from simulations with multiphasic pulses, namely pseudomonophasic pulses with R ps equal to 8, 4 and 2 (blue circles, red upward pointing triangles and green squares, respectively), symmetric biphasic pulses (yellow diamonds) and triphasic pulses (purple downward pointing triangles). The grey diagonal represents the line where simulations with monophasic pulses result in PE values that are identical to their multiphasic equivalents. erage 0.1 dB higher than monophasic PE values. At the other extreme end, PE values of biphasic pulses deviated from the diagonal substantially, with the data points arranged almost horizon- contacts in all model geometries, plotted separately for each neural condition; panels a1 and a2 show data from simulations with intact neurons, panels b1 and b2 show data from short terminal neurons and panels c1 and c2 contain results from neurons with missing peripheral processes. The top row of panels (a1-c1) shows the relationship between normalised PE values and normalised average thresholds, plotted analogously to the psychophysical data in Carlyon et al. (2018) . Simulated perceptual thresholds are obtained with anodic and cathodic triphasic pulses; the average values of these two are normalised and plotted on the abscissa, while the resulting PE values are also normalised and plotted on the ordinate. The data are normalised for each combination of cochlear geometry, array type and neural condition, by calculating the average value across the entire array and subtracting the result from that of each individual electrode contact. In each panel, blue circles represent lateral electrode contacts, red triangles represent mid-scalar contacts and green squares correspond to medial contacts; linear regression fits are plotted as black lines (panel a1: r 2 = 0.26, p << 0.001, panel b1: r 2 = 0.09, p << 0.001, panel c1: r 2 = 0.01, p = 0.09). The bottom row of panels (a2-c2) shows the same data as in the top row, but now plotted against electrode insertion angle (measured in degrees from the round window) along the abscissa. Data in each panel is again split by array type and black lines are the linear regression fits (panel a2: r 2 = 0.62, panel b2: r 2 = 0.32, panel c2: r 2 = 0.50; p << 0.001 for each regression line). tally in Fig. 7 (yellow diamonds), indicating that there was a much smaller difference between I 1mm values of BC and BA pulses than there was between those of MC and MA pulses. The average PE value obtained with biphasic pulses was −0.7 dB, showing that the model tended to be slightly more sensitive to BC pulses than to BA pulses. Since a pseudomonophasic pulse with R ps = 1 is identical to a biphasic pulse, the plots for pseudomonophasic 4 and pseudomonophasic 2 pulses in Fig. 7 can be seen as a stepwise transition from pseudomonophasic 8 pulses to biphasic pulses, with the data points going from arranged close to the diagonal to being arranged along a decreasing slope.
With triphasic pulses, PE values were consistently higher than those obtained with monophasic pulses, as evidenced by the fact that the data points for triphasic pulses are all located above the diagonal in Fig. 7 (purple downward pointing triangles). This shift was on average 1.2 dB and was mostly due to an increase in cathodic thresholds; I 1mm levels for TC pulses were on average 2.2 dB higher than those of MC pulses, while I 1mm levels of TA pulses were only 1 dB higher than those of MA pulses.

Normalised simulated loudness levels and pe values
To remove the variance between subjects, recent studies have normalised their data by subtracting average values for each subject ( Mesnildrey, 2017 ;Carlyon et al., 2018 ;Mesnildrey et al., 2020 ). The same was done for the modelling results of this study in Fig. 8 , where normalised PE values of triphasic pulses are plotted against their normalised average thresholds (panels a1-c1), analogous to the manner the psychophysical data was plotted in Carlyon et al. (2018) . For this figure, every combination of cochlear geometry, array type and neural condition was treated as a 'virtual subject'; the data points in the plot correspond to individual electrode contacts of these subjects. The data was normalised by calculating the average values across the entire array and subtracting them from the values of each individual electrode contact for each virtual subject. The ordinates show normalised PE values, while the abscissae indicate the average of the I 1mm levels of TA and TC pulses, with I 1mm expressed in dB relative to 1 mA. In each panel of Fig. 8 the data is split by array type: blue circles for lat-eral arrays, red triangles for mid-scalar arrays and green squares for medial arrays. The three panels illustrate the effects of different neural conditions; the data for intact neurons is plotted in panels a1 and a2, the results of the short peripheral terminals are shown in panels b1 and b2, and panels c1 and c2 corresponds to complete loss of peripheral processes. The black line in each panel shows a linear regression fit of the plotted data points. Fig. 8 shows significant positive correlations between PE and average threshold level for intact neurons and short terminal neurons (panels a1 and b1), as was apparent from the linear regression lines (panel a1: r 2 = 0.26, p << 0.001, panel b1: r 2 = 0.09, p << 0.001). When there is a complete loss of peripheral processes, the correlation was negative but non-significant (panel c1: r 2 = 0.01, p = 0.09). Splitting the data by array type revealed no significant differences between their regression slopes, but there was an effect of insertion angle. To illustrate this, the normalised PE values have been replotted in panels a2-c2, but here the abscissa indicates the insertion angles of the electrode contacts, rather than normalised average I 1mm levels. For each neural condition there was a negative correlation between normalised PE value and insertion angle, as indicated by the linear regression lines (panel a2: r 2 = 0.62, panel b2: r 2 = 0.32, panel c2: r 2 = 0.50; p << 0.001 for all three regression lines).
To expand on this analysis a multivariate regression was performed with the normalised PE values as response variables and with the predictor variables consisting of the normalised average I 1mm threshold levels of TA and TC pulses, the insertion angles, the interaction between those two variables, and an intercept term (i.e. PE ∼ 1 + Angle + Mean_I 1mm + Angle:Mean_I 1mm , in Wilkinson-Rogers notation ( Wilkinson and Rogers, 1973 ), as it is generally implemented in software packages such as MATLAB and R). Performing this multivariate regression separately for each neural condition reiterated that normalised PE was significantly negatively correlated with insertion angle ( p << 0.001 for each neural condition). However, in contrast to the results of the simple linear regressions of Fig. 8 , correlations with normalised average threshold levels were weak and non-significant for intact and short terminal neurons ( p = 0.6 and p = 0.4, respectively), while there was a significant negative correlation for neurons without peripheral processes ( p << 0.001). Significant but weak positive correlations were found between normalised PE and the interaction between normalised average threshold and insertion angle, but only for intact fibres and neurons without peripheral processes ( p = 0.03 and p << 0.001, respectively); for short terminal neurons the correlation was negative but not significant ( p = 0.1). The r 2 values of the multivariate regression fits were 0.64 for the intact neurons, 0.32 for the short terminal neurons and 0.69 for the neurons without peripheral processes.

Discussion
This computational modelling study was performed to test the hypothesis that polarity sensitivity of CI stimulation can be seen as an indicator of neural health. Neural responses to anodic and cathodic stimuli were simulated in a realistic model of a human cochlea implanted with an electrode array, which consisted of a volume conduction model and an active nerve fibre model. The latter was updated to better simulate the human auditory neurons by employing a new morphology and using the Schwarz-Reid-Bostock neural kinetics scheme ( Schwarz et al., 1995 ), both of which were based on human data. Modelled stimuli consisted of monophasic, biphasic, triphasic and pseudomonophasic pulses, which were applied to lateral, mid-scalar and medial electrode contacts in five different representations of human cochlea. Neural responses were simulated for three different conditions of neural health: intact, with shortened peripheral terminal nodes and with missing peripheral processes.
The modelling of neural health was based on observations of gradual retrograde degeneration in histological studies ( Kujawa and Liberman, 2009 ;Lin et al., 2011 ;Wu et al., 2019 ). To limit the scope of this study, only uniform degeneration (or lack thereof) was modelled and no attempt was made to implement more lifelike neural conditions, where the extent and type of degeneration varies along the cochlear duct. Furthermore, truly realistic representations of neural degeneration would not only be limited to a shortening of the peripheral processes and would include different types of degeneration or neural dysfunction that were not accounted for in this study. For instance, there are reports of demyelination of auditory neurons after hearing loss in animal studies ( Tagoe et al., 2014 ;Wan and Corfas, 2017 ). A recent modelling study looked at the effects of such axonal myelin loss and concluded that it caused changes in the polarity sensitivity of the auditory neurons ( Resnick et al., 2018 ), which is pertinent to the present study. As another example, certain (genetic) diseases or cochlear trauma may affect the density or functioning of sodium channels in the auditory neurons ( Fryatt et al., 2011 ;Heffner et al., 2019 ;Meisler et al., 2021 ). As a preliminary investigation, we attempted to gauge the effect of sodium channel degeneration/dysfunction by repeating a subset of our simulations with the voltage dependant sodium permeability parameter P Na reduced by 50%. We found that this raised PE values for neurons with (mostly) intact peripheral processes, but that for neurons without peripheral processes the PE values were largely unaffected (data not shown); in other words, the polarity sensitivity changes due to reducing P Na interacted with the effects of degeneration of the peripheral processes. This illustrates how complex it would be to rigorously model neural health, since the different types of neural degeneration/dysfunction discussed here are not mutually exclusive and any combination could have some interaction effect. So, while they were relevant to the present study, alternative types of neural degeneration were ignored here in favour of simplicity. The way neural health was represented in this study should therefore not be seen as comprehensive or lifelike, but rather as an exploration of specific extreme cases to gain insight into how polarity sensitivity can change.

Polarity dependency of the site of excitation
Model results showed that anodic pulses predominantly excited the auditory neurons in their central axons, while cathodic stimuli tended to excite them in their peripheral processes or near their cell bodies. This was in line with previous modelling work in literature ( Rattay, 1999 ;Rattay et al., 2001aRattay et al., , 2001bPotrusil et al., 2020 ); as in those studies, the polarity-dependant site of excitation was determined by the way the model geometries shaped the electrical potential fields and by the trajectories of the modelled neurons in them. Site of excitation was therefore not an innate property of the modelled neurons, but mainly a result of spatial factors, including electrode position (i.e. both the distance from the modiolus as well as the cochlear insertion angle of the stimulating contact). Despite this, changes in lateral to medial electrode positioning did not meaningfully affect the 'targeted areas' of the two polarities, though they did affect anodic and cathodic thresholds, as well as the resulting PE.
The tendency of the two polarities to excite the auditory neurons in different locations also meant that, in general, thresholds of cathodic stimuli were more affected by neural degeneration than anodic thresholds were. As a result, changes in PE values due to neural degeneration were mostly the result of changes in cathodic thresholds. However, the effects of degeneration were inconsistent at perceptual threshold and sometimes counter intuitive. For ex-ample, simulations showed that shortening the terminal node on the peripheral processes of the auditory neurons or removing the peripheral processes altogether often lowered simulated thresholds, rather than raising them as one would expect. The explanation for this is that the terminal node and the peripheral process itself both have a capacitance that must be overcome in order to generate an action potential at the terminal node or near the cell body, respectively. Shortening the terminal node or removing the peripheral process reduces that capacitance, which makes it easier to excite the peripheral end of the neuron (or the remainder thereof). However, this effect and its magnitude depended not only on stimulus polarity, but also on the stimulating contact's distance from the modiolus and its insertion angle.

Using multiphasic approximations of monophasic pulses
When comparing the model's results to clinical data, one should keep in mind that the model predicts that multiphasic pulses designed to approximate monophasic pulses will not generally behave exactly as true monophasic pulses. The resulting difference in obtained PE values will depend on the multiphasic pulse used; pseudomonophasic pulses with long, shallow second phases will act similarly to monophasic pulses (in Fig. 7 , the simulated PE values using PSA8 and PSC8 pulses were on average only about 0.1 dB higher than those determined with MA and MC pulses), but the range of measured PE values will become increasingly narrow as the second phase becomes shorter and its amplitude higher. For triphasic pulses, PE values are consistently shifted upwards relative to monophasic pulses, by about 1.2 dB on average ( Fig. 7 ). This is important since most recent studies have used triphasic pulses in their polarity sensitivity experiments and even a 1.2 dB shift can have a considerable effect on the number of electrode contacts that have positive or negative PE values.
At perceptual threshold (I 1mm ), the model's average polarity sensitivity to monophasic stimuli is mostly in agreement with psychophysical data. The absolute values of model thresholds are much higher than thresholds of actual CI subjects, but this is a known problem of the type of model used in this study . Despite this, the average range of simulated PE values is similar to those found in human test subjects ( Mesnildrey, 2017 ;Goehring et al., 2019 ;Jahn and Arenberg, 2019a ;Mesnildrey et al., 2020 ), though the model leans more towards negative PE values, whereas the literature reports mostly positive PE values. However, model results are brought more in line with the psychophysical data if one considers that these studies all used triphasic pulses. For individual geometries, model simulations showed erratic changes in PE values along the electrode array ( Fig. 5 ), comparable to psychophysical data ( Mesnildrey, 2017 ;Carlyon et al., 2018 ;Goehring et al., 2019 ;Jahn and Arenberg, 2019a ;Mesnildrey et al., 2020 ). Notably though, these kinds of erratic patterns are present in simulations with completely intact neurons and idealised electrode array positioning, so they arise mostly from inhomogeneities in cochlear anatomy and neural trajectories.

Polarity sensitivity at high stimulus levels
Unfortunately, simulated PE values at MCL (I 4mm ) are much less consistent with data from literature. While findings in human subjects are inconsistent, most studies report that their human test subjects were more sensitive to high amplitude anodic stimuli than cathodic ones ( Macherey et al., 20 06 , 20 08 ;van Wieringen et al., 2008 ;Macherey et al., 2010 ;Undurraga et al., 2010 ;Carlyon et al., 2013Carlyon et al., , 2013Carlyon et al., 2017 ;Hughes et al., 2018 ). However, the model more often shows the reverse, which cannot entirely be explained by the use of monophasic versus multiphasic stimuli, so there may be some shortcoming of either the model itself or the interpretation of its output. Both the volume conduction model and the neural model are based on well-known laws of physics, but there remains some uncertainty in the precise details and electrochemical properties of the cochlea and the auditory neurons, despite the available data in literature.
However, there is one important unknown involved in interpreting the model's data, namely the question of how loudness is coded; in both the present study as well as past ones, loudness in the model has been assumed to be linearly proportional to the number of excited neurons. This has previously led to model results that were comparable to psychophysical data, particularly when simulating loudness growth curves of multipolar current focussing strategies ( Kalkman et al., 2015 ). Nevertheless, this is a very simplistic way to model loudness, since it ignores temporal information such as spike rate and spike timing. Furthermore, in this study only single pulse stimuli have been simulated, while experiments with real life test subjects are usually performed with pulse trains, which means that for a more thorough comparison, stochastic neural effects should also be included in the model. If there were large enough differences in stochasticity between different parts of the neurons or polarity dependant effects, then that may explain the discrepancy in high stimulus amplitude polarity sensitivity between model and clinical data. Van Gendt et al. have made steps to expand the computational model to include stochasticity and enable pulse train simulations in a practical manner ( van Gendt et al., 2016 ;2017 ;, but its application lies outside the scope of the present study. To further expand on interpreting loudness in the model, the present study found no occurrences of the kind of non-monotonic loudness growth that has been reported in psychophysical studies with quadriphasic and triphasic pulses Mesnildrey, 2017 ). The authors of these studies speculated that their findings may have been the result of a 'cathodal block' that was observed in physiological experiments ( Ranck, 1975 ) and simulated in a previous modelling study ( Frijns et al., 1996 ). A cathodal block occurs when a stimulus depolarises one node, but (nearly) simultaneously hyperpolarises another node that is located more centrally along the nerve fibre, thereby preventing the action potential from propagating. Although cathodal blocks were also observed in the present study (data not shown), they mainly happened just below the threshold levels of individual nerve fibres. In other words, the presence of a cathodal block effectively raised the thresholds of some model neurons, and once stimulus amplitude was high enough to overcome the block, the action potential would propagate along the rest of the fibre. There were rare cases of cathodal blocks occurring at higher current levels, but only for isolated fibres and within small ranges of amplitudes, which meant that they did not meaningfully affect the simulated loudness growth curves. The model is therefore currently unable to confirm the cathodal block hypothesis posed by Macherey et al. and Mesnildrey as an explanation for their findings. This could indicate that the model is not behaving realistically enough yet, but it is also possible that the explanation for non-monotonic loudness growth must be found elsewhere, for example in temporal effects such as spike timing.

Differences in polarity sensitivity between animals and humans
As mentioned in the introduction, experiments in humans have revealed different polarity sensitivity than animal experiments. Although properly investigating these differences in our model would require repeating the simulations of this study in an equivalent animal model, some preliminary insight could be gained from the present study by running additional simulations with modified neural morphologies. Since the most notable difference between the auditory neurons of humans and experimental animals is the degree of myelination of the cell bodies, the number of myelin layers around the modified neurons' cell bodies was increased from 1 to 100 . In all other aspects these neurons were identical to the neurons used in the rest of this study, including being modelled in three stages of neural degeneration. The results of these simulations are compared to those of the default (human-like) model neurons in Fig. 9 , which shows PE values from the modified neurons plotted along the ordinate, against the corresponding PE values from the default model neurons along the abscissa.
Adding myelin around the cell bodies of the model neurons universally decreased the PE of all model geometries, electrode contacts and neural conditions, which is apparent from the fact that all data points in Fig. 9 are located below the diagonal. The degree to which the PE decreased was dependant on the presence of the peripheral processes; PE values from neurons without peripheral processes were lowered considerably more than those from intact neurons and neurons with shortened terminal nodes. Examining the I 1mm levels of the anodic and cathodic monophasic pulses individually revealed that adding myelin around the cell bodies increased the threshold levels of anodic pulses and decreased those of cathodic pulses; both effects therefore contributed to the decrease in PE values.
These results suggest that differences in the degree of myelination around the cell bodies of the auditory neurons are an important factor underlying the differences in polarity sensitivity between humans and animals. Similar observations were made recently by Potrusil et al. ( Potrusil et al., 2020 ), though they did not observe a change in anodic thresholds. Additionally, they only varied the amount of cell body myelination of neurons without peripheral processes, while results of the present study indicate that this effect occurs with healthy neurons and neurons with short-ened peripheral terminal nodes as well, though the changes in PE values were more severe when peripheral processes were missing entirely.
It should be reiterated that these are inconclusive findings, as there are significant differences between animals and humans that have not been accounted for (chiefly cochlear anatomy, precise neural morphology and neural kinetics), so they should be the subject of future research (cf. Frijns et al., 2001 ).

The neural health hypothesis
The results of this study cast doubt on the neural health hypothesis, which is exemplified by the fact that neural degeneration in the model had counter-intuitive and often inconsistent effects. The electrode array's distance from the modiolus greatly affected the direction in which the simulated PE changed when comparing different states of neural health (though not at higher stimulus levels). Furthermore, in the most extreme case of neural degeneration included in this study, the model average PE values were universally negative, which meant that cathodic thresholds were lower than those of anodic pulses. This is the opposite of what would be expected based on the neural health hypothesis, which assumes that loss of peripheral processes would make anodic stimulation more favourable and thus lead to positive PE values. The important observation here is that, for the average lateral and mid-scalar electrode arrays, severe neural degeneration tended to decrease PE rather than increase it, while for the average medial array neural degeneration did increase PE at threshold, but never to positive values. Although again, using triphasic stimuli instead of true monophasic pulses produced higher PE values (on average 1.2 dB) than the model averages shown in the error bars along the axes in Fig. 6 . This was a large enough shift to lead to positive PE values in a number of cases, particularly in the most basal half turn of the cochlea, so it provides at least a partial explanation for why model PE values tended to be lower than those found in psychophysical studies.
As stated above, erratic changes in PE values along electrode arrays could be seen in simulations involving completely intact neurons, meaning that similar patterns found in psychophysical studies may not necessarily arise from localised neural degeneration, also known as dead regions. Similarly, results from a psychophysical study by Carlyon et al. (2018) were replicated while simulating healthy neurons ( Fig. 8 a1). Carlyon et al. found a significant positive correlation between normalised PE values and normalised average threshold of anodic and cathodic stimuli, which, at first glance, seemed consistent with the idea that absence of the peripheral processes would increase the cathodic thresholds but not the anodic ones, which would increase their combined average as well as their resulting PE value. However, the present study complicates that interpretation, as it shows that this correlation can also arise from cochleae with completely healthy neurons. Moreover, Fig. 8 b1 and c1 show that in the model the correlation disappears when neural degeneration becomes more severe.
Plotting the same data from our model against insertion angle, rather than normalised average thresholds, revealed that the normalised PE was negatively correlated with the insertion angle of the stimulating contact, regardless of array type or neural degeneration ( Fig. 8 a2-c2). The same could also be seen for nonnormalised model average PE values, in which both the thresholds and the PE trended downwards for more deeply inserted electrode contacts. Subsequent multivariate regression analysis of the normalised PE values reinforced the observation that normalised PE mainly correlates with electrode contact insertion angle. The effect of normalised average thresholds (the mean of the I 1mm levels of TA and TC pulses) on normalised PE values is less clear; it likely interacts with the insertion angle and only showed a negative significant correlation when all peripheral processes were removed.
Though Carlyon et al. reported a significant decrease in anodic and cathodic thresholds from base to apex, they did not find a significant decrease in PE values for more deeply inserted contacts. It should be noted, however, that they only looked at electrode numbers and not cochlear insertion angles; there is considerable variability in the cochlear insertion angles of a given contact number between users of the same device due to anatomical differences and surgical techniques ( van der Marel et al., 2014 ;, so it is unclear if a correlation would have emerged if insertion angle had been taken into account. Other studies have not explicitly reported on the relationship between PE and electrode contact insertion angle, though gauging from their plotted data, there does not seem to be a universal downwards trend ( Mesnildrey, 2017 ;Goehring et al., 2019 ;Jahn and Arenberg, 2019a ;Mesnildrey et al., 2020 ). However, in real-life subjects neither neural degeneration nor electrode array position would be as uniform or idealised as represented in the model; this means that the downwards trend in PE values from base to apex could be obscured in cases where the type and degree of neural degeneration varies along the cochlea, or where the array's distance from the modiolus is inconsistent, as it would make the change in PE value along the array even more erratic.
Regarding lateral to medial positioning of the electrode array, several studies also included estimations of the electrode-tomodiolus distance (EMD) based on CT scans ( Mesnildrey, 2017 ;Jahn and Arenberg, 2019a ;Mesnildrey et al., 2020 ) and found significant correlations between perceptual thresholds and (normalised) EMD, which is consistent with the findings of this study. Mesnildrey (2017) found a significant correlation between normalised EMD and PE, but this finding was not replicated in the other studies. However, Mesnildrey et al. (2020) did conclude that EMD and PE were both contributing factors to variance in measured thresholds. This is also consistent with the results of this study, though to reiterate, the model implies that electrode contact insertion angle is also an important factor. Due to the scala tympani becoming narrower towards the apex, the EMD will generally decrease with increasing insertion angle, so ideally, they should both be accounted for. A more thorough analysis of clinically obtained PE measurements, EMD's and their relation to electrode insertion angles might be able to confirm the model's findings, but currently these predictions cannot be directly compared to literature.
Finally, although the model results imply that the state of the neurons cannot reliably be deduced from measuring the PE on a single contact, especially when important spatial factors such as electrode insertion depth and distance from the modiolus are unknown, it must be noted that there were statistical differences between the modelled neural conditions when looking at the combined data from all electrode contacts across all model geometries. In Figs. 6 and 8 there are notable differences between averages and ranges/standard deviations of (normalised) PE values between different neural conditions, so it may be possible to observe similar statistical differences between groups of CI subjects with comparable states of neural health.

Conclusions
The main conclusion of this study is that polarity sensitivity was not a reliable measure of neural health. While neural degeneration did affect simulated polarity sensitivity, its effect was not consistent; polarity sensitivity was not simply a product of the state of the neurons, but also depended on spatial factors (particularly array type and electrode contact insertion angle). Moreover, these spatial factors are at best poorly quantifiable in imaging methods available today, so in light of this study's findings, it seems difficult or impossible to determine with certainty what the underlying reasons for clinically measured polarity sensitivity data are. However, the results do suggest that there may be statistical differences between groups of subjects/electrodes with similar types and degrees of neural damage, but estimating neural health from the polarity sensitivity of an individual electrode contact seems unfeasible.
Additionally, the other notable findings of this study are: (I) In agreement with previous studies, anodic pulses mostly excited auditory neurons in their central axons, while cathodic stimuli generally excited neurons in their peripheral processes or near their cell bodies.
(II) As a consequence, cathodic thresholds were generally more affected by neural degeneration than anodic thresholds were.
(III) Measuring polarity sensitivity using charge-balanced multiphasic pulses as approximations of monophasic stimuli produced different results than those obtained with true monophasic pulses. The degree to which this occurred depended on the specific pulse shape, but in general the lower the amplitudes of the charge balancing phases were, the closer the polarity sensitivity was to that of true monophasic pulses.
(IV) Differences in polarity sensitivity between humans and animal can be explained by differences in the amount of myelin around the cell bodies of the auditory neurons. Increasing the number of myelin wraps around the cell bodies of human neurons increased sensitivity to cathodic pulses while decreasing anodic sensitivity.