function task3_stepwise(model_type)
% model_type
%1: Ras
%2: Rap1

time_basal = 1:10000; 
time = 1:1000;
fold_change = 5;
basal = 0.1;

figure;
hold on

if model_type ==1 % Ras model
    param = [0.5,5,0.0005,0.005,0.05,100]; % [k1,k2,k3,k4,k5,k6]
    
    % basal刺激後、定常に達するまでシミュレーション
    y0 = [basal,0,0,0]; % [S,GEFa,GAPa,Rasa]
    [t,time_course] = ode15s(@(t,y) ODE_step_Ras(t,y,param),time_basal,y0);

    % basal刺激後のシミュレーションの最後の値をt=0での初期値にする
    y0 =　;
    
    %t=0の刺激
    y0(1) = 　;
    
    % step刺激時の結果を作図
    [t,time_course] = ode15s(@(t,y) ODE_step_Ras(t,y,param),time,y0);
    plot(　,　, 'b'); %Rasaの応答の時間変化を図示
    
    % stepwise刺激時の結果を作図
    [t,time_course] = ode15s(@(t,y) ODE_stepwise_Ras(t,y,param, basal, fold_change),time,y0);
    plot(　,　, 'r'); %Rasaの応答の時間変化を図示
    
else % Rap1 model
    param = [0.5,0.5,0.2,5];
    
    % basal刺激後、定常に達するまでシミュレーション
    y0 = [basal,0,0.005,0]; % [S,GEFa,GAPa,Rap1a]
    [t,time_course] = ode15s(@(t,y) ODE_step_Rap1(t,y,param),time_basal,y0);

    % basal刺激後のシミュレーションの最後の値をt=0での初期値にする
    y0 =　;
    
    %t=0の刺激
    y0(1) = 　;
    
    % step刺激時の結果を作図
    [t,time_course] = ode15s(@(t,y) ODE_step_Rap1(t,y,param),time,y0);
    plot(　,　, 'b'); %Rap1aの応答の時間変化を図示
    
    % stepwise刺激時の結果を作図
    [t,time_course] = ode15s(@(t,y) ODE_stepwise_Rap1(t,y,param, basal, fold_change),time,y0);
    plot(　,　, 'r'); %Rap1aの応答の時間変化を図示
    
end

xlabel('time');
ylabel('response');
legend({'step', 'stepwise'})%凡例
if model_type == 1
    title('Ras model')
    ylim([0,0.004])%y軸の範囲を修正
else
    title('Rap1 model')
    ylim([0,1])%y軸の範囲を修正
end
end


%% Rasモデルでのstep刺激
function dydt = ODE_step_Ras(t,y,param)
S = y(1);
GEFa = y(2);
GAPa = y(3);
Rasa = y(4);

k1 = param(1);
k2 = param(2);
k3 = param(3);
k4 = param(4);
k5 = param(5);
k6 = param(6);

dydt(1,:) = 0;
dydt(2,:) = k1*S*(1-GEFa)-k2*GEFa;
dydt(3,:) = k3*S*(1-GAPa)-k4*GAPa;
dydt(4,:) = k5*GEFa*(1-Rasa)-k6*GAPa*Rasa;
end

%% Rasモデルでのstepwise刺激
function dydt = ODE_stepwise_Ras(t,y,param, basal, fold_change)
%Sを使ってステップワイズ刺激を表現
% 時間が500以降は刺激をfold change倍する
　　　　　
　　　　　
　　　　　
GEFa = y(2);
GAPa = y(3);
Rasa = y(4);

k1 = param(1);
k2 = param(2);
k3 = param(3);
k4 = param(4);
k5 = param(5);
k6 = param(6);

dydt(1,:) = 0;
dydt(2,:) = k1*S*(1-GEFa)-k2*GEFa;
dydt(3,:) = k3*S*(1-GAPa)-k4*GAPa;
dydt(4,:) = k5*GEFa*(1-Rasa)-k6*GAPa*Rasa;
end

%% Rap1モデルでのstep刺激
function dydt = ODE_step_Rap1(t,y,param)
S = y(1);
GEFa = y(2);
GAPa = y(3);
Rap1a = y(4);

k1 = param(1);
k2 = param(2);
k3 = param(3);
k4 = param(4);

dydt(1,:) = 0;
dydt(2,:) = k1*S*(1-GEFa)-k2*GEFa;
dydt(3,:) = 0;
dydt(4,:) = k3*GEFa*(1-Rap1a)-k4*GAPa*Rap1a;
end

%% Rap1モデルでのstepwise刺激
function dydt = ODE_stepwise_Rap1(t,y,param, basal, fold_change)
%Sを使ってステップワイズ刺激を表現
% 時間が500以降は刺激をfold change倍する
　　　　　
　　　　　
　　　　　
GEFa = y(2);
GAPa = y(3);
Rap1a = y(4);

k1 = param(1);
k2 = param(2);
k3 = param(3);
k4 = param(4);

dydt(1,:) = 0;
dydt(2,:) = k1*S*(1-GEFa)-k2*GEFa;
dydt(3,:) = 0;
dydt(4,:) = k3*GEFa*(1-Rap1a)-k4*GAPa*Rap1a;
end

