A phenomenological computational model of the evoked action potential fitted to human cochlear implant responses

There is a growing interest in biomedical engineering in developing procedures that provide accurate simulations of the neural response to electrical stimulus produced by implants. Moreover, recent research focuses on models that take into account individual patient characteristics. We present a phenomenological computational model that is customized with the patient’s data provided by the electrically evoked compound action potential (ECAP) for simulating the neural response to electrical stimulus produced by the electrodes of cochlear implants (CIs). The model links the input currents of the electrodes to the simulated ECAP. Potentials and currents are calculated by solving the quasi-static approximation of the Maxwell equations with the finite element method (FEM). In ECAPs recording, an active electrode generates a current that elicits action potentials in the surrounding auditory nerve fibers (ANFs). The sum of these action potentials is registered by other nearby electrode. Our computational model emulates this phenomenon introducing a set of line current sources replacing the ANFs by a set of virtual neurons (VNs). To fit the ECAP amplitudes we assign a suitable weight to each VN related with the probability of an ANF to be excited. This probability is expressed by a cumulative beta distribution parameterized by two shape parameters that are calculated by means of a differential evolution algorithm (DE). Being the weights function of the current density, any change in the design of the CI affecting the current density produces changes in the weights and, therefore, in the simulated ECAP, which confers to our model a predictive capacity. The results of the validation with ECAP data from two patients are presented, achieving a satisfactory fit of the experimental data with those provided by the proposed computational model.


