Efficient Time-domain Simulation of Nonlinear, State-space, Transmission-line Models of the Cochlea (l)

Nonlinear models of the cochlea are best implemented in the time domain, but their computational demands usually limit the duration of the simulations that can reasonably be performed. This letter presents a modified state space method and its application to an example nonlinear one-dimensional transmission-line cochlear model. The sparsity pattern of the individual matrices for this alternative formulation allows the use of significantly faster numerical algorithms. Combined with a more efficient implementation of the saturating nonlinearity, the computational speed of this modified state space method is more than 40 times faster than that of the original formulation. V C 2015 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 Unported License.


I. INTRODUCTION
Active nonlinearity is an essential feature of cochlear mechanics that is generally best modelled in the time domain.A range of time domain numerical methods have been proposed and applied to a number of models over the years (Diependaal et al., 1987;Elliott et al., 2007).One widely adopted two-stage strategy (Diependaal et al., 1987) first solves the model boundary value problem using the finite difference approximation, in order to compute the pressure difference, and then uses numerical integration techniques to solve the remaining initial value problem.More recently, a state space matrix formulation has been proposed in Elliott et al. (2007), which has also been applied in several studies such as Liu and Neely (2010), Sisto et al. (2010), Rapson et al. (2012), and Ayat et al. (2014).It expresses the model states as a single set of coupled first-order ordinary differential equations (ODEs), the time domain solution of which can then be obtained using well-established ODE solvers.An extra benefit of this approach is that it allows rigorous inspection of the stability of linearized models, which, for instance, facilitates the study of spontaneous otoacoustic emissions (Ku et al., 2009).The relation between the state space formulation and Diependaal's method has previously been discussed in Rapson et al. (2012).
Despite the advantages of time domain solutions, they are usually very computationally challenging and memory demanding, because a sufficient spatial and time domain resolution is required to ensure simulation accuracy.Taking the standard state space method as an example, it is typically thousands of times slower than real time when implemented in the time domain on a desktop computer.A hybrid directiterative solver was developed in Bertaccini and Sisto (2011) to accelerate the simulation of a nonlinear feed-forward model within the standard state space framework.This letter proposes a modified state space (MSS) method for the time domain simulation of nonlinear one-dimensional (1D) transmission-line cochlear models.Numerical efficiency of this MSS is compared with both the standard state space (SS) and Diependaal's method, and it is found to be more than 40 times more efficient than the SS and marginally faster than Diependaal's method.

II. THE EXAMPLE COCHLEAR MODEL
The example model adopted here for comparison of algorithm complexity is that reported in Ku et al. (2009), which is a nonlinear adaptation to the linear active model originally proposed by Neely and Kim (1986) for a cat cochlea, but with parameters re-tuned to match human cochlear physiology.The micromechanics of this model is represented by an array of 500 discrete elements, with the first one modelling the middle ear, the last one the helicotrema and the remaining 498 elements modelling the cochlear partition, as shown in Fig. 4 of Elliott et al. (2007).The cochlear partition elements are modelled as two degree-of-freedom (DOF) lumped-parameter oscillators [Fig. 1 of Elliott et al. (2007)], coupled solely by the cochlear fluid.The saturating active mechanism is simulated by placing a compressive nonlinearity, a first-order Boltzmann function, before the active impedance in the micromechanical feedback loop, so that the input to the Boltzmann function is the difference between the displacements of the two masses of each oscillator.The dimensionless factor, c, which regulates the strength of the active feedback force, is then determined as the absolute value of the ratio of the output and input of the Boltzmann function, as shown in Eq. ( 8) of Ku et al. (2009).A two-port network model of the ear canal and middle ear, based on the model of Kringlebotn (1988) and programmed by Ku (2008), is also included.A thorough description of the model and its parameter values can be found in Ku (2008).

