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 5
% first the open-loop, then closed-loop, then gtonic forcing

clear all

scale_factors=[0.8 1 1.5 2 2.5 3 3.5 4 4.5];

%% initialize open and closed-loop

global taumax_h gtonic_open

v0=-55; n0=0.001; h0=0.74; alpha0=0.0001; vollung0=2; PO2lung0=100; PO2blood0=100; 

inits1_open_closed=[v0 n0 h0 alpha0 vollung0 PO2lung0 PO2blood0];

t0=0;
tf1=60*1000;
tf2=60*1000*2;
tf3=60*1000*3;
tf4=60*1000*4;
tf5=60*1000*5;
tf6=60*1000*6;

options=odeset('RelTol',1e-10,'AbsTol',1e-10,'Events',@events_figure5);

gtonic_open=0.3;


%% initialize forced gtonic

inits1_forced=[-56.031664239999998 0.001160025770000 0.753151682700000 0.000095924685630 2.024680252000000 97.636469689999998 96.467834690000004];

tf=(6.372270214e3)*20;
tin=0:0.01:tf;

[t1,u1,te1,ye1,ie1]=ode15s(@closedloop,[0 tf],inits1_forced,options);

data=[t1 u1];


%% loop through scale factors

for scalex=1:length(scale_factors)
    
    scale_factor=scale_factors(scalex)
       
    taumax_h=scale_factor*10000;

    % open-loop

    [t1_open,u1_open,te1_open,ye1_open,ie1_open]=ode15s(@openloop_tauh,[t0 tf1],inits1_open_closed,options);

    inits2_open=u1_open(end,:);
    [t2_open,u2_open,te2_open,ye2_open,ie2_open]=ode15s(@openloop_tauh,[t0 tf2],inits2_open,options);

    inits3_open=u2_open(end,:);
    [t3_open,u3_open,te3_open,ye3_open,ie3_open]=ode15s(@openloop_tauh,[t0 tf3],inits3_open,options);

    inits4_open=u3_open(end,:);
    [t4_open,u4_open,te4_open,ye4_open,ie4_open]=ode15s(@openloop_tauh,[t0 tf4],inits4_open,options);

    inits5_open=u4_open(end,:);
    [t5_open,u5_open,te5_open,ye5_open,ie5_open]=ode15s(@openloop_tauh,[t0 tf5],inits5_open,options);

    inits6_open=u5_open(end,:);
    [t6_open,u6_open,te6_open,ye6_open,ie6_open]=ode15s(@openloop_tauh,[t0 tf6],inits6_open,options);

    % closed-loop

    [t1_closed,u1_closed,te1_closed,ye1_closed,ie1_closed]=ode15s(@closedloop_tauh,[t0 tf1],inits1_open_closed,options);

    inits2_closed=u1_closed(end,:);
    [t2_closed,u2_closed,te2_closed,ye2_closed,ie2_closed]=ode15s(@closedloop_tauh,[t0 tf2],inits2_closed,options);

    inits3_closed=u2_closed(end,:);
    [t3_closed,u3_closed,te3_closed,ye3_closed,ie3_closed]=ode15s(@closedloop_tauh,[t0 tf3],inits3_closed,options);

    inits4_closed=u3_closed(end,:);
    [t4_closed,u4_closed,te4_closed,ye4_closed,ie4_closed]=ode15s(@closedloop_tauh,[t0 tf4],inits4_closed,options);

    inits5_closed=u4_closed(end,:);
    [t5_closed,u5_closed,te5_closed,ye5_closed,ie5_closed]=ode15s(@closedloop_tauh,[t0 tf5],inits5_closed,options);

    inits6_closed=u5_closed(end,:);
    [t6_closed,u6_closed,te6_closed,ye6_closed,ie6_closed]=ode15s(@closedloop_tauh,[t0 tf6],inits6_closed,options);

    % forced gtonic
    
    taumax_h=10000;

    tdata=scale_factor*data(:,1);
    vdata=data(:,2);
    ndata=data(:,3);
    hdata=data(:,4);
    po2blood=data(:,end);
    gtonic_orig=0.3*(1-tanh((po2blood-85)./30));

    % parameters
    gnap=2.8; gna=28; gk=11.2; gl=2.8; 
    Ena=50; Ek=-85; El=-65; Esyn=0;

    C=21;  %pF

    theta_mp=-40;  %mV
    sigma_mp=-6;  %mV
    theta_h=-48;  %mV
    sigma_h=6;  %mV

    theta_m=-34;  %mV
    sigma_m=-5;  %mV

    taumax_n=10;  %ms
    theta_n=-29; %mV
    sigma_n=-4;  %mV

    mp_inf=@(v) 1/(1+exp((v-theta_mp)/sigma_mp));
    m_inf=@(v) 1/(1+exp((v-theta_m)/sigma_m));
    h_inf=@(v) 1/(1+exp((v-theta_h)/sigma_h));
    tau_h=@(v) taumax_h/cosh((v-theta_h)/(2*sigma_h));
    n_inf=@(v) 1/(1+exp((v-theta_n)/sigma_n));
    tau_n=@(v) taumax_n/cosh((v-theta_n)/(2*sigma_n));

    tf=max(tdata);
    dt=0.01;

    t=0:dt:tf;

    gtonic=interp1(tdata,gtonic_orig,t);

    numSteps=length(t)-1;

    v=zeros(length(t),1);
    n=zeros(length(t),1);
    h=zeros(length(t),1);

    v(1)=vdata(1);
    n(1)=ndata(1);
    h(1)=hdata(1);

    for ix=1:numSteps

        Inap=gnap*mp_inf(v(ix))*h(ix)*(v(ix)-Ena);
        Ina=gna*(m_inf(v(ix))^3)*(1-n(ix))*(v(ix)-Ena);
        Ik=gk*(n(ix)^4)*(v(ix)-Ek);
        Il=gl*(v(ix)-El);
        Itonic=gtonic(ix)*(v(ix)-Esyn);

        k1v=(-Inap-Ina-Ik-Il-Itonic)/C;
        k1n=(n_inf(v(ix))-n(ix))/tau_n(v(ix));
        k1h=(h_inf(v(ix))-h(ix))/tau_h(v(ix));

        av=v(ix)+k1v*dt;
        an=n(ix)+k1n*dt;
        ah=h(ix)+k1h*dt;

        Inap=gnap*mp_inf(av)*ah*(av-Ena);
        Ina=gna*(m_inf(av)^3)*(1-an)*(av-Ena);
        Ik=gk*(an^4)*(av-Ek);
        Il=gl*(av-El);
        Itonic=gtonic(ix+1)*(av-Esyn);

        k2v=(-Inap-Ina-Ik-Il-Itonic)/C;
        k2n=(n_inf(av)-an)/tau_n(av);
        k2h=(h_inf(av)-ah)/tau_h(av);

        v(ix+1)=v(ix)+(k1v+k2v)*dt/2;
        n(ix+1)=n(ix)+(k1n+k2n)*dt/2;
        h(ix+1)=h(ix)+(k1h+k2h)*dt/2;

    end
    
    t_forced=t;
    v_forced=v;


    %%% find spike times and burst properties

    %% open-loop

    t=t6_open;
    v=u6_open(:,1);

    numT=length(t);
    spikeTimes=[];
    tlap=0;

    for ix=2:numT
        if v(ix)>0
            if v(ix-1)<0
                if (t(ix)-tlap)>10
                    spikeTimes=[spikeTimes t(ix)];
                    tlap=t(ix);
                end
            end
        end
    end

    numSpikes=length(spikeTimes);
    isi=spikeTimes(2:end)-spikeTimes(1:(end-1));

    max_isi=max(isi);

    ibi_inds=find(isi>0.8*max_isi);
    ibi=isi(ibi_inds);
    num_ibi=length(ibi);
    mean_ibi=mean(ibi/1000);

    first_spike_inds=ibi_inds+1;
    first_spikeTimes=spikeTimes(first_spike_inds);

    last_spike_inds=ibi_inds;
    last_spikeTimes=spikeTimes(last_spike_inds);

    period_first=first_spikeTimes(2:end)-first_spikeTimes(1:(end-1));
    period_last=last_spikeTimes(2:end)-last_spikeTimes(1:(end-1));

    mean_period=mean([period_first period_last])/1000;
    mean_burst_dur=mean_period-mean_ibi;

    spikes_per_burst=length(first_spike_inds(end-1):last_spike_inds(end));

    openIBI(scalex)=mean_ibi;
    openPeriod(scalex)=mean_period;
    openBurstDur(scalex)=mean_burst_dur;
    openSPB(scalex)=spikes_per_burst;

    %% closed-loop

    t=t6_closed;
    v=u6_closed(:,1);

    numT=length(t);
    spikeTimes=[];
    tlap=0;

    for ix=2:numT
        if v(ix)>0
            if v(ix-1)<0
                if (t(ix)-tlap)>10
                    spikeTimes=[spikeTimes t(ix)];
                    tlap=t(ix);
                end
            end
        end
    end

    numSpikes=length(spikeTimes);
    isi=spikeTimes(2:end)-spikeTimes(1:(end-1));

    max_isi=max(isi);

    ibi_inds=find(isi>0.8*max_isi);
    ibi=isi(ibi_inds);
    num_ibi=length(ibi);
    mean_ibi=mean(ibi/1000);

    first_spike_inds=ibi_inds+1;
    first_spikeTimes=spikeTimes(first_spike_inds);

    last_spike_inds=ibi_inds;
    last_spikeTimes=spikeTimes(last_spike_inds);

    period_first=first_spikeTimes(2:end)-first_spikeTimes(1:(end-1));
    period_last=last_spikeTimes(2:end)-last_spikeTimes(1:(end-1));

    mean_period=mean([period_first period_last])/1000;
    mean_burst_dur=mean_period-mean_ibi;

    spikes_per_burst=length(first_spike_inds(end-1):last_spike_inds(end));

    closedIBI(scalex)=mean_ibi;
    closedPeriod(scalex)=mean_period;
    closedBurstDur(scalex)=mean_burst_dur;
    closedSPB(scalex)=spikes_per_burst;


    %% gtonic forcing
    
    t=t_forced;
    v=v_forced;

    numT=length(t);
    spikeTimes=[];
    tlap=0;

    for ix=2:numT
        if v(ix)>0
            if v(ix-1)<0
                if (t(ix)-tlap)>10
                    spikeTimes=[spikeTimes t(ix)];
                    tlap=t(ix);
                end
            end
        end
    end

    numSpikes=length(spikeTimes);
    isi=spikeTimes(2:end)-spikeTimes(1:(end-1));

    max_isi=max(isi);

    ibi_inds=find(isi>0.8*max_isi);
    ibi=isi(ibi_inds);
    num_ibi=length(ibi);
    mean_ibi=mean(ibi/1000);

    first_spike_inds=ibi_inds+1;
    first_spikeTimes=spikeTimes(first_spike_inds);

    last_spike_inds=ibi_inds;
    last_spikeTimes=spikeTimes(last_spike_inds);

    period_first=first_spikeTimes(2:end)-first_spikeTimes(1:(end-1));
    period_last=last_spikeTimes(2:end)-last_spikeTimes(1:(end-1));

    mean_period=mean([period_first period_last])/1000;
    mean_burst_dur=mean_period-mean_ibi;

    spikes_per_burst=length(first_spike_inds(end-1):last_spike_inds(end));

    forcedIBI(scalex)=mean_ibi;
    forcedPeriod(scalex)=mean_period;
    forcedBurstDur(scalex)=mean_burst_dur;
    forcedSPB(scalex)=spikes_per_burst;

