Respiratory control model with brainstem CPG and sensory feedback (Diekman, Thomas, and Wilson 2017)

 Download zip file 
Help downloading and running models
Accession:229640
This is a closed-loop respiratory control model incorporating a central pattern generator (CPG), the Butera-Rinzel-Smith (BRS) model, together with lung mechanics, oxygen handling, and chemosensory components. The closed-loop system exhibits bistability of bursting and tonic spiking. Bursting corresponds to coexistence of eupnea-like breathing, with normal minute ventilation and blood oxygen level. Tonic spiking corresponds to a tachypnea-like state, with pathologically reduced minute ventilation and critically low blood oxygen. In our paper, we use the closed-loop system to demonstrate robustness to changes in metabolic demand, spontaneous autoresuscitation in response to hypoxia, and the distinct mechanisms that underlie rhythmogenesis in the intact control circuit vs. the isolated, open-loop CPG.
Reference:
1 . Diekman CO, Thomas PJ, Wilson CG (2017) Eupnea, tachypnea, and autoresuscitation in a closed-loop respiratory control model. J Neurophysiol 118:2194-2215 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Neuron or other electrically excitable cell;
Brain Region(s)/Organism: Brainstem;
Cell Type(s): Respiratory column neuron; PreBotzinger complex neuron;
Channel(s): I Na,p; I Na,t; I K;
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: MATLAB; XPP;
Model Concept(s): Pacemaking mechanism; Respiratory control;
Implementer(s): Diekman, Casey O. [casey.o.diekman at njit.edu];
Search NeuronDB for information about:  I Na,p; I Na,t; I K;
% This script reproduces Figure 3 by simulating the open-loop system for a 
% range of gtonic values and the closed-loop static h system for a range of
% h values. For the open-loop simulations it finds the minimum and maximum h
% values observed and for the closed-loop static h simulations it finds the minimum
% and maximum gtonic values observed. It also plots the boundaries of
% quiscence, bursting, and beating regimes. These boundaries can be found
% by inspecting the spike time statistics of the simulations.

clear all; close all

tf=12e4;

options=odeset('RelTol',1e-8,'AbsTol',1e-8);

%% open-loop

global gtonic_open

gtonic_vals=0.05:0.01:0.7;

% Initial conditions [V n h alpha vollung PO2lung PO2blood]
inits0 = [-60 0 0 0 2 110 110];

for gx=1:length(gtonic_vals)
    
    gtonic_open = gtonic_vals(gx)

    [t0,u0]=ode15s('openloop',[0 tf],inits0,options);

    inits1=u0(end,:);
    [t1,u1]=ode15s('openloop',[0 tf],inits1,options);
    
    t=t1;
    v=u1(:,1);
    h=u1(:,3);
    
    % find spike time statistics
      
    spikeTimes=[];
    thresh=0;

    for jx=2:length(t)
        if v(jx)>=thresh
            if v(jx-1)<thresh
                spikeTimes=[spikeTimes t(jx)];
            end
        end
    end
    
    numSpikes_open(gx)=length(spikeTimes);
    
    if length(spikeTimes)>1
    
        isi=diff(spikeTimes);
        maxisi_open(gx)=max(isi);
        cv_open(gx)=std(isi)/mean(isi);
    
    else
        
        maxisi_open(gx)=NaN;
        cv_open(gx)=NaN;
        
    end
    
    % find range of h values traversed
    
    hrange_open(gx,:)=[min(h) max(h)];
    
end

%% closed loop static h

hvals=0.2:0.01:0.9;