Introduction
The best-known and most successful implantable neurostimulator or neuroprosthesis is the CI, designed to provide a sense of sound to an adult or child with severe to profound hearing loss [1]. CI is composed of a microphone, a speech processor, a transmitter and receiver/stimulator, and an electrode array, which is inserted in the cochlea. A recent survey of the state-ofthe-art of CIs could be found in [2].
In the past two decades, assessment of CI functionality is clinically recognized to be widely used through the electrically compound action potential (ECAP) (see for example [3]), which represents the synchronous firing of a population of electrically stimulated ANF. A clinical applicability of ECAPs is to obtain the threshold current to evoke acoustical perception and the maximum comfortable loudness levels before producing any kind of discomfort [4]. A recent overview of recording methodologies, response characteristics and potential applications of the ECAP is published in [5], highlighting the importance of this measure in CI clinical practice.
A crucially important aspect in neuronal modeling is to know the spread of excitation of nerve activity that is excited when it is subjected to neuronal stimulation produced by a CI. It is essential to predict the activation of nerve fibers out of the target region of each electrode of the CI (crosstalk) [6][7][8] and to design electrodes that improve the focalization [9][10][11][12]. Also, it could be useful to predict dead regions of the cochlea [13,14]. Specifically, we want to know the relationship between the current applied by each electrode and the set of neurons of the auditory nerve that are excited.
The complexity of biological systems makes it impossible to build a deterministic model that provides a precise response to an external stimulus. As it is pointed in [15], "Most models of neural response to electrical stimulation, such as the Hodgkin-Huxley equations, are deterministic, despite significant physiological evidence for the existence of stochastic activity. For instance, the range of discharge probabilities measured in response to single electrical pulses cannot be explained at all by deterministic models." As far as we know, the current computational cochlear models have a deterministic approach, that is, they try to accurately reproduce all the physical and biological processes involved in cochlear stimulation [6,[16][17][18][19][20][21][22][23]. This approach would be the appropriate if we knew precisely all the variables (conductivities of the media, patient-specific cochlear geometry, neuron positions in the nerve and its synapses positions of the electrodes, healthy state of the auditory nerve, possible fibrosis on the electrode, etc.) implicated in the problem. In addition, it would be needed to know the precise excitation mechanisms of the neurons when they are stimulated with extracellular electrodes. However, at present these variables and mechanisms are not accurately known to faithfully reproduce the electrical behavior of the cochlea. In contrast, our work approaches the problem "phenomenologically", that is, adjusting the model to reproduce the clinical data.
Note that it is not currently possible to reproduce the clinical ECAP with a "deterministic" model due to, among other things, the great variability of clinical ECAP. The novelty of the model is the phenomenological approach of fitting simulated ECAPs by using weights that are obtained from a probability distribution.
Among the most relevant data to fit the electrical parameters of a computational model of a cochlea are the impedance matrix and the registers of the ECAP. The transimpedance matrix allows us to fit the electrical behavior of the model. The ECAP is used to conform the model to reproduce the patient's neural response.
There are a lot of works dealing with the modelization of CIs. For example, in [6,[18][19][20][21]24] the authors develop a FEM volume conduction model (VC) to predict the electrical fields inside a stimulated cochlea and to analyze the neural response of the ANFs. In [25] it is analyzed the effect of the position of the electrodes and the conductivities of different tissues on the potential distribution. The work of Hanekom [26] presents an extensive review of the 3D modelling of CI, and its applications.
Currently, there are different computational models of CIs that have been fitted using clinical data, such as transimpedances, cochlear dimensions, electrode locations, and loudness and pitch perception [10,25,27,28], and speech understanding [29]. Additionally, there are other models that simulate ECAPs responses [14,16,22,[30][31][32] as the one proposed in this work. In [30], ECAP is modeled using transfer functions that multiply the temporal derivatives of the membrane potential and are calculated using FEM. Moreover, [22] proposes a circuit model whose parameters are calculated by FEM.
The amplitudes of the ECAPs can vary significantly from one patient to another. The main difference between the above mentioned models and ours is the ability of our model to fit the ECAPs of a specific patient. Another fundamental difference is that this fitting is performed using weights that depend on the current densities. This dependency is what gives our model a predictive character. For example, a variation in the design of the CI could lead to a variation in the current densities that reach the neurons and, therefore, a change in the response of the model. The weights are calculated using a computational intelligence optimization method (evolutionary algorithm).
The aim of our computational model is to reproduce the clinical ECAP generated by the auditory nerve fibers (ANFs) when they are stimulated by a CI.

Method
Our model fits the clinical data of a specific patient. The prototype consists of a set of line sources, the virtual neurons (VNs), propagating the membrane currents given by Hodgkin-Huxley-type (H-H) models, (see e.g., [33][34][35]). The extracellular potential generated by the VN is calculated by FEM. Each VN contributes with a suitable weight to the simulated potential. The weights are assessed using a cumulative beta distribution characterized by two shape parameters. These parameters are chosen to fit the model to ECAP of the patient. The number of VNs must be large enough to determine, with sufficient precision, the extent of the auditory nerve stimulated around the active electrode. The fitting of the model to the clinical data is done using a differential evolution algorithm. As a result, we obtain a customized model of the neural response produced by the electric stimulation of the CI. The electrodes are stimulated in monopolar mode.
The density of real neurons in the auditory nerve is far greater than the density of VNs in our computational model, therefore, each VN represents a large number of real neurons. The weight of a VN can be understood as the number of real neurons, represented by each VN, that are effectively activated.
In ECAP recording, an active electrode E s produces a current stimulus that elicits action potentials in the surrounding ANFs. These action potentials are registered by the recording electrode E r . Our model establishes a connection between the current injected by E s and the potential registered by E r . Our model connects the stimulation current at E s with the computed ECAP, the potential generated by the VN and registered at E r . The procedure involves two types of FEM computations and a parameter adjustment by using differential evolution (DE): • Electrode mode: Finite Element computation of current densities at the VNs: A 3D FEM model is constructed to calculate the electrostatic potential and the currents produced by the active electrode E s . This potential determines the current densities at the VNs.
• Neuron mode: Finite Element computation of electric potential at E r : The FEM model is similar to the previous one, but in this case the sources are the VNs. That is, line sources propagating the membrane currents given by Hodgkin-Huxley-type models.
The main stages of the computational model are shown in Fig 1.

Ethics statement
The ethics committee of the Complejo Hospitalario Universitario Insular Materno Infantil approved the study prior to the data collection process. All patients agrees to participate in the study by signing the informed consent.

FEM model
Domain description. Our study is focused on the small region around the stimulating electrode. For that reason, it is not essential to construct a complete spiral-shaped representation of the cochlea. Nevertheless, as each tissue has a different conductivity, it is important that the model represents correctly the different structures that conform the cochlea.
The cochlear model is delimited by a 20 mm radius sphere containing a medium whose conductivity, σ ext , is chosen to fit the values of the transimpedances of the model to the ones provided by clinical data, thus guaranteeing a reasonable electrical behavior of our simplified model. The transimpedance matrix is obtained by measuring the potential profile recorded along all contacts when a single contact is stimulated in monopolar mode. The element Z ij of the transimpedance matrix is given by the ratio V i /I j , where V i is the voltage at the electrode i respect to the reference electrode (ground) and I j is the current delivered by the stimulating electrode j. The transimpedance matrix was recorded from a specific patient using the following configuration of the Custom Sound Evoke Potentials software from Cochlear (Cochlear Ltd. Sydney, Australia). A biphasic stimulus with a phase width of 25 μs, an interphase gap width of 8 μs and a frame period of 122 μs was set. The recording time was set at the end of the first phase. The stimulation current level was 0.636 mA. The ball electrode (MP1) was used as reference in the stimulation circuit and the implant case (MP2) as reference in the recording circuit. The neural responses were carried out with Custom Sound Evoke Potential Software tool (version 5.2).
The simulated impedances depend on the conductivity of the medium, so a change in σ ext causes a variation in the computed values of Z ij . Our objective is to determine the value of σ ext that adjusts the simulated impedances to the real ones.
Our computational model represents a section of the cochlea with a total length of 9 mm. In this case, we focus on the basal portion of the cochlea selecting the first six values Z patient = {Z 11  The reference electrode is placed close to the boundary of the delimiting sphere and its radius is r = 0.5 mm.
The computational model of the cochlea has been generated by histological data from temporal bone shown in Fig 2. The histological technique used is described in reference [36]. Fig 3A shows the straight section of the cochlea, the VNs and the electrode array used in this work. The mesh of this geometry is shown in Fig 3B. The sphere delimiting the domain is not represented in this figure.
The computational model contains 11 electrodes. They are embedded in a silicone carrier, that is a good electrical insulator, so the current inside the carrier is negligible and it can be removed from the domain. In this case the surface of the silicone constitutes the internal boundary of the domain and it is considered as an isolating boundary. Thus, the boundary of the domain O is formed by an exterior boundary, the sphere, @O ext , and the interior boundary, @O int , that is, The conductivities of the different tissues considered in this work are: σ endolymph = 1.67 S m −1 , σ perilymph = 1.42 S m −1 , σ auditorynerve = 0.3 S m −1 , σ organ of corti = 0.012 S m −1 , σ modiolar wall bone = 0.2 S m −1 and σ ext = 0.3 S m −1 . A review of the resistivities given by some authors is shown in appendix E of [37]. Additionally, a detailed study of impedances of the cochlear structures can be found in [38].
The number of VNs of our model is N = 43 and they are separated by 0.21 mm from each other.
The radius of the electrode array is 0.3 mm, the dimensions of the electrodes are 0.7 mm × 0.3 mm and the inter-electrode distance is 0.7 mm. This design is based on the CI512 electrode array from Cochlear (Cochlear Ltd. Sydney, Australia). We have used Comsol Multiphysics 5.6 with quadratic tetrahedral elements for FEM calculations. The number of tetrahedral elements of the complete domain is 866 309, the average quality (mean ratio) is 0.735 and the minimum element quality is 0.043. The number of degrees of freedom is 1171 692.
It is important to highlight that the size of the elements at the VNs must be small enough to accurately reproduce the potential abruptly fluctuating in time and space. The mesh around the VNs has been refined 3 times, giving rise to a mesh with edges of 0.018 mm at the VNs.
Formulation of the problem. Although we are dealing with a dynamic problem, the involved frequencies are very low. The stimulus rate of CI is around 1 kHz and the conduction velocity of the ANF is around 15 m/s [39,40]. Therefore, the problem can be solved by using a steady current approach [41]. The potential ϕ is given by where σ is the conductivity and S is the volume current source (see for example [42]). This approach, widely used, is known as the volume conduction model. Here, O is the domain formed by the cochlea and the surrounding sphere and it is constituted by different structures characterized by their conductivities σ i , that we assume to be constant. The current density in each structure is given by The source depends on whether the FEM is actuating as electrode or neuron mode. Electrode mode: The active electrode E s is a surface delivering a net current I 0 . As E s is part of the inner boundary, it is modeled as a Neumann boundary condition and the source S = 0 is null. The other electrodes are considered as floating potentials. Although in an ECAP recording the input current takes several values, it is not necessary to solve a FEM problem for each current. The linearity of the problem allows us to accomplish a unique FEM simulation for a given I 0 . Thus, if I k = c k I 0 is the input current at E s and c k is the proportionality constant between both currents, the solution of Eq (1) for I k is ϕ(I k ) = c k ϕ(I 0 ), being ϕ(I 0 ) the solutions of (1) for I 0 . The currents considered here are the amplitudes of the stimulus pulse of the CI, so the problem is strictly steady.
Neuron mode: In this case, we can consider each VN as a line current source with a spatio-temporal dependence given by λ(t, s) = λ(t − v −1 s), where s is the arc length of VN and v is the conduction velocity in the ANF. The line current density λ(t) is inferred from the Hodgking-Huxley model for an ANF from [35]. The determination of λ(t − v −1 s) will be detailed in a subsection below. The line current source is given by where C is the parametric curve (X(s), Y(s), Z(s)) describing the VN and s the arc length. Each VN is represented by a source S n (x, y, z, t) with a line current λ n (t, s).

Boundary conditions
Active electrode. The active electrode E s is implemented as the Neumann boundary condition where J is the current density, n is the unitary normal vector to E s , I 0 is the current delivered by the active electrode, A its surface, and σ the conductivity of the surrounding medium. Disconnected electrodes. The electrodes not actuating as terminals are disconnected. They must be considered as floating potentials because they are equipotential surfaces with an unknown potential. This type of boundary condition can be implemented using the linearity of the problem. As example, let us consider one active electrode, E 1 , injecting a current I 0 , and other electrode, E 2 , disconnected (I = 0). Now, let ϕ 1 be the solution of the potential problem (1) when the electrode E 1 is at 1 V and electrode E 2 is at zero potential. Reciprocally, let ϕ 2 be the solution of (1) when the electrode E 1 is at zero potential and electrode E 2 is at 1 V. Taking into account the linearity of the problem we can write the general solution of (1) as: The coefficients α 1 and α 2 can be calculated imposing the restrictions of our particular problem. Thus, integrating Eq (3) on the electrodes E 1 (active) and E 2 (floating), and imposing the conditions À s being I nm ¼ À s H E n n � r� m da. The solution is given by (4) after solving the unknown coefficients α 1 and α 2 of the above equation system.

Boundary condition at the surrounding sphere and reference electrode
Electrode mode. In this case the currents flow from the active electrode, E s , to the reference one, taken as ground (ϕ = 0). The surrounding sphere is considered to be an isolating surface (J � n = 0), so that the boundary condition is Neuron mode. Now, the neurons act as potential sources and the small currents emanating from them are diffused in the medium. In this situation we consider the reference electrode as a floating potential and the surrounding sphere is modeled as an asymptotic boundary condition that mimics an unbounded domain with null potential at infinity (see e.g., [43,44]). We have imposed a first order asymptotic boundary condition (ABC) on the sphere @O ext .
This is a homogeneous Robin boundary condition, where R = 20 mm is the radius of the sphere limiting the domain. The error introduced by this approximation is O(R −3 ) and it is accurate enough for our purpose.

Isolating surfaces
The boundary of the silicone carrier is taken as an isolating surface, and therefore, in the part of @O int not occupied by the electrodes, we have Modeling of the neural response Line source approximation for computing the extracellular potential. The line source approximation (LSA) is an easy and accurate approach to compute the extracellular potential at any point in space using the values of the membrane currents [45][46][47][48][49]. Usually, LSA is used to calculate an analytic expression of the electric potential produced by line segments in an homogeneous medium. In this work it is applied to deal with lines having a continuous spatiotemporal variation in a non homogeneous medium.
In a myelinated nerve fiber the highest current concentration takes place in the nodes of Ranvier (see, for example, chapter 6 of [50]). Thus, the modelization of a unique nerve fiber could be performed by placing punctual current sources in the positions of the nodes of Ranvier. Nevertheless, each VN of our model represents a great number of myelinated nerve fibers (about 100) whose nodes of Ranvier can be considered as randomly distributed along the VN. For this reason, an approach closer to reality is to take the VN as a line source current with a continuous linear density. Furthermore, this approach presents few difficulties when generating the FEM mesh.
The current density normal to the surface at the nodes of Ranvier for an ANF, J r (t), is the membrane current (ionic and leak) density, obtained from the type H-H (Wang-Buzsaki) model [51] implemented in [35] (see Fig 4A).
This current density allows us to calculate an equivalent linear current density λ(t) of a VN. Let us consider an ANF with the following morphology: the total length is l = 2.7 mm and it has N r = 135 nodes of Ranvier with a diameter d = 2 μm and length h = 1 μm (see [40] for a discussion about the ANF morphology). The internodes length is 200 μm. The extracellular current source is located at 1 mm from Ranvier node 20 with a negative stimulating current of 2 mA during 100 μs. The membrane current density, J r (t), is calculated at node 22. Then, the total charge passing across the membrane in the nodes of Ranvier per unit of time is J r (t) N r πdh, and the equivalent linear current density is λ(t) = J r (t)N r πdh/l = 4.65 × 10 −8 J r (t) A m −1 .
The propagation throughout the neuron is explicitly introduced by using the expression where v is the propagation velocity and s is the arc length of the trajectory of the neuron. A similar approach is used in [52], but considering discrete sources, placed at the Ranvier nodes, instead of a current line.
In order to show the influence of the propagation velocity in the shape of the extracellular potential, let us consider the following academic example consisting in a line extended from (0, 0, −1) mm to (0, 0, 1) mm. In this case, the arc length is s = 1 + z. The explicit calculation of the potential at any point in space due to the linear current density λ(t − v −1 s) = 6.27 × 10 −8 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ðz À zÞ where J r (t) is the membrane current density (see Fig 4A), r and z are cylindrical coordinates of the the field point, z is the source point along the line, and σ = 1 S m −1 is the conductivity of the medium. The resulting potentials at cylindrical coordinates r = 1 mm and z = 0 mm, for velocities v = 2 m s −1 , and v = 15 m s −1 , are shown in Fig 4B. Note that an increase in the propagation velocity implies a greater amplitude of the ECAP. In addition, the increase in the velocity causes the duration of the ECAP to decrease, being more similar to that of the stimulus. In our FEM model we have implemented N = 43 VNs (see Fig 3A). The linear current density of the nth VN, described by the curve C n , is: where w n is the weight assigned to this VN. We have taken a propagation velocity v = 15 m s −1 in the ANF [39,40]. This velocity value is also in consonance with the results obtained in [53].
In general, it is not possible to get an analytic expression for the arc length s of C n . To overcome this difficulty, we have constructed a parameterization s(x) by interpolating the numerical computation of true length for a set of x coordinates. In Fig 5 it is shown the potential generated at t = 0.4 ms when all the VN are turned off except the central one.
Determination of the weights. Our model must reproduce the clinical data as accurately as possible. The only measure of neuronal activity that can be recorded with a CI is the action potential registered at the recording electrode E r , that is, the ECAP. The objective of the ECAPs is double. Firstly, we want to know if there is a neural response of the region of the auditory nerve close to the stimulating electrode. The second aim is to detect the current intensity threshold that evokes an action potential. For this last reason, the ECAPs are carried out for several current intensities, I k . We have to determine the weight w n that must be assigned to each VN so that the amplitude of the potential provided by the model coincides with the amplitude of the action potential registered by the ECAP. The amplitudes are the difference between negative N1, and positive P2 peaks for both real and simulated action potentials.
First of all, a remark about the enumeration of the electrodes should be pointed out. As our computational model only represents a portion of the real CI, we always take the stimulating electrode E s as the central one and it will be considered as electrode 0. The rest of electrodes are enumerated with positive or negative numbers, depending on the relative position with respect to the central one. The local numbering adopted for the real CI is similar to the numbering of the model (considering E s as electrode 0). The numeration of the VNs follows an analogous criterion: the VN over the stimulating electrode is taken as 0 and the rest of VNs are enumerated with positive or negative numbers. Fig 6 shows three real ECAPs corresponding to three stimulation currents delivered by a basal electrode. In this case, the recording electrode is placed two electrodes away from the stimulating one (electrode 2 in our local numbering).
We denominate � e ECAP ðI k Þ to the amplitude of the action potential registered by the ECAP at the eth electrode when the stimulating electrode E s supplies a current I k . � e n is the amplitude of the simulated potential produced by the nth VN and measured at the electrode eth, when the nth VN is turned on and the other VNs are turned off. We assume in this situation that all the active VNs are fed with the same linear current density. As an example, in Fig 7 it is shown the potential computed at electrode 2 (two electrodes away from E s ) when all the VNs, except the central one, are turned off. Comparing Figs 6 and 7 it can be seen that the shapes of both ECAPs are similar. Furthermore, the time lapse between N1 and P2 is also similar, being around 0.4 ms for real ECAPs and 0.5 ms for the computed ECAP.
The weight w e n ðI k Þ by which � e n must be multiplied to be equal to � e ECAP ðI k Þ is given by: This equation reflects that we only need one VN to reproduce the amplitude of the ECAP measured at a given electrode e. However, a computational model with a unique VN is not able to specify what region of the auditory nerve is affected after the stimulation of the active electrode. To be able to discriminate which part of the auditory nerve is activated, we need a considerable number of VNs. If we have N VNs and a unique ECAP measure, there are N − 1 extra degrees of freedom (weights) that must be determined. This is the usual situation, in which the ECAP measure is carried out at the electrode +2 (two electrodes away from the stimulating one). If we multiply Eq (12) by an arbitrary parameter δ n and sum for all the VNs we obtain: where N 0 ¼ NÀ 1 2 , being N the number (odd) of VNs. Writing D ¼ P N 0 n¼À N 0 d n we deduce from the equation above: where � w e n ¼ d n D w e n are the new weights. Eq (14) shows that the potential registered at electrode e is the sum of the weighted contributions of the potentials produced for each VN. The arbitrary parameters δ n represent the extra degrees of freedom, introduced by the additional VNs, that will be related with a proper physical magnitude.
Suppose that for a selected number of patients we have a number of ECAPs samples N d > 1 for each current I k . We intend to train our model with these samples in order to be valid in the usual situations, where we only have a unique ECAP. Our objective is twofold. Using the expression (14), we want to interpolate the amplitude of the ECAP for a given electrode, named anchor electrode. This electrode is the one for which the computed and measured amplitudes of the ECAPs are enforced to be identical. In addition, we want the difference between the ECAP amplitudes provided by the clinical data and those calculated by the model to be as small as possible in the remaining electrodes. The achievement of the latter objective lies in the proper choice of the δ n coefficients.
Specifically, let a be the anchor electrode. Following the Eq (12) we calculate the weights w a n ðI k Þ using the data � a ECAP ðI k Þ. Introducing these weights in the second term of Eq (14) we have: where � w a n ¼ d n D w a n . This is the expected amplitude of the potential at eth electrode, but computed from the ECAP given at the anchor electrode a. Obviously, if e = a we have � e ECAP;comp ðI k Þ ¼ � e ECAP ðI k Þ, but this is not true for the rest of electrodes. The objective is to find the parameters δ n that minimize the difference between � e ECAP;comp ðI k Þ and � e ECAP ðI k Þ using the remaining N d − 1 data. In the next subsection we will explain how to calculate these parameters.
Relation of δ n with the current density. The usefulness of a computational model is linked to its ability to predict what happens when electrical or structural changes occur in the CI, such as the type of electrode stimulation (bipolar, tripolar, etc.), the geometric design of the electrodes [12], the distance between the electrode array and the auditory nerve (perimodiolar vs lateral implant), etc. The only way for the model to be sensitive to these changes is to link the weights to some physical variable that varies when the design changes. For the reasons outlined below, we consider this variable to be the current density norm.
The input current establishes a distribution of potential and current density along the VNs. In Fig 8 it is shown the current densities at the first eleven VNs (from 0 to 10) when the stimulating electrode supplies a current of −1 mA. Most of the models for neural simulation attribute the activation of the neuron to the difference of potentials along the axon [18,54]. Nevertheless, the most accurate electrodiffusion models based on the Poisson-Nernst-Planck attribute the trigger of an action potential to the ion currents normal to the surface of the fibers [47,50]. An experimental study of how the orientation of the electric field affects the central nervous system neurons is presented in [55]. In myelinated fibers the ionic flux takes place in the nodes of Ranvier. It is extremely difficult to know with precision which is the normal vector to the Ranvier nodes in the cochlear nerve. As our VNs represent a set of real ANF, we consider more realistic to link the probability of excitation of a nerve fiber to the norm of the current density, J = |J|, generated by E s . However, J variates along the VN, but we need to assign a unique weight to each VN.
It seems reasonable to think that the higher the current density reaching the VN, the greater should be δ n . This motivates us to define the parameters δ n as the probability P n that an ANF, represented by a VN, is excited. This probability P n must increase when the current density on the VN increases. This idea is inspired by the paper [15], where the authors use an integratedgaussian distribution to calculate the probability of discharge of an ANF in terms of stimulus intensity. To calculate, P n , firstly we define the probability per unit length, p n , of a section ANF to be excited with a current density less than or equal to J: where L is the length of the VN and is the cumulative distribution function of a beta distribution probability defined in the interval [J min , J max ] (see e.g., [56]) given by: where Γ is the gamma function, J min is the current density threshold below which there is no response to the input stimulus, and J max a maximum current density that will be adjusted to the particular problem. Thus, the probability of excitation of the complete ANF, and our definition of δ n , is: where C n is the curve described by the nth VN and s is the arc length along this curve. Note that, since our model is symmetric with respect to the stimulating electrode, so are the current densities. As a consequence of this symmetry, the delta parameters verify that δ n = δ −n . However, in general, this is not true for weights � w a n . The choice of the beta distribution is due to the fact that it constitutes a family of continuous probability distributions supported on a finite interval. In practice, we take J min = 0 and α, β and J max are calculated to fit the clinical data using a differential evolution algorithm. The flexibility of the beta distribution makes it very suitable for this purpose.
Once the connection between δ n and J is established, we could interpret the weights � w a n as the number of real neurons of each VN that have been activated.
Outline of the model. Let X ECAP ¼ fe 1 ; e 2 ; . . . ; e N d g be the set of electrodes where we have an ECAP measure. The flowchart of the program used to fit the model is shown in Fig 9. 1. FEM (electrode mode) is used to calculate the electric potential originated by the stimulating electrode E s when it supplies the intensities I k , k = 1, . . ., N lev . With these results, the current density at each VN, J k,n , n = 1, . . ., N is calculated.
2. The electric potential ϕ n (x, y, z) produced by each VN acting alone is evaluated with FEM (neuron mode). Afterwards, the amplitudes � e n of the potentials at each electrode where we have ECAP data e = 1, . . ., N d are computed.

3.
A particular anchor electrode a to interpolate the amplitude is chosen. Additionally, the weights w a n are obtained. 4. The differential evolution algorithm selects new values of α, β and J max .
5. The parameters δ n are evaluated using Eq (19). Then, the weights � w a n ðI k Þ and the computed ECAP � e ECAP;comp ðI k Þ Eq (15) are calculated.  This function evaluates the error between the computed and real ECAPs. A power p = 4 has been selected for our applications.

PLOS COMPUTATIONAL BIOLOGY
Steps 4 to 6, belonging to the DE algorithm, are iteratively repeated until the stopping criterion (either E < tol or iter > maxiter) is satisfied.
As a result of this algorithm we obtain the parameters α, β and Jmax that minimize the difference between � e ECAP;com ðI k Þ and � e ECAP ðI k Þ. Differential evolution algorithm for fitting the clinical data. Differential evolution (DE) is a population based stochastic search technique successfully and widely used as global optimizer inspired by the Darwinian principle of natural selection and genetic reproduction [57,58]. It is currently being considered as one of the most powerful stochastic real-parameter optimization algorithms [59]. Its application in engineering and applied sciences has contributed to optimize and solve problems in many fields (e.g.: in the bioinformatics and biomedical engineering field, [60]); particularly, in relation with the type of problem required to be solved in this work, it has been applied in parameter estimation (e.g.: [61-63]).
Our optimization problem consists of a chromosome of three variables: α, β and J max , whose minimum and maximum initial constraint limits for the variables searched are (0, 0, 0) and (2000,2000,1500), respectively, although the DE could attain values out of those limits if they improve the fitness function. After trying different DE configurations, we have observed that the problem rapidly converges and the result hardly depends on the selected configuration.
Attained weights used in the Results section were obtained with a population size of 10 individuals, with crossover probability of 1, and stopping criterion set to a maximum of 300 generations (maxiter), or alternatively, achieving a value of the fitness function (as exposed in step 6 of the above section) of tol = 10 −30 . An often used DE/rand/1/bin strategy has been set (as in the aforementioned references of parameter estimation problems: [62,63]), particularly using per-vector-dither (of F weight). As an example, in Table 1 it is shown the parameters obtained after the optimization for the first test presented in the Results section (Fig 10).

ECAP recording
The clinical data provided from the Hospital Insular de Gran Canaria correspond to two patients meeting the following selection criteria: adults (> 18 years) implanted with a Cochlear CI532 half-band perimodiolar electrode array [64], having a full insertion of the electrode array, having five or more adjacent active electrodes for the ECAPs recording and being CI holders with more than 6 months of use.
The ECAPs have been recorded by means of the software Custom Sound Evoke Potentials (Cochlear Ltd. Sydney, Australia). The Neural Response Telemetry (NRT) technique records, in a neighboring electrode, the action potential resulted from the stimulus applied on a given electrode. Forward masking [65] was used to reduce the electrical artifact produced by stimulation. From now on, instead of ECAP, the term NRT is used in next sections because it was the procedure used for recording. The CI532 is numerated from most basal to most apical. The selected stimulating electrodes are 6 (basal) and 18 (apical). The standard NRT measurement typically uses the +2 electrode to record the action potential. A broader set of data is needed in order to adjust the proposed computational model. Therefore, in this work, we have extended the number of recording electrodes until 8 (the 4 consecutive ones on each side).

Results
We have several options to choose the electrode that interpolates the NRT: the anchor electrode. We could either select any of the recording electrodes, or interpolate the value of the NRT in the stimulating electrode (electrode 0), since we do not have NRT measure in this electrode. This last choice has the particularity of producing symmetrical results with respect to the stimulating electrode, that is � w 0 n ¼ � w 0 À n . It is especially suitable when the clinical data present this type of symmetry.
The first test consists in adjusting the NRT values of the two patients using all available data. The anchor electrode chosen in this case has been a = 0. Fig 10A (upper row) shows � e ECAP ðI k Þ and the fitted values of � e ECAP;comp ðI k Þ for the three input currents I 1 = 0.41 mA, I 2 = 0.65 mA and I 3 = 1.02 mA for patient 1. The lower part of Fig  10A shows the corresponding weights � w 0 n ðI k Þ, n = {−21, . . ., 0, . . ., 21}. Fig 10B shows the same result than A, for patient 2. In general, the results show a good agreement between the real and simulated data. The greatest difference occurs at the electrode 6 of patient 2. This difference is due to the great asymmetry of real data.
The selection of the anchor influences the NRT fitting, which is evaluated using the root mean square error (RMSE). The effect of choosing the optimum anchor, defined as the one that gives rise to the minimum RMSE, is analyzed here using the real data of Fig 10. Left and right parts of Fig 11 show, for patient 2, the NRTs fitting using electrodes 6 and 18 as the stimulating ones and being electrodes +4 and +1 the optimum anchors, respectively. The optimum anchor for patient 1 is the electrode 0, which is the same as in Fig 10A  It can be appreciated that the fitting amplitude of the simulated NRT to the clinical data has improved compared to Fig 10B. Note that a central anchor imposes a symmetrical fitting for asymmetric data. It is logical that the optimal anchors (+4 and +1) are not the central ones since the clinical data are asymmetric (see Fig 11).
In Fig 12 it is shown the real and computed NRTs from the patient 1, after stimulating electrode 6 and recording at electrode 8 (+2 in local numbering). Note that the values of peaks N1 and P2 of the real and computed graphics are similar and so is the time lapse between them.
The third test consists in discarding some real data in the fitting procedure of the model. Fig 13 shows  ±1 and ±3. The anchor electrode is the optimum from the available data. It can be seen that, despite the lack of data, the agreement between the real and simulated data is satisfactory. The last test tries to reproduce the usual situation in which we only know the NRT of the electrode +2 of a certain patient. Lack of data does not allow calculating the δ n parameters for this patient. A possibility to calculate the weights � w þ2 n is to take δ n from other patient and to use the NRT of the patient under study to calculate w þ2 n . Thus, the final weights are � w þ2 n ¼ d n D w þ2 n . Fig 14 reproduces this situation. In particular, Fig 14 (left) shows the amplitudes of the NRT for the patient 1 and electrode 6 in which the δ n have been taken from the patient 2. Fig  14 (right) reproduces the same experiment, but now the δ n of the electrode 18 and patient 1 have been transferred to patient 2. Note that, even in this case, the agreement between real data and computed result is very good.
The RMSE values between clinical and simulated data are shown in Figs 10, 11, 13 and 14.

Discussion
We have developed a computational model to simulate the neural response to electrical stimulus of CIs. The fitting parameters of our model (the weights) depend on the current density distribution. This dependence gives our model the ability to predict variations in response to changes in the design or configuration of the CI. One of the most relevant characteristics of our model is its ability to reproduce the NRT values of each patient and electrode. These data can vary widely. For example, the values recorded at electrode 20 (+2 in the local numbering) when we stimulate electrode 18 with a current I 3 = 1.02 mA are 202.18 μV for patient 1 and 806.35 μV for patient 2, as can be seen in the upper right corners of Fig 10A and 10B. That is, the variation in the amplitude of the NRT is 604.17 μV. Our model reproduces the NRTs of each patient with an RMSE less than 46.46 μV (see caption of Fig 10). As far as we know, no other cochlear model published in the literature has this ability to adapt to clinical data.
The example in Fig 14 shows the ability to predict NRT values by simulating the usual situation where only the NRT at the +2 electrode is known. The NRT measurements of the other electrodes have been used only to contrast the model, but not to fit it. In this example, the δ n parameters of patient 1 are used in patient 2, and vice versa. The greatest discrepancy between the clinical and simulated NRT amplitudes is 205, 27 μV, that is, 23.7%, and it occurs in patient 2 and electrode 1 with a current I 3 = 1.02 mA (Fig 14 right). This experiment shows, to some extent, the predictive capacity of our model. Nevertheless, it would be necessary to collect data from more patients to validate such predictive capacity. We will intend to complete a larger database in the future.
The model could help in the classification of patients according to the type of pathology. For example, it could be useful to determine dead regions of the auditory nerve. Our model can also be used to determine the focusing capacity of an electrode [12]. A higher focusing capacity is linked to a lower spreading of the current and, therefore, a higher concentration of the weights around the stimulating electrode. Our straight cochlea model has certain limitations. For instance, it produces symmetric transimpedances when the central electrode is stimulated. The symmetry of the model implies that Z i0 = Z −i0 . This means that the potential registered at the electrode i is equal to the potential registered in the electrode −i, when the electrode 0 (central) is stimulated. The construction of a more realistic 3D geometric model could improve the concordance between the clinical data and its modeling. Furthermore, a patient-adapted model can be built from computerized tomography scans [66].
The conductivity σ ext has been fitted using the transimpedances of patient one, but the model could be fitted by using the specific values of each patient.
Our model considers that all the VN are activated at the same time. The biological reality is more complex and there is not a perfect synchronization in the activation of neurons. This lack of synchronization may have some effect on the recorded ECAP and could be taken into account in future work.
In this work we have considered that the propagation velocity of the action potential is uniform throughout the neuron. However, this velocity is different depending on the section of the neuron considered, due to the change in the diameters between axon and dendrite [53]. This question could be considered in future work.
The input membrane current used in this work was provided by [35]. However, another H-H-type model could be taken for providing the input current, although the weights obtained after adjusting the computed ECAP to the clinical data would be different, given that the amplitude of the membrane current could vary according to the type of H-H model selected.
In this study only ionic and leak currents were used. This is because clinical records, obtained by the NRT technique, only present the N1-P2 response due to the electrical artifact cancellation procedure (see Fig 12). If we use a multi-compartment model it is possible to include the capacitive current, but in this case appears the P1-N1-P2 complex. If we wanted to include the capacitive current, we would have to remove the part of the wave that contains P1 in order to adjust the computed ECAP to the real recording.
Our phenomenological model, by construction, will always find the optimal solution to adjust the ECAPs. As indicated in Method section, the weights represent the number of activated neurons that each virtual neuron represents. The interpretation of the weights will be more or less realistic depending on the data provided to the model (more realistic geometry, conductivities, density of the auditory nerve, etc.). The more the model conforms to reality, the more credible the interpretation of the weights will be.