III. MODIFIED STATE SPACE FORMULATION
The original state space formulation proposed by Elliott et al. (2007) used the finite difference approximation for the 1D wave propagation inside the cochlea, Eq. ( 1), and state space formulation for the element micromechanics, Eqs. ( 2) and (3), so that the final fluid-coupled state space equation can be derived as Eq. ( 4), FpðtÞ À € wðtÞ ¼ qðtÞ; (1) (2) where F is the second-order finite difference matrix, representing the coupling between the cochlear partition acceleration vector € wðtÞ and fluid pressure difference vector p(t); q(t) is the stapes acceleration vector, the only nonzero entry of which is the first one; x(t) is the complete state vector; is the input matrix.Complete derivations and definitions of these matrices can be found in Elliott et al. (2007).An alternative method for deriving the final state space equation, with much better computational efficiency, is to first combine Eqs. ( 2) and (3) to eliminate variable _ xðtÞ, and then combine the result with Eq. ( 1) to eliminate dependence on € wðtÞ, Equating the right hand side of Eqs. ( 5) and ( 6) leads to the pressure difference vector as Substituting Eq. ( 7) back into Eq.( 2), the modified state space equation is obtained as where the new system matrix becomes . This new set of state space matrices A 0 , B 0 can be shown to be analytically identical to the original ones A, B, using the Woodbury matrix identity for the inverse of a sum of matrices (Henderson and Searle, 1981).Applying this identity to the original system matrix A, we have Similarly, for the original input matrix B, we have

IV. NUMERICAL ADAVANTAGE OF THE MODIFIED STATE SPACE METHOD
As shown in Eqs. ( 10), ( 21)-( 23) of Elliott et al. (2007), F is a tri-diagonal matrix (i.e., having nontrivial elements only on the diagonal, first sub-diagonal and first super-diagonal), A E , B E , and C E are block diagonal matrices (in which the diagonal elements are square matrices of any size and off-diagonal elements are zero).It is therefore expected that considerable memory and calculation savings can be achieved by employing sparse matrix storage and computation methods, especially when several hundred micro-elements are used.The sparsity patterns of the state space matrices are shown in Fig. 1, where the upper two rows show the individual matrices and the bottom two rows both those inverted and the final matrices in the state space formulation.Each square represents the structure of a matrix and each dot denotes a nonzero element at the corresponding location.Only the first 20 Â 20 elements of each matrix are depicted for better visualization of their internal patterns.The size of each matrix for a 500-element model, the number of nonzero entries (nnz) and sparse matrix density, d, defined as the ratio of nnz to the total number of entries, are also given in the figure caption.It can be seen that more than 99% of the matrices A E , B E , C E , F, and F À C E B E are empty, whereas this number reduces to about 75% for the final set of state space matrices A and B. In MATLAB, sparse algorithms (Gilbert et al., 1992) can be employed when these matrices are defined to be sparse and the resultant sparse state space method is about twice as efficient as the original formulation.The modified state space formulation, however, can be arranged to be even faster with sparse algorithms, because the intermediate matrix F À C E B E is tri-diagonal and the calculation of the pressure difference vector p(t) in Eq. ( 7) is identical to solving a sparse tri-diagonal linear system of the form Ax ¼ b.There exist very rapid algorithms for computing the solution of such linear system using either iterative or direct methods.Direct solvers based upon sparse matrix factorization are preferable in this case due to the dimension of the model and the nature of matrix F À C E B E .In MATLAB, the matrix left divide operator, mldivde or backslash (\), encapsulates a host of algorithms for solving sparse linear systems and invokes the most suitable one according to the sparsity pattern of involved matrices (Davis, 2006).One solver is specifically designed for tri-diagonal systems, and direct use of this algorithm leads to a speed improvement of at least 20 times compared to the original formulation.But since the tri-diagonal matrix, F À C E B E , is fixed during the entire simulation, the speed of solving this linear system can be further increased by performing a sparse LU decomposition beforehand.Essentially, this converts the tri-diagonal system into a combination of a lower and an upper triangular system which can be easily solved using forward and backward substitution.According to (Davis, 2004), this decomposition for the matrix F À C E B E takes the form of where L is a lower triangular matrix, U is an upper triangular matrix, P and Q are row and column permutation matrices used to reduce the fill-in problem during the actual factorization process, and R is a diagonal row scaling matrix which can lead to a sparser and more stable factorization.Thus, in addition to using sparse matrix multiplication, the modified state space equation is implemented as the following: qðtÞ, all of which can be computed before calling the ODE solver.Such computational advantages are not possible with the original state space method, as its intermediate matrix, I À B E F À1 C E contains significantly more nonzero entries (more than 167 times) and is of larger size (about 16 times) than those of F À C E B E .The reduction of sparsity mainly results from F À1 , which is a completely full matrix.

