Human auditory periphery model: cochlea, IHC-AN, auditory brainstem responses (Verhulst et al 2018)

 Download zip file 
Help downloading and running models
Accession:246535
The human auditory periphery model can simulate single-unit response of basilar-membrane vibration, IHC receptor potential, instantaneous AN/CN/IC firing rates, as well as population responses such as otoacoustic emissions, auditory brainstem responses. The neuron models (IHC, AN,CN,IC) can be run independently to relate their responses to electrophysiology, or be simulated as part of the human auditory periphery.
References:
1 . Verhulst S, Altoè A, Vasilkov V (2018) Computational modeling of the human auditory periphery: Auditory-nerve responses, evoked potentials and hearing loss. Hear Res 360:55-75 [PubMed]
2 . Altoè A, Pulkki V, Verhulst S (2018) The effects of the activation of the inner-hair-cell basolateral K+ channels on auditory nerve responses. Hear Res 364:68-80 [PubMed]
3 . Altoè A, Pulkki V, Verhulst S (2014) Transmission line cochlear models: improved accuracy and efficiency. J Acoust Soc Am 136:EL302-8 [PubMed]
4 . Verhulst S, Dau T, Shera CA (2012) Nonlinear time-domain cochlear model for transient stimulation and human otoacoustic emission. J Acoust Soc Am 132:3842-8 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type:
Brain Region(s)/Organism:
Cell Type(s): Cochlea hair inner GLU cell; Cochlear nucleus bushy GLU cell; Auditory nerve; Brainstem neuron;
Channel(s):
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: Python; MATLAB;
Model Concept(s):
Implementer(s): Verhulst, Sarah [s.verhulst at ugent.be]; Altoé, Alessandro ;
Search NeuronDB for information about:  Cochlear nucleus bushy GLU cell; Cochlea hair inner GLU cell;
    # -*- coding: utf-8 -*-
import numpy as np
import time
from scipy.integrate import ode
from scipy import signal
import ctypes
import os
import sys

DOUBLE = ctypes.c_double
INT = ctypes.c_int
PINT = ctypes.POINTER(ctypes.c_int)
PLONG = ctypes.POINTER(ctypes.c_long)
PDOUBLE = ctypes.POINTER(ctypes.c_double)


class tridiag_matrix(ctypes.Structure):
    _fields_ = [("aa", ctypes.POINTER(ctypes.c_double)),
                ("bb", ctypes.POINTER(ctypes.c_double)),
                ("cc", ctypes.POINTER(ctypes.c_double))]

# load C library
libName = 'tridiag.dll' if 'win' is sys.platform else 'tridiag.so'
libtrisolv = np.ctypeslib.load_library(libName,os.path.dirname(os.path.abspath(__file__)))

# load tridiagonal solver function and defines input
libtrisolv.solve_tridiagonal.restype = None
libtrisolv.solve_tridiagonal.argtypes = [ctypes.POINTER(tridiag_matrix),  # aa
                                         PDOUBLE,  # vv
                                         PDOUBLE,  # solution
                                         INT,  # nrows
                                         ]

libtrisolv.delay_line.restype = None  # TODO SPEEDUP W POINTERS!
libtrisolv.delay_line.argtypes = [PDOUBLE,  # in_matrix
                                  PINT,  # delay1
                                  PINT,  # delay2
                                  PINT,  # delay1
                                  PINT,  # delay1
                                  PDOUBLE,  # dev
                                  PDOUBLE,  # YZweig
                                  INT,  # delay_buffer_length
                                  INT  # n
                                  ]

# definition of the function