end

openRec=[openIBI' openPeriod' openBurstDur' openSPB'];
closedRec=[closedIBI' closedPeriod' closedBurstDur' closedSPB'];
forcedRec=[forcedIBI' forcedPeriod' forcedBurstDur' forcedSPB'];

dlmwrite('open_IBI_BurstDur_SPB.csv',openRec,'precision',10)
dlmwrite('closed_IBI_BurstDur_SPB.csv',closedRec,'precision',10)
dlmwrite('forced_IBI_BurstDur_SPB.csv',forcedRec,'precision',10)


%% make plot

figure(1)

lw=2;

subplot(3,1,1)
hold on
plot(scale_factors,openIBI,'bs-','Linewidth',lw,'MarkerFaceColor','b')
plot(scale_factors,closedIBI,'ko-','Linewidth',lw,'MarkerFaceColor','k')
plot(scale_factors,forcedIBI,'gs-','Linewidth',lw,'MarkerFaceColor','g')
ylabel('IBI (s)','Interpreter','Latex')
xlabel('Scale Factor $\gamma$','Interpreter','Latex')
xlim([0.8 4.5])
set(gca,'box','off','XTick',[1:4],'YTick',[0:5:30])
h=legend(' $\bar{\tau}_h = \gamma \times 10^4 \quad \mbox{(open loop)}$',' $\bar{\tau}_h = \gamma \times 10^4 \quad \mbox{(closed loop)}$',' $\tau_{P_\mathrm{a}\mathrm{O}_2} =  \gamma$','Location','Northwest');
legend('boxoff')
set(h,'Interpreter','latex','fontsize',24)
grid on