V. INCLUSION OF NONLINEARITY IN THE STATE SPACE EQUATIONS
Only the active feedback force needs to be compressed with stimulus level, so the matrix A E can be decomposed into a time-invariant passive part, A E pas and a time-varying active part, A E act (c), which consists of all the components that are functions of c.One way of implementing the level-dependent nonlinearity in the time domain is to update A E act (c) every time c changes.This involves several modifications to a 1996 Â 1996 matrix for every time step when using the ODE solver and can be extremely time-consuming.By observing that each individual block matrix inside A E act (c) only contains nonzero entries on its first row, which also have a common factor c, an equivalent implementation is to scale all the components of the state vector for each cochlear partition element by a factor c(n, t), which is given by where n is the element index; x d (n, t) is the relative displacement between the BM and TM; and f(x) is the Boltzmann function.The active part of the state space equation, A E act (c)x(t), can now be written as A E full x scaled ðtÞ, where , is a constant matrix, x scaled ðtÞ ¼ cðtÞ xðtÞ, cðtÞ is a column vector including all of the scaling factors, c(n, t), arranged in the same order as are the elements of the complete state vector x(t) and the symbol denotes element-by-element multiplication of the two column vectors c(t) and x(t).Therefore, the final fluid-coupled state space equation for the nonlinear cochlear model is realized as the following: This method is considerably more effective than the first one, as none of the state space matrices is changed during the numerical integration.The original SS method for the nonlinear model is implemented in a similar way for comparison of numerical efficiency.

VI. RESULTS
This section presents a comparison of the computational efficiency of four time domain numerical algorithms: the SS, sparse state space (SSS), modified state space (MSS), and Diependaal's two-stage method after extending it to allow two DOF micromechanics following the Appendix of Diependaal et al. (1987).All of these were programmed and simulated in MATLAB R2015a using a desktop computer with a 3.40 GHz, quad-core Intel Core i5-3570 processor and 4 GB DDR3 RAM.The final component of all four algorithms is numerical integration.We implemented an explicit adaptive solver, which is a modified version of the MATLAB ode45 function, removing unnecessary functionalities such as event function, mass matrix and non-negative solutions.The relative error tolerance was set as 10 -5 in all experiments, whereas the absolute error tolerance was determined individually for each stimulus to be one order of magnitude lower than the corresponding scale of the model displacement response.This reflects the flexibility of the adaptive solver, which can solve the model with the required accuracy and without unnecessary effort.The effect of the stimulus spectrum and level on program runtime was investigated by using white noise and sinusoids having frequencies of 2, 4, 8, 12, and 16 kHz, each of which was sampled at 100 kHz, and had a duration of 30 ms with a 10 ms half-Hanning window onset ramp and levels varying from 30 to 100 dB sound pressure level (SPL) in step of 10 dB.
For each stimulus, the simulation time gradually rises with increasing input level and frequency for sinusoidal stimuli.But the speedup factor, defined as the ratio of the SS simulation time to that of any other method, is roughly constant across most stimulus types and levels, as shown in Fig. 2 for the MSS.Similarly, the speedup factor is found to be relatively independent of the number of microelements in the model (not shown).Table I shows the average runtime across different signal levels and frequencies of sinusoidal excitations for each type of stimulus and algorithm.The overall average runtime of all stimuli are taken as the final metric for comparison of the computational efficiency of each method and the overall speedup factor is shown in the last row of Table I.It can be seen that the sparse algorithms alone are able reduce the runtime by almost a half, but the MSS yields a speed improvement of more than a factor of 40 and is also about 40% quicker than the two-stage method.

VII. CONCLUSIONS
The main contribution of this paper is the description of the modified state space method, which has been applied here to the nonlinear and active cochlear model developed in Ku et al. (2009).Although analytically identical to the original state space equation, the sparsity pattern of the constituting matrix of this alternative formulation offers the opportunity for considerably more efficient numerical algorithms, producing a speedup factor of more than 40.Its computational efficiency is on a similar scale to that of the Diependaal's method.The approach presented here can be readily applied to various other 1D cochlear models, such as the nonlinear and active ones based on an array of one DOF lumped-parameter oscillators (Sisto et al., 2010).

FIG. 2
FIG. 2. (Color online) Dependence of speedup factor for the modified state space method on stimulus SPL for white noise and sinusoidal stimuli with frequencies of 2, 4, 8, 12, and 16 kHz.

TABLE I .
Comparison of average run-time in seconds, used by four algorithms to solve the nonlinear 1D cochlear model with 500 discrete elements in response to white noise and sinusoidal stimuli with frequencies of 2, 4, 8, 12, and 16 kHz, each of which was sampled at 100 kHz, had a duration of 30 ms with a 10 ms half-Hanning window onset ramp and levels varying from 30 to 100 dB SPL in step of 10 dB.SS: original dense state space; SSS: sparse state space; MSS: modified state space; two-stage: Diependaal's method.