An Analogue VLSI Implementation of the Meddis Inner Hair Cell Model

The Meddis inner hair cell model is a widely accepted, but computationally intensive computer model of mammalian inner hair cell function. We have produced an analogue VLSI implementation of this model that operates in real time in the current domain by using translinear and log-domain circuits. The circuit has been fabricated on a chip and tested against the Meddis model for (a) rate level functions for onset and steady-state response, (b) recovery after masking, (c) additivity, (d) two-component adaptation, (e) phase locking, (f) recovery of spontaneous activity, and (g) computational e ﬃ ciency. The advantage of this circuit, over other electronic inner hair cell models, is its nearly exact implementation of the Meddis model which can be tuned to behave similarly to the biological inner hair cell. This has important implications on our ability to simulate the auditory system in real time. Furthermore, the technique of mapping a mathematical model of ﬁrst-order di ﬀ erential equations to a circuit of log-domain ﬁlters allows us to implement real-time neuromorphic signal processors for a host of models using the same approach.


INTRODUCTION
Inner hair cells (IHCs) are mechanical to neural transducers located within the cochlea and play an important role in biological sound processing.Sound captured by the eardrum is translated into movement of the cochlear fluid, which in turn causes the basilar membrane to vibrate (see Figure 1).This vibration is converted into a neural signal by the IHCs and results in the firing of the auditory nerve cells.The transduction of the sound signal by the IHC is nonlinear and exhibits several time constants of adaptation.To mimic the IHC processing, past silicon cochleae (see [1,2]) have used nonlinear lowpass filters to produce the adaptation characteristic of the IHC.These IHC circuits responded favourably to a simple set of stimuli but failed with more complex stimuli [1].We present measurement results of an analogue VLSI implementation of the Meddis IHC model [3], which is the most descriptive computational model for the function of the IHC [4].
In the circuit presented here, we use log-domain lowpass filters [5] to map the differential equations of the Meddis IHC model to circuits on a silicon chip.This technique allows the Meddis model to run in real time.In our model, we have furthermore increased the flexibility of the Meddis model by maintaining control of all time constants while introducing independent gain controls of the signals between the filters.The circuit implements a more accurate electronic model of the IHC function than any previous circuit.It can be shown that it exhibits the correct time constants of adaptation over a large range of stimulus conditions.Direct measurements of real-time signals on a fabricated silicon chip confirm this behaviour.

THE MEDDIS IHC MODEL
The IHC function is characterised in the Meddis model by describing the dynamics of neurotransmitter at the hair cell synapse (i.e., the membrane-cleft boundary, see Figure 1).In the model, transmitter is transferred between three reservoirs in a reuptake and resynthesis process loop (see Figure 2).
The first reservoir is a transmitter factory that releases neurotransmitter at the hair cell boundary and delivers it to a second reservoir, the free transmitter pool.The amount of neurotransmitter released from the pool into the cleft is controlled by changes in the permeability of the cell membrane.This fluctuates as a function of the intracellular voltage,  which is directly related to the instantaneous mechanical stimulus amplitude.Some transmitter is lost in the cleft through diffusion and a further fraction is taken back up into the cell.Some of the remaining transmitter in the cleft stimulates the postsynaptic afferent fibre of an auditory nerve cell.The level of transmitter in the cleft dictates the probability of the nerve cell firing (spiking).The transmitter taken back up into the cell is initially reprocessed and stored in a third reservoir in preparation for delivery to the free transmitter pool.Incorporation of this third reservoir enables the model to display the type of two-component adaptation typical of real IHC.
The five equations representing the Meddis IHC model are presented as follows.
(a) In equation ( 1), the permeability of the cell membrane is represented by k(t) and A, B, and g are constants of the model.In the absence of sound, k(t) = gA/(A + B) which represents the spontaneous response of hair cells at rest: (b) The level of available transmitter in the pool q(t) depends on the rate at which transmitter is manufactured y[1 − q(t)], the rate at which it is reprocessed xw(t), and the rate at which it is lost to the cleft k(t)q(t): (c) The cleft receives neurotransmitter at a rate k(t)q(t), where some of it is lost through diffusion at a rate lc(t) and some is actively returned to the reprocessing store at a rate rc(t): (d) The reprocessing store receives transmitter at a rate rc(t) and returns it to the free transmitter pool at a rate xw(t): (e) The remaining level of transmitter in the cleft c(t)dt determines the probability of the afferent nerve firing.The constant h is used to scale the output for comparison with empirical data: prop(event) = hc(t)dt. (5)

