% This script reproduces Figure 7 by simulating the averaged system for % a range of M values and finding the fixed points (gbar=0) clear all %% Parameters global po2blood taulb R Temp R=62.364; Temp=310; taulb=500; po2bloodvals=110:-1:10; Mvals=2e-6:0.1e-6:18e-6; %% Run fast subsystem (closed-loop with constant PO2 blood) over a range of PO2 blood values for px=1:length(po2bloodvals) po2blood=po2bloodvals(px) po2bloodTimes10=round(10*po2blood); po2bloodTimes10plus1=round(10*(po2blood+.1)); % Initial Conditions v0=-55; n0=.001; h0=.75; alpha0=.0001; vollung0=2; po2lung0=po2blood; inits0=[v0 n0 h0 alpha0 vollung0 po2lung0]; t0=0; tf=60000*5; options=[]; [t0,u0]=ode15s('fastsubsystem',[t0 tf],inits0,options); inits=u0(end,:); %% Determine whether CPG is quiescient, spiking, or bursting options=odeset('RelTol',1e-8,'AbsTol',1e-8); newSol=ode15s(@fastsubsystem,[0 30000],inits,options); t=newSol.x; v=newSol.y(1,:); vollung=newSol.y(5,:); po2lung=newSol.y(6,:); Jlb=(1/taulb)*(po2lung-po2blood).*(vollung/(R*Temp)); numT=length(t); spikeInds=[]; tlap=0; for ix=2:numT if v(ix)>0 && v(ix-1)<0 if t(ix)>(tlap+1) spikeInds=[spikeInds ix]; tlap=t(ix); end end end spikeTimes=t(spikeInds); numSpikes=length(spikeTimes); if numSpikes<=1 numBursts=0; end if numSpikes>1 ISI=spikeTimes(2:end)-spikeTimes(1:(end-1)); largeISIinds=find(ISI>0.5*max(ISI)); numBursts=length(largeISIinds); if numBursts>1 largeISIs=ISI(largeISIinds); lastBurstStartInd=largeISIinds(end); secondToLastBurstStartInd=largeISIinds(end-1); lastBurstStartTime=spikeTimes(lastBurstStartInd+1); secondToLastBurstStartTime=spikeTimes(secondToLastBurstStartInd+1); lastBurstInd=find(t==lastBurstStartTime); secondToLastBurstInd=find(t==secondToLastBurstStartTime); end end %% Obtain one period of quiescent, spiking, or bursting solution if numBursts>1 period=t(lastBurstInd)-t(secondToLastBurstInd); tSub=t(secondToLastBurstInd:lastBurstInd); JlbSub=Jlb(secondToLastBurstInd:lastBurstInd); end if numBursts<=1 if numSpikes>1 period=t(spikeInds(end))-t(spikeInds(end-1)); tSub=t(spikeInds(end-1):spikeInds(end)); JlbSub=Jlb(spikeInds(end-1):spikeInds(end)); end if numSpikes<=1 period=1; tSub=t(end); JlbSub=Jlb(end); end end %% Calculate gbar for a range of M values for mx=1:length(Mvals) M=Mvals(mx); Hb=150; volblood=5; eta=Hb*1.36; gamma=volblood/22400; betaO2=0.03; c=2.5; K=26; SaO2=(po2blood^c)/(po2blood^c+K^c); partial=(c*po2blood^(c-1))*(1/(po2blood^c+K^c)-(po2blood^c)/((po2blood^c+K^c)^2)); CaO2=eta*SaO2+betaO2*po2blood; Jbt=M*CaO2*gamma; dxdt=(JlbSub-Jbt)/(gamma*(betaO2+eta*partial)); if numSpikes>1 intFy=trapz(tSub,dxdt); gbar(px,mx)=intFy/period; end if numSpikes<=1 gbar(px,mx)=dxdt; end end end %% Find location of fixed points (gbar=0) fprec=zeros(length(Mvals),3); for mx=1:length(Mvals) fp=[]; for px=2:length(po2bloodvals) if gbar(px-1,mx)>0 if gbar(px,mx)<0 fp=[fp mean(po2bloodvals((px-1):px))]; end end if gbar(px-1,mx)<0 if gbar(px,mx)>0 fp=[fp mean(po2bloodvals((px-1):px))]; end end end fprec(mx,1:length(fp))=fp'; end %% Make plots set(0,'DefaultAxesFontSize',24) % figure 8A figure(1) hold on plot(po2bloodvals,gbar(:,21),'g','Linewidth',3) plot(po2bloodvals,gbar(:,61),'c','Linewidth',3) plot(po2bloodvals,gbar(:,141),'m','Linewidth',3) plot(po2bloodvals,0*po2bloodvals,'k--','LInewidth',3) plot(fprec(21,1),0,'ko','MarkerSize',10,'MarkerFaceColor','g') plot(fprec(21,2),0,'ko','MarkerSize',10,'MarkerFaceColor','g') plot(fprec(21,3),0,'ko','MarkerSize',10,'MarkerFaceColor','g') plot(fprec(61,1),0,'ko','MarkerSize',10,'MarkerFaceColor','c') plot(fprec(61,2),0,'ko','MarkerSize',10,'MarkerFaceColor','c') plot(fprec(61,3),0,'ko','MarkerSize',10,'MarkerFaceColor','c') plot(fprec(141,1),0,'ko','MarkerSize',10,'MarkerFaceColor','m') h=legend('$M = 0.4 \times 10^{-5}$','$M = 0.8 \times 10^{-5}$','$M = 1.6 \times 10^{-5}$','Location','Southwest'); legend('boxoff') set(h,'Interpreter','latex','FontSize',20) xlabel('$P_\mathrm{a}\mathrm{O}_2$','Interpreter','latex','Fontsize',24) ylabel('$\bar{g}$','interpreter','latex','FontSize',24) ylim([-12.5e-3 3.5e-3]) xlim([10 100]) set(gca,'box','off','YTick',-12e-3:2e-3:4e-3) grid on % figure 8B figure(2) hold on plot(Mvals,fprec(:,1),'o','Color',[1 .5 0],'MarkerFaceColor',[1 .5 0],'MarkerSize',6) plot(Mvals,fprec(:,2),'o','Color',[1 .5 0]','MarkerFaceColor',[1 .5 0],'MarkerSize',6) plot(Mvals,fprec(:,3),'o','Color',[1 .5 0],'MarkerFaceColor',[1 .5 0],'MarkerSize',6) xlim([.2e-5 18e-6]) ylim([1 100]) h=legend('$\bar{g} = 0$','Location','Northeast'); set(h,'interpreter','latex','FontSize',24) legend('boxoff') xlabel('$M$','Interpreter','latex','FontSize',24) ylabel('$P_\mathrm{a}\mathrm{O}_2$','Interpreter','latex','Fontsize',24) set(gca,'XTick',[0.4e-5:.4e-5:1.6e-5]) grid on