def TLsolver(t, y, model):  # y''=dv/dt y'=v
    n = model.n + 1
    frac = (t - model.lastT) / model.dt
    a = model.interplPoint1
    b = model.interplPoint2
    c = model.interplPoint3
    d = model.interplPoint4
    cminusb = c - b
    #fast cubic interpolation
    F0 = b + frac * \
        (cminusb - 0.1666667 * (1. - frac) *
         ((d - a - 3.0 * cminusb) * frac + (d + 2.0 * a - 3.0 * b)))
    model.Vtmp = y[0:n]
    model.Ytmp = y[n:2 * n]
    if(model.non_linearity):  # non-linearities here
        factor = 100
        Vvect = np.abs(model.Vtmp) / model.RthV1
        Sxp = (Vvect - 1.) * model.const_nl1
        Syp = model.Sb * np.sqrt(1 + (Sxp / model.Sa) ** 2)
        Sy = Sxp * model.sinTheta + Syp * model.cosTheta
        SheraP = model.PoleS + Sy / factor
        SheraP = np.fmin(SheraP, model.PoleE)
        # update non-linear parameters here only if the pole displacement is
        # larger than 1%
        if(np.max(abs(SheraP[1:n] - model.SheraP[1:n]) /
           abs(model.SheraP[1:n])) > 0.01):
            model.SheraP = SheraP
            model.SheraParameters()
            model.ZweigImpedance()
            model.current_t = t
    model.Dev[0:n] = model.Dev[0:n] + frac
    libtrisolv.delay_line(
        model.Ybuffer_pointer, model.Zrp_pointer, model.Zrp1_pointer,
        model.Zrp2_pointer, model.Zrp3_pointer, model.Dev_pointer,
        model.YZweig_pointer, ctypes.c_int(model.YbufferLgt),
        ctypes.c_int(model.n + 1))
    model.Dev[0:n] = model.Dev[0:n] - frac
    model.calculate_g()
    model.calculate_right(F0)
    #compute q
    libtrisolv.solve_tridiagonal(
        ctypes.byref(model.tridata), model.r_pointer, model.Qpointer,
        ctypes.c_int(n))
    zero_val = (model.RK4_0*model.Qsol[0] + model.RK4G_0*(model.g[0] + model.p0x * F0))

    Vderivative = (model.Qsol - model.g)
    Vderivative[0] = zero_val;
    solution = np.concatenate([Vderivative, model.Vtmp])
    return solution


class cochlea_model ():

    # init constants
    def __init__(self):
        self.ttridiag = 0
        self.calling_function = 0
        self.cochleaLength = .035
        self.bmMass = 0.5
        self.bmImpedanceFactor = 1
        self.scalaWidth = 0.001
        self.scalaHeight = 0.001
        self.helicotremaWidth = 0.001
        self.rho = float(1e3)
        self.Normal_Q = 20
        self.Greenwood_A = 20682
        self.Greenwood_alpha = 61.765
        self.Greenwood_B = 140.6
        self.stapesArea = float(3e-6)
        self.EardrumArea = float(60e-6)
        self.MiddleEarResonanceFrequency = float(2e3)
        self.MiddleEarQualityFactor = 0.4
        self.SpecificAcousticImpedanceOfAir = 415
        self.middleEarTransformer = 30
        self.damping_coupler = float(140e5)
        self.mass_coupler = float(43.2e2)
        self.stiffness_coupler = 1. / float(2.28e-11)
        self.p0 = float(2e-5)
        self.ZweigQ = 1 / 0.0606
        self.ZweigFactor = 1.7435
        self.ZweigQAtBoundaries = 20
        self.ZweigBeta = 10000
        self.ZweigGamma = 6200
        self.ZweigN = 1.5
        self.SheraMuMax = 3
        self.RMSref = 0.6124
        self.Rme = float(0.3045192500000000e12)  # TODO setRme function
        #variable to check if the model is intialize before calling the solver
        self._is_init = 0
        self.interplPoint1 = 0
        self.interplPoint2 = 0
        self.interplPoint3 = 0
        self.interplPoint4 = 0

