function selkov_parameter_dependency()
    % selcovモデル
    % ODEを数値的に解く
    %% パラメータの指定
    s0 = [1, 1]; % ADP(x)の初期値とF6P(y)の初期値
    param = [0.06, 0.6]; % パラメータ a, b
    time = 0.01:0.1:400; %シミュレーションを行う時間
    b = 0.10:0.2:1.1; % パラメータを振る
    col = jet; % 色指定に使用
    
    %% パラメータを振り分けて計算
    for i = 1:numel(b);
        param(2) = b(i); %% 振ったパラメータを代入
        
        %% 数値計算
        [t, time_course] = ode45(@(t, s) ODE(t, s, param), time, s0); % ODEを解く
        
        %% figure 1 を描く
        figure(1);
        hold on;
        subplot(numel(b), 2, 2*i-1);
        plot(　　　　　　　　　　　　　　　　　　　); % 横軸:時間、縦軸:濃度でプロット
        axis([350 time(end) 0 3]);
        
        %% figure 2 を描く
        %figure(2);
        hold on;
        subplot(numel(b),2,2:2:12);
        %axis([0.0 2.5 0 3]);
        %plot_nullcline(param); % ヌルクラインも描きたいときは、先頭の%をはずす
        ph(i)=plot(　　　　　　,　　　　　　　, 'color', col(floor(56/numel(b)*i),:)); % 横軸:x、縦軸:yでプロット
        str_tmp = 'b = b_tmp';
        leg1(i,:) = strrep(str_tmp, 'b_tmp', num2str(b(i)));
        waitforbuttonpress
        
    end
    
    %% figure 2 の凡例
  %  figure(2);
    legend([ph],leg1);
    xlabel('ADP');
    ylabel('F6P');
end

function dsdt = ODE(t, s, param)
    x = s(1); % ADP
    y = s(2); % F6P
    a = param(1);
    b = param(2);
    
    v1 = 　　　　　　　　　;
    v2 = 　　　　　　　　　;
    v3 = 　　　　　　　　　;
    
    dsdt(1, :) = 　　　　　　　　　; % ADPの時間変化
    dsdt(2, :) = 　　　　　　　　　; % F6Pの時間変化
end

function plot_nullcline(param)
    a = param(1);
    b = param(2);
    x = 0:0.05:2.5;
    y1 = 　　　　　　　;
    y2 = 　　　　　　　;
    plot(　　　　　　　　);
end