CURRENT MODE DESCRIPTION
In order to implement the Meddis IHC model on an integrated circuit, the model equations need to be written as electrical equations.We have mapped the equations of the Meddis model to electric currents to allow the use of log-domain filters for the implementation of the differential equations (see Figure 3).Each of the signals k(t), q(t), c(t), and w(t) is now represented by currents and written as I k , I q , I c , and I w , respectively.A voltage V s is used to represent the stimulus s(t).
The first equation of the Meddis model can be approximated by a half-wave rectification (HWR) function that exhibits the spontaneous bias of the IHC and saturates at some maximum value.In our implementation, a differential pair of transistors is used to create a sigmoid function with a shape similar to (1) where n is the slope factor and U T = kT/q is the thermal voltage parameter of the MOS transistor.The three differential equations (2), (3), and (4) are rewritten as difference equations by taking the Laplace transform.The resulting equations are given as follows: where Each constant of the model is represented by a time constant with the two constants l and r combined into a single time constant t c (see (10)).The model was further simplified by replacing time constant ratios with dimensionless gains (see (11)).This step makes the model easier to manipulate and time constants can now be changed without affecting the gains and vice versa.The product I k × I q is normalised by a constant current I max to ensure that the currents remain within the same order of magnitude.

Log-domain filters
Log-domain circuits are dynamic translinear circuits [6] that use the exponential transconductance of devices such as bipolar and weak inversion MOS transistors to compress and expand current mode signals.The relationship between input and output currents is linear while a logarithmic (compressive) voltage-current relationship provides these filters with high dynamic range.We use a first-order log-domain filter, first investigated by Frey [5], implemented with MOS transistors operating in the weak inversion mode (Figure 4).This circuit can be analysed using the translinear principle.The product of drain currents of transistors facing clockwise (M1 and M3) is equal to the product of drain currents of transistors facing anticlockwise (M2 and M4) [6].For the circuit of Figure 4, this gives The summation of currents at node V c also defines Furthermore, the drain currents of M3 and M4 are related by Differentiating (14) gives or Substituting ( 16) into ( 13) then gives us or where τ = CU T /I τ .This circuit thus implements a first-order lowpass filter with a time constant determined by the value of the capacitor C, the thermal voltage U T , and a current I τ , which can be used to control the time constant after fabrication.

Translinear multiplier/divider
The translinear circuit shown in Figure 5 functions as a onequadrant multiplier/divider is used to generate the product of I q and I k normalised by I max (see (7) and ( 8)).The inputs to the circuit are I k , I q , and I max and the output is I d .The transistors M1, M2, M3, and M4 form the translinear loop and M5 is an adaptive bias transistor that actively biases M2 and M3.

The IHC circuit
Three log-domain lowpass filters and the multiplier circuit are used in our IHC implementation (Figure 6).Variable gain current mirrors connect these stages to provide off-chip control of the gains g a , g b , g c , and g d .The HWR block in Figure 6 implements (6) of the current domain model.This provides conversion from voltage (the output of the silicon cochlea described in [7]) to current and HWR.The MULT circuit implements the generation of (I k × I q )/I max .The three lowpass filters, CLEFT LPF, STORE LPF, and POOL LPF, implement (7), (8), and (9), respectively.
Mismatch analysis of these circuits [8] has revealed that they are susceptible to mismatch in the threshold voltage parameter v th .By using cascoded mirrors, their mismatch can be reduced.Section 5 shows that this mismatch may be overcome by adjusting the gain and time constant parameters and does not prevent the model from working.

RESULTS
Our aim here is to show that the Meddis IHC model [3] has been implemented on an analogue chip.To achieve this, we must find a similar parameter set to that used in Meddis' own tests.The parameters given in [9] were used in the current mode model.However, as the permeability function is slightly different (compare (1) with ( 6)), some parameters had to be adjusted.Although this was very time intensive in simulation, the chip could be tuned in a "hands-on" approach by controlling the parameters using off-chip voltages, while having a real-time response seen on the oscilloscope.The hands on approach proved to be a very fast way to find a close parameter set.
The first six tests are a subset of well-reported auditory nerve properties in response to tone-burst stimuli for which electrophysiological data exists.In [10], Meddis tests eight computational models of mammalian IHC function and finds that none replicates the IHC in all tests.The Meddis IHC model is favoured due to its good agreement with physiological data and its computational efficiency.
The soundcard output of a PC was used to create a voltage signal that represents the tone-burst inputs to the Meddis model.All tones were 1 kHz sine waves except where stated otherwise.These tones correspond to the pattern of vibration at a particular point along the cochlear partition.The output of the chip and the Meddis model represents the instantaneous probability of a spike event in a postsynaptic auditory nerve fibre, and thus are independent of any postsynaptic effects.The chip-output signal was a current below 100 nA that was amplified using a current sense amplifier to a voltage which was measured on an oscilloscope.