# function to intitialize all the parameters
    def init_model(self, stim, samplerate, sections, probe_freq, sheraPo,
                   compression_slope=0.4, Zweig_irregularities=1,
                   non_linearity_type='vel', KneeVar=1.,
                   low_freq_irregularities=1, subject=1,IrrPct=0.05):
        self.low_freq_irregularities = low_freq_irregularities
        self.SheraPo = np.zeros(sections+1)
        self.SheraPo = self.SheraPo+sheraPo  # can be vector or single value, line changed so it can work with both single and vector
        self.KneeVar = (KneeVar)
        self.IrrPct = IrrPct
        self.non_linearity = 0
        self.use_Zweig = 1
        if(Zweig_irregularities == 0):
            self.use_Zweig = 0
        if(non_linearity_type == 'disp'):
            self.non_linearity = 1
        elif(non_linearity_type == 'vel'):
            self.non_linearity = 2
        else:
            self.non_linearity = 0  # linear model
        self.n = sections
        self.fs = samplerate
        self.dt = 1. / self.fs
        self.probe_freq = probe_freq
        self.initCochlea()
        self.initMiddleEar()
        self.SetDampingAndStiffness()
        self.initZweig()
        self.initGaussianElimination()
        self.compression_slope_param(compression_slope)
        self.is_init = 1
        self.lastT = 0
        self.seed = subject  # change here the seed
        np.random.RandomState(self.seed)
        np.random.seed(self.seed)
        self.Rth = 2 * (np.random.random(self.n + 1) - 0.5)
        self.Rth_norm = 10 ** (self.Rth / 20. / self.KneeVar)
        lf_limit = self.ctr
        if(self.use_Zweig==0):
            lf_limit=0
            print('No irregularities')
        factor = 100
        n = self.n + 1
        Rth = self.Rth
        Rth_norm = self.Rth_norm
        #Normalized RTH, so save a bit of computation
        self.RthY1 = self.Yknee1 * Rth_norm
        self.RthY2 = self.Yknee2 * Rth_norm
        self.RthV1 = self.Vknee1 * Rth_norm
        self.RthV2 = self.Vknee2 * Rth_norm

        Rndm = self.IrrPct * Rth / 2.
        self.PoleS = (1 + Rndm) * self.SheraPo

        self.RthY1[lf_limit:n] = self.Yknee1
        self.RthY2[lf_limit:n] = self.Yknee2
        self.RthV1[lf_limit:n] = self.Vknee1
        self.RthV2[lf_limit:n] = self.Vknee2
        self.PoleS[lf_limit:n] = self.SheraPo[lf_limit:n]
        Theta0 = np.arctan(
            ((self.PoleE - self.PoleS) * factor) /
            ((self.RthV2 / self.RthV1) - 1.))
        Theta = Theta0 / 2.
        Sfoc = (self.PoleS * factor) / (self.RthV2 / self.RthV1)
        Se = np.cos((np.pi - Theta0) * 0.5)
        self.Sb = Sfoc * Se
        self.Sa = Sfoc * np.sqrt(1. - ((Se ** 2)))
        self.const_nl1 = np.cos(Theta) / np.cos(2 * Theta)
        self.cosTheta = np.cos(Theta)
        self.sinTheta = np.sin(Theta)

        #
        # PURIAM1 FILTER             ###
        #
        puria_gain = 10 ** (18. / 20.)
