Algorithm for multi-curve-fitting with shared parameters and a possible application in evoked compound action potential measurements
© Spitzer et al; licensee BioMed Central Ltd. 2006
Received: 29 November 2005
Accepted: 22 February 2006
Published: 22 February 2006
Experimental results are commonly fitted by determining parameter values of suitable mathematical expressions. In case a relation exists between different data sets, the accuracy of the parameters obtained can be increased by incorporating this relationship in the fitting process instead of fitting the recordings separately.
An algorithm to fit multiple measured curves simultaneously was developed. The method accounts for parameters that are shared by some curves. It can be applied to either linear or nonlinear equations. Simulated noisy "measurement results" were created to compare the introduced method to the "straight forward" way of fitting the curves separately.
The analysis of the simulated measurements confirm, that the introduced method yields more accurate parameters compared to the ones gained by fitting the measurements separately. Therefore it needs more computer time. As an example, the new fitting algorithm is applied to the measurements of the evoked compound action potentials (ECAP) of the auditory nerve: This leads to promising ideas to reduce artefacts generated by the measuring process.
The introduced fitting algorithm uses the relationship between multiple measurement results to increase the accuracy of the parameters. Its application in the field of ECAP measurements is promising and should be further investigated.
It's very common to analyse a system by making measurements and trying to fit a mathematical equation (i.e. the model of the system) to the results. That way the system is described by the fitted parameters. If the system is complex and the equation has many (M) parameters a 1 to a M (described by the vector a), receiving usable values from only one measurement/fit is difficult. In this case, one possibility is to make more than one (N) measurement and alter some of the test-conditions which should appear in one or more parameters of the fittings. If the results of these measurements were fitted separately, the situation would not improve much. Here an algorithm is introduced that fits those N measurements simultaneously to N equations which may be (but need not) different and may share some of the parameters a m .
The fit method we are introducing here is able to fit these curves simultaneously and takes into account that the parameters k and c are shared. Therefore it returns only one value for k and one value for c (figure 1(d)). The following equations were used for fitting:
y 3(x) = g3 e-kx + c (1)
y 2(x) = g 2 e-kx + c (2)
y 1(x) = g 1 e-kx + c (3)
In this example, the vector of the fit-parameters a was assigned as follows:
a 1 = k a 2 = c a 3 = g 1 a 4 = g 2 a 5 = g 3
Figure 1(c) shows the result when using this algorithm with shared k only.
Principle of the algorithm
The basic idea is the following: When fitting one curve to one equation, the goodness of the fit parameter χ2 (weighted sum of the quadratic deviations of the values – definition follows) is minimised. To fit N curves to N equations simultaneously, the sum of the individual χ2 values has to be minimised. Definitions:
N ... Number of equations/curves. Each curve is represented by one equation.
A(n) ... Number of the nodes belonging to the n th measured curve.
(x ni , y ni , ) ... node number i of the n th curve (i = 1 ... A(n)). σ2 is the uncertainty and can be set to 1 if it is identical for all nodes of all measurements.
M ... Total number of parameters that should be determined with the fit.
a= (a 1, ..., a M ) ... M element vector of the fit-parameters a 1 to a M .
y n (x, a) ... n th equation. It can use all fit parameters a, but there are three possibilities the n th equation could use every single parameter a m :
The equation is independent of a m . Then the (in the following needed) partial derivation of the equation according to this parameter a m is zero.
a m is a "normal" fit-parameter, that should be optimised.
The equation depends on a m , but should not partake in the optimisation. This is done by setting the partial derivation for this parameter to zero.
The total sum of the goodness of the fit parameter χ2 is defined as sum of the goodness of the fit parameters of each curve. The bigger χ2, the worse the fit:
Method for linear functions
If all equations y n (x, a) are linear in the parameters a, we can use the least square fit method  to minimise the total χ2. In this case, the equations y n (x, a) have to be of the following form:
f nm (x) are arbitrary functions. Note that a polynomial function is a special case of this where f nm (x) has e.g. the form f nm (x) = x m . To get the values for a, the value of χ2 is minimised by setting the derivations according to each parameter a k to zero:
In this way one gets M linear equations (k = 1 ... M) to determine the M entries of a. Equation 10 may look complicated, but all values except a m are explicitly known. The equation can be solved by the method of determinants.
The advantages of the least square method over all non-linear methods are very fast processing, final results in one step and no need to specify starting parameters.
Method for nonlinear functions
Calculate χ2 with suitable starting parameters a with equations (4) and (5).
Set λ to 10-3.λ is a constant factor that controls whether the Marquardt-Fit-Method should behave more as a gradient search fit method (λ ≪ 1) or an expansion fit method (λ ≫ 1).
Calculate the vector δa, which has the same size as a(namely M elements) and describes the suggested correction to a as follows:
The M element vector β represents the first partial derivations of χ2 to the fit parameters described by a:
The M × M Matrix describes the second partial derivations to the fit parameters a k and a m (k and m vary from 1 to M respectively):
Calculate the matrix α by multiplying each diagonal element from by (1 + λ).
Determine the inverse matrix ε of α.
Calculate δa as follows:
Derive new trial-fit-parameters a' from the "old" ones a:
Determine χ'2 with the trial-fit-parameters a' using equations (4) and (5 and (5)).
If χ'2 ≥ χ2 (i.e. there is no improvement with the trial-fit-parameters): Substitute λ with 10·λ, keep a unchanged, and continue with step 3.
If χ'2 <χ2 (i.e. the new parameters are better): Substitute λ with λ/10, and a with a' and continue with step 3.
This iterative calculation may be stopped when one of the three following conditions is met:
χ2 falls below a certain predefined threshold.
The difference between χ2 and χ'2 decreases to less than a specified value.
The count of the calculation loops exceeds a maximum number.
Note: Using this method, one has to treat "constant" parameters (parameters that do not participate in the fitting process) in a special way (and not within a having derivations of zero), because otherwise their value may be biassed to minimise χ2. This could be done by an extra parameter for y n (x, a) (e.g. called b) or by incorporating the constant parameters within the functions y n (x, a).
Comparsion to single fits
Use with cochlear implant ECAP measurements
Cochlear implants are medical devices that enable deaf people to perceive hearing impressions by electrically stimulating the auditory nerve . The stimulations occur through a multichannel electrode that is placed directly into the cochlea and has contacts on different positions. This way, up to 100% speech recognition can be achieved .
Modern cochlear implants can record the evoked compound action potential (ECAP) of the auditory nerve. This is done by measuring the voltage at one of the contacts of the multichannel electrode after a stimulation pulse that evokes the action potentials. The recordings obtained show the sum of the ECAP and exponential decays resulting from residual charges (arising from capacitive components on the implant itself or from the electrode-electrolyte transition). There are several well established recording methods to eliminate or reduce this so called artefact to obtain the pure ECAP: Alternating stimulation , masker probe methods [6, 7], triphasic stimulation  and scaled template methods . However, each method has its limitations: Some of them rely on the linearity of the system (alternating stimulation), others need up to between nine to sixteen times more measurements for one result (masker probe methods).
Another issue is the influence of the used measurement system itself, which may be seen on the measured curves as offsets, drifts, additional exponential decays, ... (called system effects here).
The introduced method can be used as a tool to help separate the ECAP signal from these undesired components. One way to do so is to fit just the artefact (without the ECAP) along with the system effects and subtract this fit from the measurements. This can be done with sub-threshold measurements, recovery measurements, or in regions, without an ECAP. Another possibility is to combine this method with an artefact cancellation method mentioned above, for example to eliminate system effects.
The algorithm needs sequences with different conditions to gain advantage over single-curve-fits. In practical use this is no disadvantage, because usually two special sequences are measured: The amplitude growth sequence and the recovery sequence. The amplitude growth sequence raises the stimulation pulse amplitude from zero to the maximum comfortable loudness (MCL) level. The recovery sequence sends two MCL level pulses before the measurement – if the second pulse lies within the auditory nerve recovery time of the first pulse (about < 1 ms), the ECAP vanishes.
Example ECAP measurement
The following equation was used as "model" of the artefact plus offset for each curve. We assumed that the artefact can be described by two exponentially decreasing terms and one offset parameter.
y(t) = fe-t/τ + ge-t/T + c (15)
The final offset c and the amplitudes f and g are different for each curve, but the time constants τ and T are assumed to be global, so that the following a was used. It contains 32 parameters.
a= (τ, T, f 1, g 1, c 1, f 2, g 2, c 2, ..., f 10, g 10, c 10) (16)
As mentioned above, several ideas exist to take advantage of the algorithm. Three of them will be described in the following subsections:
Fit of the regions without an ECAP signal
The first possibility is to exclude the regions of the curves where an ECAP signal is expected from the fit. In this example, we included the whole range of the five measurements from 0 to 355 cu stimulation amplitude that are sub-threshold (figure 3), and the right part of the other five measurements, where the ECAP signal (duration of about 1 ms ) has vanished.
These patient specific time constants describe mainly the electrical properties of the electrode-electrolyte interface  and the geometry, and could be taken as diagnostic parameters along with the ECAP amplitude or the residual voltage measured in telemetry measurements.
A remark on the starting parameters: The method we used to get the starting parameters was to fit the curves with a simplified model having only one exponential term (with shared time constant) and taking into account only the right halves of the curves. In this way we obtained values for the long time constant τ, f i and c i . We then subtracted the fits of this simplified model from the curves and fitted again, using only one exponential term with a shared time constant T and taking into account the whole curve. After that we got starting values for g i and T. These additional fits needed only a few seconds calculation time.
Using σ2to fit over all regions
Instead of first excluding the parts of the curves with an ECAP signal and then fitting over all regions where some of the previously determined parameters are kept constant, we could initially fit over all regions and use the uncertainty parameter σ2 (see section "Principle of the algorithm") to characterise regions where the ECAP signal is expected: These regions should receive a high σ2 value compared with the σ2 value from the other regions, so that the fitting algorithm doesn't punish the ECAP caused deviations with a high contribution to the χ2 value.
The more of the following conditions are satisfied, the better the results from this approach:
The ECAP amplitude is small compared to the amplitude of the artefact.
There is a region before and after the ECAP-part of the curve that is described by the artefact model used.
The ECAP signal is dc-free, i.e. the integral over the ECAP signal (without any artefact) is very small.
The artefact time constants are larger than the ECAP time constants.
Using an ECAP model to fit over all regions
A further approach is to model not only the artefact, but the ECAP signal too.  would be an example for a model that could be used. This model takes into account the double-peak shape of some ECAPs as well. Doing so is part of current investigation and may be published in the future.
The introduced fit-algorithm uses the additional information of the relation between measured curves to retrieve more accurate parameters compared to the parameters extracted from single fits. In the case of ECAP measurements the algorithm could be used as (additional) artefact cancellation method with the following benefits:
There are no additional measurements necessary to apply this method.
The method can be used in combination with other artefact cancellation methods.
It can take into account implant specific system effects.
Physiological or implant specific parameters like time constants or the artefact amplitude are gained as values that can be used for diagnostic purposes.
No additional noise is added.
It is very flexible because of the possibilities to judge different regions of a curve and different curves in special ways with the σ2 parameter (e.g. curves with zero pulse amplitude fit some of the parameters very accurately), leave out data which should not be included in the fit (e.g. holes where the ECAP is expected), use different equations for different curves (e.g. combining results of amplitude growth and recovery sequences), and exclude parameters of some curves from the fit process (e.g. time constants in curves with ECAP).
Constant parameters (time constant, offset, ...) can be reused (at least as starting parameters) to save calculation time.
The drawbacks are that the calculation can be time-consuming (especially with many fit parameters) and that there have to be suitable starting parameters. The first issue can be improved by reducing the resolution of the curves to receive a rough result, and optimising the stop condition for the calculation loop to fit the problem. The second issue can be dealt with by reusing parameters from earlier runs, or by determining them by doing rough pre-fits with less complicated equations.
Further investigation is necessary to develop a final method based on the introduced ideas. ECAP signals obtained with this method should be compared in size and shape with results of traditional artefact cancelation methods.
The authors would like to thank the MedEl company for their help and for providing the new PULSARCI100 implant and Hansjörg Schösser for being the contact person for implant related questions.
- Bevington PR, Robinson KB: Data Reduction and Error Analysis. McGraw-Hill; 2003.Google Scholar
- Marquardt DW: An algorithm for Least-Square Estimation of Nonlinear Parameters. Journal of the Society for Industrial and Applied Mathematics 1963,2(2):431–441. 10.1137/0111030MathSciNetView ArticleGoogle Scholar
- Hochmair ES, Hochmair-Desoyer IJ, Zierhofer CM: Electronic Circuits for Cochlear Implants. AEÜ – International Journal of Electronics and Communications 1990., 44: Google Scholar
- Helms J, Muller J, Schon F, Moser L, Arnold W, et al.: Evaluation of performance with the COMBI40 cochlear implant in adults: a multicentric clinical study. ORL J Otorhinolaryngol Relat Spec 1997, 59: 23–35.View ArticleGoogle Scholar
- Eisen MD, Franck KH: Electrically Evoked Compound Action Potential Amplitude Growth Functions and HiResolution Programming Levels in Pediatric CII Implant Subjects. Ear & Hearing 2004,25(6):528–538. 10.1097/00003446-200412000-00002View ArticleGoogle Scholar
- Brown C, Abbas P, Gantz B: Electrically evoked whole-nerve action potentials: data from human cochlear implant users. The Journal of the Acoustical Society of America 1990,88(3):1385–1391. 10.1121/1.399716View ArticleGoogle Scholar
- Miller CA, Abbas PJ, Brown CJ: An improved method of reducing stimulus artifact in the electrically evoked whole-nerve potential. Ear & Hearing 2000,21(4):280–290.View ArticleGoogle Scholar
- Zimmerling M: Messung des elektrisch evozierten Summenaktionspotentials des Hörnervs bei Patienten mit einem Cochlea-Implantat. PhD thesis. Universität Innsbruck, Institut für Angewandte Physik; 1999.Google Scholar
- Miller CA, Abbas PJ, Rubinstein JT, Robinson B, Matsuoka A, Woodworth G: Electrically evoked compound action potentials of guinea pig and cat: responses to monopolar, monophasic stimulation. Hearing Research 1998,119(1–2):142–154. 10.1016/S0378-5955(98)00046-XView ArticleGoogle Scholar
- Klinke R, Silbernagl S: Lehrbuch der Physiologie. Stuttgart: Thieme; 2003.Google Scholar
- Geddes L: Historical evolution of circuit models for the electrode-electrolyte interface. Annals of Biomedical Engineering 1997, 25: 1–14.View ArticleGoogle Scholar
- Lai WK, Dillier N: A Simple Two-Component Model of the Electrically Evoked Compound Action Potential in the Human Cochlea. Audiology Neuro-Otology 2000, 5: 333–345. 10.1159/000013899View ArticleGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.