Sources of error
It should be noted that there are various sources of error in the results.These errors may explain discrepancies between the original model and the chip response.An estimation of the results reported by Meddis is taken from graphs presented in journal papers that were considerably small and hard to read values from.
Voltage measurements taken from the oscilloscope are subject to various forms of noise.This noise was removed as far as possible by eliminating ground loops using electromagnetic shielding and decoupling capacitors.
The number of measurements taken was increased using the averaging mode of the oscilloscope.This function displays the average of the last 16 waveforms.While this removes random noise in the chip and measurement circuit, it hides the true noise performance of the IHC circuit.
Change in response to temperature variations was not measured though there may have been some error introduced in the results due to variation in temperature during the experiments.This was reduced to a minimal level by leaving the chip turned on continuously over the days that the experiments were performed.

Rate-intensity functions
Rate-intensity functions are plots of firing rate response versus stimulus intensity and indicate the dynamic range of the model.The method of Smith and Zwislocki [11] is used to find the rate-intensity functions.Firstly, a stimulus level is found where the onset and steady-state rates are zero.This zero-dB level is the reference level.Responses are recorded for 300 ms tone bursts in steps of 10 dB to 40 dB above the reference.The rise time of the signal and the duration of the recording interval are the same as those used by Meddis, as these parameters affect the shape of the onset rate-intensity function [10].
Three rates are identified in the response, shown in Figure 7.The spontaneous rate represents the fibre response in the absence of stimuli.The onset rate is the firing rate averaged over the first 1 ms while the steady-state rate is the response averaged over the last 30 ms of a 200 ms tone burst.The rates, plotted against stimulus level, were measured directly from the oscilloscope traces using a constant gain h to convert the output voltage to a rate value for comparison with biological data (Figure 8).The onset rate is seen to increase monotonically with stimulus level and shows little or no sign of saturation.The steady-state rate is independent of stimulus level (straight line).These results agree with those reported by Meddis and with physiological results.

Recovery after masking
Tone bursts can be masked by preceding tones, depending on how the hair cell recovers after adapting to the masking tone.It has been established that auditory nerve recovery from masking stimuli follows a single exponential curve,  where the response at stimulus onset recovers at a faster rate than the total response [12,13].
The method of Westerman [14] is used in this test with 43 dB tone bursts at 1 kHz.Firstly, an unadapted response is measured in the absence of a masking tone.Then the response to a probe following a masking tone is measured as a decrement from the unadapted response (Figure 9).This is repeated for probes with increasing time delay.The onset  rate is measured within the first 1 ms and the steady-state rate after 20 ms.The masker has duration of 300 ms and the probe 30 ms.The time delay is varied between 0 and 200 ms. Figure 9 shows the chip replicating the response of the Meddis model in forward masking.

Additivity
A model is additive if changes in the firing rate caused by increases or decreases in stimulus levels are independent of the state of adaptation.Short-term and rapid adaptations have been shown to be additive in the IHC [11,15].This test uses the method by Smith [12] which is designed to emphasize the properties of rapid adaptation.This method uses short (SW) and long (LW) analysis windows as shown in Figure 10.Increments and decrements of 6 dB are applied at various delays, ∆t from the start of the pedestal.The control response is a pedestal with no increment nor decrement.For each window, the increase or decrease in firing rate from the control response is measured.Smith found that adaptation was additive in the short term (Large windows of 20 ms) for both increments and decrements.Rapid adaptation (small windows of 1 ms) was found to be additive for level increments, while decrements decreased the short-   term firing rate with increasing time delay, and in proportion to the decrease in firing rate.
Figure 11 shows that in the Meddis model, and hence in the chip, increments in the short term are not additive.This error is thought to be due to the small number of reservoirs used in the Meddis model [10].Models that use multiple-reservoir sites were shown to report adaptation trends correctly, with the penalty of decreased computational efficiency.Multiple reservoir models (e.g., [16]) contain multiple release sites that are spatially ordered by increasing threshold.This attribute gives these models a timeindependent response to time-varying stimuli.However, the results from rapid increments agree with the findings of Smith.Furthermore, the measurements of rapid and shortterm decrements (Figure 12) also agree with Smith's findings.

