function repressilator_paramater_dependency()
    % Repressilatorモデル
    % ODEを数値的に解く
    %% パラメータの指定
    s0 = [1, 0.5, 0.0, 1, 1, 1]; % 各mRNA、タンパク質の初期値
    param = [100, 1, 5, 2]; % パラメータ a, a0, n
    time = 0.01:0.1:100; %シミュレーションを行う時間
    a0 = -5:1:1; % パラメータを振る
    a0 = 10.^a0; % パラメータを振る
    col = jet; % 色指定に使用
    figure(1);
    clf;
    figure(2);
    clf;
    %% パラメータを振り分けて計算
    for i = 1:numel(a0);
        param(2) = a0(i); %% 振ったパラメータを代入
        
        %% 数値計算
        [t, time_course] = ode45(@(t, s) ODE(t, s, param), time, s0); % ODEを解く
        
        %% figure 1 を描く
        figure(1);
        hold on;
        subplot(numel(a0), 1, i);
        plot(　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　); % 横軸:時間、縦軸:濃度でプロット
        %         axis([0 time(end) -2 2]);
        
        %% figure 2 を描く
        figure(2);
        ph(i)=plot3(　　　　　　　　　　　　,　　　　　　　　　　　, 'color', col(floor(53/numel(a0)*i),:)); % 横軸:TetR、縦軸:LacI、高さ軸: λcIでプロット
        hold on;
        str_tmp = 'a0 = Input_tmp';
        leg1(i,:) = strrep(str_tmp, 'Input_tmp', num2str(a0(i), '%+6.1e'));
        
    end
    
    %% figure 2 の凡例
    figure(2);
    hleg = legend([ph],leg1);
    %     set(hleg,'Location','NorthEastOutside')
    xlabel('x');
    ylabel('y');
    zlabel('z');
end

function dsdt = ODE(t, s, param)
    x = s(1); % tetR-lite
    y = s(2); % lacI-lite
    z = s(3); % λcI-lite
    X = s(4); % TetR
    Y = s(5); % LacI
    Z = s(6); % λcI
    
    alpha = param(1);
    alpha0 = param(2);
    beta = param(3);
    n = param(4);
    
    dsdt(1, :) = 　　　　　　　　　　　　　　; % tetR-lite
    dsdt(2, :) = 　　　　　　　　　　　　　　; % lacI-lite
    dsdt(3, :) = 　　　　　　　　　　　　　　; % λcI-lite
    dsdt(4, :) = 　　　　　　　　　　　　　　; % TetR
    dsdt(5, :) = 　　　　　　　　　　　　　　; % LacI
    dsdt(6, :) = 　　　　　　　　　　　　　　; % λcI
end