for hx=1:length(hvals)
    
    h0=hvals(hx)
    
    % Initial conditions [V n h alpha vollung PO2lung PO2blood]
    inits0 = [-60 0 h0 0 2 110 110];
    
    [t0,u0]=ode15s('closedloop_hstatic',[0 tf],inits0,options);

    inits1=u0(end,:);
    [t1,u1]=ode15s('closedloop_hstatic',[0 tf],inits1,options);

    inits2=u1(end,:);
    [t2,u2]=ode15s('closedloop_hstatic',[0 15000],inits2,options);
    
    t=t2;
    v=u2(:,1);
    PO2blood=u2(:,7);
    gtonic=0.3*(1-tanh((PO2blood-85)./30));
    
    % find spike time statistics
      
    spikeTimes=[];
    thresh=0;

    for jx=2:length(t)
        if v(jx)>=thresh
            if v(jx-1)<thresh
                spikeTimes=[spikeTimes t(jx)];
            end
        end
    end
    
    numSpikes_hstatic(hx)=length(spikeTimes);
    
    if length(spikeTimes)>1
    
        isi=diff(spikeTimes);
        maxisi_hstatic(hx)=max(isi);
        cv_hstatic(hx)=std(isi)/mean(isi);
        
    else
        
        maxisi_hstatic(hx)=NaN;
        cv_hstatic(hx)=NaN;
        
    end
    
    % find range of gtonic values traversed
    
    grange_hstatic(hx,:)=[min(gtonic) max(gtonic)];
   
end

%% Initial conditions [V n h alpha vollung PO2lung PO2blood]

initsA = [-58.5754 0.0006 0.7252 0.0010 2.2665 103.3461 102.2229];
initsB = [-41.7429 0.0313 0.3442 0.0025 2.4355 23.9533 23.3940];

tf = 15000;

options = odeset('RelTol',1e-9,'AbsTol',1e-9);

[tA,uA] = ode15s('closedloop',[0 tf],initsA,options);
[tB,uB] = ode15s('closedloop',[0 tf],initsB,options);

hA=uA(:,3);
PO2bloodA=uA(:,7);
gtonicA=0.3*(1-tanh((PO2bloodA-85)./30));

hB=uB(:,3);
PO2bloodB=uB(:,7);
gtonicB=0.3*(1-tanh((PO2bloodB-85)./30));


%% make plot

figure(1)

hold on
% boundaries of Q/Be/Bu/Be for closed-loop static h
plot([0 1],[.35 .35],'r--','Linewidth',1)
plot([0 1],[.46 .46],'r--','Linewidth',1)
plot([0 1],[.78 .78],'r--','Linewidth',1)

% boundaries of Q/Bu/Be for open-loop
plot([.27 .27],[0 1],'b--','Linewidth',1)
plot([.46 .46],[0 1],'b--','Linewidth',1)

% min and max h values during open-loop
plot(gtonic_vals,hrange_open(:,1),'b','Linewidth',3)
plot(gtonic_vals,hrange_open(:,2),'b','Linewidth',3)

% min and max gtonic values during closed-loop static h
plot(grange_hstatic(:,1),hvals,'Color',[1 0 0],'Linewidth',3)
plot(grange_hstatic(:,2),hvals,'Color',[1 0 0],'Linewidth',3)

% closed-loop (dynamic h) eupneic limit cycle
plot(gtonicA,hA,'k','Linewidth',5)

% closed-loop (dynamic h) tachypneic fixed point
plot(gtonicB,hB,'ko','MarkerSize',10,'MarkerFaceColor','k')

gtonicA_max_ind=find(gtonicA==max(gtonicA));
gtonicA_min_ind=find(gtonicA==min(gtonicA));
gtonicA_0pt18_ind=find(abs(gtonicA-0.18)==min(abs(gtonicA-0.18)));

plot(gtonicA(gtonicA_max_ind),hA(gtonicA_max_ind),'go','MarkerSize',10,'MarkerFaceColor','g')
plot(gtonicA(gtonicA_min_ind),hA(gtonicA_min_ind),'co','MarkerSize',10,'MarkerFaceColor','c')
plot(gtonicA(gtonicA_0pt18_ind),hA(gtonicA_0pt18_ind),'mo','MarkerSize',10,'MarkerFaceColor','m')

ylim([0.2 0.9])
xlim([0.05 0.7])

ylabel('$h$','Interpreter','Latex')
xlabel('$g_\mathrm{tonic}$','Interpreter','Latex')
set(gca,'YTick',[0.1:.1:1])




    
    

    

Loading data, please wait...