Two-component adaptation
The adaptation curve was characterised by Smith and Westerman [14] as the sum of two exponentially decaying components (t r and t st ) (Figure 13): where t r is the decay time constant of rapid adaptation and t st is the decay time constant of short-term adaptation.The  magnitudes of the two components are given by A r and A st , respectively, and A ss is the steady-state response.
In the Meddis model, t r is largely determined by the time constant parameters τ c = 1/(l+r) associated with the cleft filter (model parameter = 2 ms).In the literature, it is reported to be between 1 and 10 ms and decreases with increases in stimulus level [14,17].The decay time constant of shortterm adaptation t st is largely determined by the time constant τ y = 1/y associated with the transmitter factory (model parameter = 50-200 ms).In the literature, it is reported to be between 20-100 ms and is independent of tone level [11,14].The third time constant of the Meddis model is the time constant of the store reservoir (model parameter = 1 ms) and it represents a lowpass filter of around 1 kHz [18], which was found to be related to phase locking which is seen in the next test.
The time constants were derived from the model response to 100 ms tone bursts varying in amplitude from 10 dB to 40 dB.For the rapid time constant, points at 1 ms and 2 ms were measured and for the short-term time constant, points at 40 ms and 80 ms were measured on the response.The steady-state response A ss was also measured as the level towards the end of adaptation at 95 ms.The time constants were calculated from this data using the straight line approximation technique reported in [9].Again the results of the chip compare well to that of the Meddis model and the data of Westerman.The rapid adaptation component is independent of the stimulus level whereas the shortterm adaptation component decreases with stimulus level (Figure 14).

Low-frequency phase locking
At low frequencies the auditory nerve response synchronises or "phase locks" with the positive half cycle of the stimulus tone.This is shown by the high magnitude AC (alternating current) component in the output that resembles the stimulus signal phase characteristics (Figure 15a).At frequencies above 1 kHz there is no longer a phase lock and the AC component is removed (Figure 15b).This corresponds to the   store filter in the Meddis model which has a time constant τ w of 1 ms.The synchronization index (SI) used by Meddis is not used here as it shows that the model has a first-order lowpass filter response above 1 kHz and this filter has already been identified [18].(The graphs in [10] show that the SI or the gain drops by 5 dB in half a decade.)

Recovery of spontaneous activity & computational efficiency
After a tone burst, the auditory nerve firing rate ceases before recovering to the spontaneous rate (Figure 16).An exponential with a time constant between 20 ms and 100 ms describes the recovery function [12,13,14].The measurements from the chip showed a recovery rate that compares well to that reported by Meddis and physiological data (Table 1).Indeed, the result is better than that found in simulation, which suggests that a closer parameter set was found using the chip.
Computational efficiency is determined by the time to process a 1 second tone.As IHC models become increasingly popular in speech processing systems, the computational efficiency is an important metric.While the computational models typically took 2-5 seconds to simulate on a computer, the chip was able to respond in real time as a real IHC would.This improvement in response time enabled a better parameter set to be found quickly.The speed of this model will be useful for physiological investigators in understanding more about the auditory systems and in the construction of sensors based on auditory physiology.

CONCLUSIONS
An improved analogue VLSI implementation of the IHC is presented, utilising the Meddis model, a widely accepted and computationally convenient model.This model is transferred to the current domain and is constructed using translinear and log-domain circuits.A chip was fabricated and was successfully tested against seven different models of IHC functions.The chip was seen to be faster than the computer simulation (i.e., it worked in real time).As a result, it was easier to find a parameter set due to the hands on tuning provided by this analogue chip.

Figure 3 :
Figure 3: The Meddis model in the current mode.

Figure 7 :
Figure 7: Plot of rate intensity functions.

Figure 8 :
Figure 8: Oscilloscope traces of IHC output for various stimulus levels.
Response to 5 kHz tone burst.

Table 1 :
Recovery time constant and computational efficiency.