subplot(3,1,2)
hold on
plot(scale_factors,openBurstDur,'bs-','Linewidth',lw,'MarkerFaceColor','b')
plot(scale_factors,closedBurstDur,'ko-','Linewidth',lw,'MarkerFaceColor','k')
plot(scale_factors,forcedBurstDur,'gs-','Linewidth',lw,'MarkerFaceColor','g')
ylabel('Burst Duration (s)','Interpreter','Latex')
xlabel('Scale Factor $\gamma$','Interpreter','Latex')
xlim([0.8 4.5])
set(gca,'box','off','XTick',[1:4],'YTick',[0:.5:2.5])
grid on

subplot(3,1,3)
hold on
plot(scale_factors,openSPB,'bs-','Linewidth',lw,'MarkerFaceColor','b')
plot(scale_factors,closedSPB,'ko-','Linewidth',lw,'MarkerFaceColor','k')
plot(scale_factors,forcedSPB,'gs-','Linewidth',lw,'MarkerFaceColor','g')
ylabel('Spikes Per Burst','Interpreter','Latex')
xlabel('Scale Factor $\gamma$','Interpreter','Latex')
xlim([0.8 4.5])
set(gca,'box','off','XTick',[1:4],'YTick',[0:10:60])
ylim([0 60])
grid on





Loading data, please wait...