#         was the orignal Puria in 2012
#        second order butterworth
#        b, a = signal.butter(
#           1, [100. / (samplerate / 2.), 3000. / (samplerate / 2)],
#            'bandpass')
#         self.stim = signal.lfilter(b * puria_gain, a, stim)


        b, a = signal.butter(1, [600 / (samplerate / 2.), 4000. / (samplerate / 2)],'bandpass')
        self.stim = signal.lfilter(b * puria_gain, a, stim)

    # from intializeCochlea.f90
    def initCochlea(self):
        self.bm_length = self.cochleaLength - self.helicotremaWidth
        self.bm_width = self.scalaWidth
        self.bm_mass = self.bmMass * self.bmImpedanceFactor
        self.ZweigMso = 2. * self.rho / (self.bm_width * self.scalaHeight)
        self.ZweigL = 1. / (2.3030 * self.Greenwood_alpha)
        self.ZweigOmega_co = 2.0 * np.pi * self.Greenwood_A-self.Greenwood_B
        self.ZweigMpo = (
            self.ZweigMso * (self.ZweigL ** 2)) / ((4 * self.ZweigN) ** 2)
        self.Ko = self.ZweigMpo * (self.ZweigOmega_co ** 2)
        self.x =np.array(np.linspace(0, self.bm_length, self.n+1), order='C') #Old
        self.dx = self.bm_length / (1. * self.n)
        self.g = np.zeros_like(self.x)
        self.Vtmp = np.zeros_like(self.x)
        self.Ytmp = np.zeros_like(self.x)
        self.Atmp= np.zeros_like(self.x)
        self.right = np.zeros_like(self.x)
        self.r_pointer = self.right.ctypes.data_as(PDOUBLE)
        self.zerosdummy = np.zeros_like(self.x)
        self.gamma = np.zeros_like(self.x)
        self.Qsol = np.zeros_like(self.x)
        self.Qpointer = self.Qsol.ctypes.data_as(PDOUBLE)

    def initMiddleEar(self):
        self.q0_factor = self.ZweigMpo * self.bm_width
        self.p0x = self.ZweigMso * self.dx/(1. * self.ZweigMpo * self.bm_width)
        self.d_m_factor = -self.p0x * self.stapesArea * self.Rme
        self.RK4_0 = -(self.bm_width * self.ZweigMpo) / (self.stapesArea)
        self.RK4G_0 = (self.ZweigMpo * self.bm_width) / (
            self.ZweigMso * self.stapesArea * self.dx)

    def SetDampingAndStiffness(self):
        self.f_resonance = self.Greenwood_A * \
            10 ** (-self.Greenwood_alpha * self.x) - self.Greenwood_B
        self.ctr = np.argmin(np.abs(self.f_resonance - 100.))
        if(self.low_freq_irregularities):
            self.ctr = self.n + 1
        self.onek = np.argmin(np.abs(self.f_resonance - 1000.))
        self.omega = 2. * np.pi * self.f_resonance
        self.omega[0]=self.ZweigOmega_co
        self.omega2 = self.omega ** 2
        self.Sherad_factor = np.array(self.omega)
        self.SheraP = np.zeros_like(self.x)
        self.SheraD = np.zeros_like(self.x)
        self.SheraRho = np.zeros_like(self.x)
        self.SheraMu = np.zeros_like(self.x)
        self.SheraP = self.SheraPo + self.SheraP
        self.c = 120.8998691636393

        #
        # PROBE POINTS               ##
        #
        if(self.probe_freq=='all'):
            self.probe_points=np.zeros(len(self.f_resonance)-1)
            for i in range(len(self.f_resonance)-1):
                self.probe_points[i]=i+1
            self.probe_points=(self.probe_points)
            self.cf=(self.f_resonance[1:len(self.f_resonance)])
        elif(self.probe_freq=='half'):
            self.probe_points=np.zeros((len(self.f_resonance)-1)/2)
            for i in range((len(self.f_resonance)-1)/2):
                self.probe_points[i]=i+1
            self.probe_points=(self.probe_points)
            self.cf=(self.f_resonance[range(1,len(self.f_resonance),2)])
        elif(self.probe_freq=='abr'):
            #            self.probe_points=np.zeros([401,1])
            self.probe_points=np.array(range(110,911,2))
            self.cf=(self.f_resonance[range(110,911,2)])
#            print(np.shape(self.probe_points))
#            print('abr')
        else:
            self.probe_points=np.zeros(self.probe_freq.size,dtype=int)
            for i in range(len(self.probe_freq)):
                idx_help=abs((self.f_resonance)-np.float(self.probe_freq[i]))
                self.probe_points[i]=np.argmin(idx_help)
            self.cf=self.f_resonance[self.probe_points]
        self.probe_points=np.array(self.probe_points)
