Neural model of frog ventilatory rhythmogenesis (Horcholle-Bossavit and Quenet 2009)

 Download zip file 
Help downloading and running models
Accession:123987
"In the adult frog respiratory system, periods of rhythmic movements of the buccal floor are interspersed by lung ventilation episodes. The ventilatory activity results from the interaction of two hypothesized oscillators in the brainstem. Here, we model these oscillators with two coupled neural networks, whose co-activation results in the emergence of new dynamics. .. The biological interest of this formal model is illustrated by the persistence of the relevant dynamical features when perturbations are introduced in the model, i.e. dynamic noises and architecture modifications. The implementation of the networks with clock-driven continuous time neurones provides simulations with physiological time scales."
Reference:
1 . Horcholle-Bossavit G, Quenet B (2009) Neural model of frog ventilatory rhythmogenesis. Biosystems 97:35-43 [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism:
Cell Type(s):
Channel(s):
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: MATLAB;
Model Concept(s): Temporal Pattern Generation; Oscillations; Synchronization;
Implementer(s):
%Biosystems. 2009 Jul;97(1):35-43.
%Horcholle-Bossavit G, Quenet B.
%Neural model of frog ventilatory rhythmogenesis.
%Models from  Izhikevich E : http://nsi.edu/users/izhikevich/publications/spikes.htm
%Modified by Bruce Land -- BioNB 222, March 2008 

%this program performs the computation of interspike time P and mdeltat
format long
nNeuron = 2 ;
nCurrent = nNeuron;
dt = 0.125 ; %millisecond -- time step  for computation accuracy dt must be (1/2^n)
tEnd =1500; %maximum simulation time
time = dt:dt:tEnd ;
pars=[0.02      0.2     -65     8       14 ];    %1-tonic spiking
nType = 1*ones(1,nNeuron) ;
a=pars(nType,1)';
b=pars(nType,2)';
c=pars(nType,3)';
d=pars(nType,4)';
Is=0*ones(nNeuron,1)';

Idep=4.8;                               % steady depolarizing current input 
Gmaxsynex=50;                            %excitation weight 

depmes=5;                               %number of first spikes suppressed in the computation of P

w = [0   0;1 0] ;
SynStrength =Gmaxsynex*abs(w').* (w'>0); ...
R=zeros(1,nNeuron);
R(1)=1;
CurrentStrength = diag(R);
CurrentInput = zeros(length(time), nNeuron);
CurrentInput(:,1) =Idep* ones(size(CurrentInput(:,1)));
SynRatex = 1*ones(1,nNeuron);                       
SynDecayex = 1 - SynRatex*dt;
v = c;                                              % reset voltage;
s = zeros(1,nNeuron);                               % spike occurance
u = zeros(1,nNeuron);                               % revocery variable
Istate = 0;                                         %state variable used for synaptic currents decay
I = Is ;
tindex = 1;                                         %time pointer for output arrays
                                                    
Vout = zeros(length(time),nNeuron) ;
Init=zeros(nNeuron,1);
firings=[zeros(sum(Init),1),find(Init==1)];
for t=time
    I = Is + CurrentInput(tindex,:)*CurrentStrength ;
    v = v + dt * (0.04*v.^2+5*v+140-u+I);
    u = u + dt * a .* (b.*v-u);
    Vout(tindex,:) = v ;
    Iout(tindex,:) = I ;
    s = zeros(1,nNeuron) ;
    fired=find(v>=30);                                  % indices of cells that spiked
    if ~isempty(fired)
        firings=[firings; t+0*fired', fired'];
        v(fired) = c(fired) ;
        u(fired) = u(fired)+d(fired) ;
        
    end
tindex = tindex+1;                                      %update time index
end
indquand=find(firings(:,2)==1);
indquand=indquand(depmes:end);          %suppression of the first spikes
tempsfire1=firings(indquand,1);
interval1=diff(tempsfire1);
P=mean(interval1);                 %tonic response interspike
erreurP=std(interval1);
tindex = 1;                             %pointer for output arrays
Istate = 0;                             %state variable used for synaptic currents decay
I = Is ;
v = c;                                  % reset voltage;
s = zeros(1,nNeuron);                   % spike occurrence
u = zeros(1,nNeuron);                   % recovery variable

%init plotting variables
Vout = zeros(length(time),nNeuron) ;
Init=zeros(nNeuron,1);
firings=[zeros(sum(Init),1),find(Init==1)];
for t=time
    format short
    tempsavant=t-P;
    indqui=find(firings(:,1)<=tempsavant & firings(:,1)>(tempsavant-dt));
    qui=firings(indqui,2);
    s(qui)=1;
  Istate =  s*SynStrength+Istate .* SynDecayex.*(Istate>0);
  I = Istate +  Is + CurrentInput(tindex,:)*CurrentStrength ;
  v = v + dt * (0.04*v.^2+5*v+140-u+I);
    u = u + dt * a .* (b.*v-u);
    Vout(tindex,:) = v ;
    Iout(tindex,:) = I ;
    s = zeros(1,nNeuron) ;
    fired=find(v>=30);                      % indices of cells that spiked
    if ~isempty(fired)
        firings=[firings; t+0*fired', fired'];
        v(fired) = c(fired) ;
        u(fired) = u(fired)+d(fired) ;
    end;

    tindex = tindex+1;
end 
indquand2=find(firings(:,2)==2);
indquand2=indquand2(depmes:end);
tempsfire2=firings(indquand2,1);
indi=find(Iout(:,2)>=Gmaxsynex);
indi=indi(depmes:end)*dt;
if length(tempsfire2)==length(indi)
        deltat=tempsfire2-indi;
        Delex=mean(deltat);               %mean time delay between pre and post synaptic spikes 
        errdelex=std(deltat);
else
    disp('size problem for deltat')
    Delex=0;
end