Model of cochlear membrane adapted (Peterson, Bogert 1950)

 Download zip file 
Help downloading and running models
This model, adapted from Peterson and Bogert (1950), simulates the response of the gerbil basilar membrane to a pure tone stimulus. This model does not attempt to simulate the effect of outer hair cell motility. The program prompts the user for the stimulus frequency and the Q (quality factor) for the basilar membrane impedance. It then plots cochlear partition volume velocity, the pressure difference across the partition and the cochlear partition impedance as a function of cochlear location. More information on the actual computations are contained in comments within the m-file.
1 . Peterson LC, Bogert BP (1950) A dynamical theory of the cochlea. J Acoust Soc Am 22:369-381
Model Information (Click on a link to find other models with that property)
Model Type:
Brain Region(s)/Organism:
Cell Type(s):
Gap Junctions:
Simulation Environment: MATLAB;
Model Concept(s):
Implementer(s): Naidu, R ; Mountain, David [dcm at];
% Pbmodel1.m
%   		1D Linear Cochlear Model for the gerbil
%					R. Naidu and D. Mountain
%				Cochlear Biophysics Laboratory
%			Boston University Hearing Research Center
% The model computes Pressure difference and the BM volume
% velocity by picking a number for the fluid volume velocity
% through the helicotrema and then working back towards
% the base.  At each section the fluid volume velocity is computed
% by adding the fluid volume velocity from the previous section to
% the basilar membrane volume velocity.  The pressure drop across the
% fluid due to this velocity is then computed and added to the pressure
% difference across the cochlear partition.  Since the model is linear,
% the results above can then be scaled to the pressure difference at the
% stapes.

%The following are for gerbil, other frequency-place maps are available at:

L   = 1.21;		% length of cochlea (cm)
a   = 398.;		%Greenwood coefficient
b   = 2.2;		%Greenwood coefficient
k   = 0.631;	%Greenwood coefficient
N = 512; 		% No. of sections
T = L/N; 		% spatial step size

rho = 1.02;		% Fluid density (gm/cm^3)

P  = zeros(1,N);	% Pressure difference array (dyne/cm)
Zp = zeros(1,N);	% Cochlear Partition Impedance array (cgs Ohms)
Vbm= zeros(1,N);	% Basilar membrane volume velocity array (cm^3/s)
						% To get to mechanical velocity (cm/s), this needs to be mulitiplied
						% by 1/(T*Wbm) where Wbm is the effective width of the BM

f = input('Frequency of tone ==>');
q = input('Quality factor (try 10) ==>');

w = 2*pi*f;

for i = 1:N					% Compute model parameters
	Wbm(i) = 0.0175;		%Width of the BM (cm)
	y(i)   = 0.01;			%Thickness of the organ of Corti (cm)
	As(i)  = .0002;		%Cross-sectional area of scalae (cm^2)
	Q(i)   = q;				%Quality factor of the RLC circuit
	x(i)   = i*T;			%Distance from base (cm)

	Ls(i)  = rho*T/As(i);						%Mass of scalae
	CF(i)  = a*(10^(((L-x(i))/L)*b) - k);	%Frequency-place map
	Cbm(i) =  exp(-30.9 + 4.04*x(i));		%Volume compliance of cochlear partition
	Lbm(i) = 1/(Cbm(i)*(2*pi*CF(i))^2);		%Effective Mass of cochlear partition
	Rbm(i) = (1/Q(i))*(Lbm(i)/Cbm(i))^0.5;	%viscous resistance of cochlear partition
%Cochlear partition impedance
	Zp(i)  = (-w^2*Lbm(i)*Cbm(i) + 1 + Rbm(i)*j*w*Cbm(i))/(j*w*Cbm(i));   %j = sqrt(-1)
% define the conditions at the helicotrema modeling it as a resistance:
Rh  = (Ls(N)/Cbm(N))^0.5;		%Match the cochlear impedance	
Vh  = 1;								%Pick a velocity at the helicotrema
Ph = Rh*Vh;
for i = N:-1:1						% Start from the apex and solve for P and vBM
    if i == N
       P(i) = Ph;    			% Boundary conditions at helicotrema
       V(i) = Vh + P(i)/Zp(i);
       P(i) = P(i+1) + V(i+1)*j*w*Ls(i+1);
       V(i) = V(i+1) + P(i)/Zp(i);

    Vbm(i)= P(i)/Zp(i);
norm_Vbm = abs(Vbm)/abs(P(1));	% Scale for input pressure of 1 dyne/cm^2 inside stapes
norm_P  = abs(P)/abs(P(1));		% which corresponds to approximately 40 dB SPL
norm_Zp = abs(Zp);


MagFig=figure;				%Start generating the magnitude figure (assume 1024x768 display
set(MagFig, 'Position',[5  200   500   500])	%Position it on the left side of screen

semilogy(x, norm_Vbm);
axis([0 1.4 1e-10 1e-6]);
ylabel('BM Volume Velocity');
title('Magnitude (cgs units)');

semilogy(x, norm_P);
axis([0 1.4 0.01 1]);
ylabel('Pressure Difference');

semilogy(x, norm_Zp);
axis([0 1.4 1e6 1e10]);
xlabel('Distance from Base (cm)');
ylabel('CP Impedance');

PhaseFig=figure;				%Start generating the magnitude figure (assume 1024x768 display
set(PhaseFig, 'Position',[512  200   500   500])	%Position it on the left side of screen

plot(x, theta_Vbm);
%axis([0 1.4 -5 5]);
ylabel('BM Velocity');
title('Phase (radians)');

plot(x, theta_P);
%axis([0 1.4 -4 0 ]);
ylabel('Pressure Difference');

plot(x, theta_Zp);
axis([0 1.4 -2 2]);
xlabel('Distance from Base (cm)');
ylabel('CP Impedance');

Loading data, please wait...