#        print(np.shape(self.probe_points))


    def initZweig(self):
        n = self.n + 1
        self.exact_delay = self.SheraMuMax / (self.f_resonance * self.dt)
        self.delay = np.floor(self.exact_delay) + 1
        self.YbufferLgt = int(np.amax(self.delay)) #quick fix, maybe can be shorter but it works so better not touch
        self.Ybuffer = np.zeros([n, self.YbufferLgt])
                                # Ybuffer implemented here as a dense matrix
                                # (python for cycles are slow...)
        self.Ybuffer = np.array(self.Ybuffer, order='C', ndmin=2, dtype=float)
        self.Ybuffer_pointer = self.Ybuffer.ctypes.data_as(PDOUBLE)
        self.ZweigSample1 = np.zeros_like(self.exact_delay)
        self.Zwp = int(0)
        self.ZweigSample1[0] = 1.
        self.ZweigSample2 = self.ZweigSample1 + 1
        # init buffers etc...
        self.Dev = np.zeros_like(self.x)
        self.Dev_pointer = self.Dev.ctypes.data_as(PDOUBLE)
        self.YZweig = np.zeros_like(self.x)
        self.YZweig_pointer = self.YZweig.ctypes.data_as(PDOUBLE)
        self.Zrp = np.array(np.zeros(n), dtype=np.int32, order='C')
        self.Zrp_pointer = self.Zrp.ctypes.data_as(PINT)
        self.Zrp1 = np.array(np.zeros(n), dtype=np.int32, order='C')
        self.Zrp1_pointer = self.Zrp1.ctypes.data_as(PINT)
        self.Zrp2 = np.array(np.zeros(n), dtype=np.int32, order='C')
        self.Zrp2_pointer = self.Zrp2.ctypes.data_as(PINT)
        self.Zrp3 = np.array(np.zeros(n), dtype=np.int32, order='C')
        self.Zrp3_pointer = self.Zrp3.ctypes.data_as(PINT)

    #set tridiagonal matrix values for transmission line
    def initGaussianElimination(self):
        n = self.n + 1
        self.ZweigMs = (self.ZweigMso * self.ZweigOmega_co) / self.omega # TAPERING self.omega*
        self.ZweigMp = self.Ko / (self.ZweigOmega_co * self.omega)
        #self.ZweigMs = self.ZweigMs/ self.omega[1]
        #self.ZweigMp = self.ZweigMp/self.omega[1]
        self.ZASQ = np.zeros_like(self.x)
        self.ZASC = np.zeros_like(self.x)
        self.ZAL = np.zeros_like(self.x)
        self.ZAH = np.zeros_like(self.x)
        # init values of transimission line
        self.ZASQ[0] = 1.
        self.ZASC[0] = 1 + self.ZweigMso * self.dx
        self.ZAH[0] = -1*self.ZweigOmega_co/self.omega[1]
        self.ZAL[1:n] = -self.ZweigMs[1:n]*self.omega[1:n]/self.omega[0:n-1]
        self.ZAH[1:n-1] =-self.ZweigMs[0:n-2]*self.omega[1:n-1]/self.omega[2:n]
#        self.ZAH[0] = -1.
#        self.ZAL[1:n] = -self.ZweigMs[1:n]
#        self.ZAH[1:n - 1] = -self.ZweigMs[0:n - 2]
        self.ZASQ[1:n] = self.omega[1:n] * self.ZweigMs[1:n] * self.ZweigMs[0:n - 1] * (self.dx ** 2) / (self.ZweigOmega_co *self.ZweigMpo)
        self.ZASC[1:n] = self.ZASQ[1:n] +self.ZweigMs[1:n] + self.ZweigMs[0:n - 1]
        self.tridata = tridiag_matrix()
        self.tridata.aa = self.ZAL.ctypes.data_as(PDOUBLE)
        self.tridata.bb = self.ZASC.ctypes.data_as(PDOUBLE)
        self.tridata.cc = self.ZAH.ctypes.data_as(PDOUBLE)

    def calculate_g(self):  # same as in fortran
        n = self.n + 1
        self.g[0] = self.d_m_factor * self.Vtmp[0]
        dtot = self.Sherad_factor * self.SheraD
        stot = (self.omega2) * (self.Ytmp + (self.SheraRho * self.YZweig))
        self.g[1:n] = (dtot[1:n] * self.Vtmp[1:n]) + stot[1:n]

    def calculate_right(self, F0):  # same as in fortran
        n = self.n + 1
        self.right[0] = self.g[0] + self.p0x * F0
        self.right[1:n] = self.ZASQ[1:n] * self.g[1:n]

    def SheraParameters(self):  # same as in fortran
        a = (self.SheraP + np.sqrt((self.SheraP ** 2.) +
             self.c * (1.0 - self.SheraP ** 2))) / self.c
        self.SheraD = 2.0 * (self.SheraP - a)
        self.SheraMu = 1. / (2.*np.pi*a)
        self.SheraRho = 2. * a * \
            np.sqrt(1. - (self.SheraD / 2.) ** 2.) * np.exp(-self.SheraP / a)

    def ZweigImpedance(self):
        n = self.n + 1
        MudelayExact =(2*np.pi)*self.SheraMu / (self.omega * self.dt)
        Mudelay = np.floor(MudelayExact) + 1.
        self.Dev[:] = Mudelay - MudelayExact
        self.Zrp1[0:n] = (
            (self.Zwp + self.YbufferLgt) - Mudelay[0:n]) % self.YbufferLgt
        const = self.YbufferLgt - 1
        self.Zrp[0:n] = (self.Zrp1[0:n] + const) % self.YbufferLgt
        self.Zrp2[0:n] = (self.Zrp1[0:n] + 1) % self.YbufferLgt
        self.Zrp3[0:n] = (self.Zrp2[0:n] + 1) % self.YbufferLgt

    def compression_slope_param(self, slope):
        self.Yknee1 = float(1.0*(6.9183e-10))
#self.Yknee1 = float(2.0*(6.9183e-10))
        self.Yknee2 = float(1.5488e-8)
    #    #       THdB=10.0 #SARAH's Style
#        THdB=10.0
#        Ax=1
#        Bx=20*np.log10(8.461e-11)
#        Vknee1 = float(1.0*(2.293e-7))
#        self.PoleE = np.zeros_like(self.x)+0.3
#        BoffsetV=-slope*THdB+20*np.log10(Vknee1) #find the offset of the compression curves
#        Vint=(Bx-BoffsetV)/(slope-Ax) #is the intersection in dB on xaxis
#        Vknee2=Ax*Vint+Bx #what it corresponds to in dB on y axis
#        self.Vknee2=10.0 ** (Vknee2/20.0)
#        self.Vknee1=Vknee1

#       Ale Style
        self.PoleE = np.zeros_like(self.x)+0.31 #saturating pole
        v1=0.6807e-08/3/np.sqrt(2); # velocity at -10 dB with starting Pole
        v2=26.490e-11/3/np.sqrt(2); # velocity at -10 dB with saturating pole
        K1dB=20; # Knee point of the first linear regime in dB (you can select it from here now)
        #but it does not work precisely for K1dB<20...?? So, by using v1 and v2 peak velocities at -10 dB it is possible to impose the desired Knee down to 10 dB.
        K1dB=K1dB+20; #fix for the -10 dB of v1 and v2
        K1L=10**(K1dB/20) #knee point in linear scale
        self.Vknee1=K1L*v1
        vst1dB=20*np.log10(v1)+K1dB #velocity with the two poles when the compression starts
        vst2dB=20*np.log10(v2)+K1dB
        K2dB=(vst1dB-vst2dB)/(1-slope) #intersection in dB re Knee 1
        self.Vknee2=v2*10**(K2dB/20)*K1L
    
    def polecalculation(self):  # TODO

        factor = 100.
        # lf_limit = self.ctr
        # n = self.n + 1
        if(self.non_linearity == 1):  # To check
            # non-linearity DISP cost about three times more than in
            # fortran (Not implemented now)
            Yknee1CST = self.RthY1 * self.omega[self.onek]
            Yknee2CST = self.RthY2 * self.omega[self.onek]
            Yknee1F = Yknee1CST / self.omega
            Yknee2F = Yknee2CST / self.omega
            Yvect = np.abs(self.Ytmp / Yknee1F)
            Theta0 = np.arctan(
                ((self.PoleE - self.PoleS) / ((Yknee2F / Yknee1F) - 1.)))
            Theta = Theta0 / 2.
            # save 2 call to trigonometric function on vector by storing
            # some data
            cos_Theta = np.cos(Theta)
            sin_Theta = np.sin(Theta)
            cos_Theta0 = 2 * cos_Theta ** 2 - 1
            Sfoc = self.PoleS * factor * (Yknee2F / Yknee1F)
            Se = sin_Theta
            Sb = Sfoc * Se
            Sa = Sfoc * np.sqrt(1. - (1. * (Se ** 2)))
            Sxp = (Yvect - 1.) * cos_Theta / cos_Theta0
            Syp = Sb * np.sqrt(1 + (Sxp / Sa) ** 2)
            Sy = Sxp * sin_Theta + Syp * cos_Theta
            self.SheraP = self.PoleS + Sy / factor

        elif(self.non_linearity == 2):  # non-linearity VEL
            Vvect = np.abs(self.Vtmp) / self.RthV1
            Sxp = (Vvect - 1.) * self.const_nl1
            Syp = self.Sb * np.sqrt(1 + (Sxp / self.Sa) ** 2)
            Sy = Sxp * self.sinTheta + Syp * self.cosTheta
            self.SheraP = self.PoleS + Sy / factor
        else:
            print('linear')
            self.SheraP = self.PoleS
        self.SheraP = np.fmin(self.SheraP, self.PoleE)

    def solve(self):
        n = self.n + 1
        tstart = time.time()
        if not(self.is_init):
            print("Error: model to be initialized")
        length = np.size(self.stim) - 2
        time_length = length * self.dt
        #each probe point signal in a row
        self.Vsolution = np.zeros([length + 2, len(self.probe_points)])
        self.Ysolution = np.zeros([length + 2, len(self.probe_points)])
        self.Asolution= np.zeros([length + 2, len(self.probe_points)])
        self.oto_emission = np.zeros(length + 2)
        self.time_axis = np.linspace(0, time_length, length)
        r = ode(TLsolver).set_integrator('dopri5', rtol=1e-2, atol=1e-13)
        r.set_f_params(self)
        r.set_initial_value(
            np.concatenate([np.zeros_like(self.x), np.zeros_like(self.x)]))
        r.t = 0
        j = 0
        self.last_t = 0.0
        self.current_t = r.t
        self.polecalculation()
        self.SheraParameters()
        self.ZweigImpedance()
        self.V1=np.zeros_like(self.x)
        while(j < length):
            if(j > 0):
                self.interplPoint1 = self.stim[j - 1]
            # assign the stimulus points and interpolation parameters
            self.interplPoint2 = self.stim[j]
            self.interplPoint3 = self.stim[j + 1]
            self.interplPoint4 = self.stim[j + 2]
            r.integrate(r.t + self.dt)
            self.lastT = r.t
            self.V1 = r.y[0:n]
            self.V1[0]=self.Vtmp[0]
            self.Y1 = r.y[n:2 * n]  # Non linearities HERE
            self.Atmp=self.Qsol-self.g
            self.Zwp = (self.Zwp + 1) % self.YbufferLgt  # update Zweig Buffer
            self.Ybuffer[:, self.Zwp] = self.Y1
            self.ZweigImpedance()
            self.current_t = r.t
            if(self.probe_freq=='all'):
                self.Vsolution[j,:] = self.V1[1:n]  #
                self.Ysolution[j,:] = self.Y1[1:n]
            elif(self.probe_freq=='half'):
                self.Vsolution[j,:]=self.V1[range(1,n,2)]
                self.Ysolution[j,:] = self.Y1[range(1,n,2)]
            elif(self.probe_freq=='abr'):
               self.Vsolution[j,:]=self.V1[range(110,911,2)]
               self.Ysolution[j,:] = self.Y1[range(110,911,2)]
            else:
                self.Vsolution[j,:] = self.V1[self.probe_points]  # storing the decided probe points
                self.Ysolution[j,:] = self.Y1[self.probe_points]
            self.oto_emission[j] = self.Qsol[0]
            j = j + 1
    # filter out the otoacoustic emission ####
        samplerate = self.fs
        b, a = signal.butter(1, [600 / (samplerate / 2.), 4000. / (samplerate / 2)],'bandpass')
        self.oto_emission = signal.lfilter(b * self.q0_factor, a, self.oto_emission)
        elapsed = time.time() - tstart
#        print(elapsed)
# END

Loading data, please wait...