function [dConcentration]= pathwayIGEM(~,y) 
%% ODE set for the simulation
%
%% set global variable
global par;

%% rename the parameter vector to match the names in the equation
K_l=par(1);                 % concentration of lipases
K_e=par(2);                 % concentration of esterases
K_k=par(3);                 % concentration of keratinases

K_m_lip = 1.03*10^(-3);      % Michaelis-Menton constant(Molar)
k_lip = 8.8*60*60;            %k_cat(per hour)

K_m_ester = 1.03*10^(-3);    % Michaelis-Menton constant (Molar)
k_ester = 8.8*60*60;          %k_cat(per hour)

K_m_kera = 1.03*10^(-3);     % Michaelis-Menton constant(Molar)
k_kera = 8.8*60*60;           %k_cat(per hour)

%**************************%
%**Initial Concentrations**%
%**************************%
% Initial_Concentration = Initial_Concentration_vector;
Initial_Concentration = y ;
%****************************************************%
%********************ODE FUNCTION********************%
%****************************************************%

dConcentration(1) = -(K_l*k_lip*Initial_Concentration(1)/(Initial_Concentration(1)+K_m_lip));
dConcentration(2) = (K_l*k_lip*Initial_Concentration(1)/(Initial_Concentration(1)+K_m_lip))-(K_e*k_ester*Initial_Concentration(2)/(Initial_Concentration(2)+K_m_ester));
dConcentration(3) = +(K_e*k_ester*Initial_Concentration(2)/(Initial_Concentration(2)+K_m_ester))-( K_k*k_kera*Initial_Concentration(3)/(Initial_Concentration(4)+K_m_kera));
dConcentration(4) = +(K_k*k_kera*Initial_Concentration(3)/(Initial_Concentration(4)+K_m_kera));
dConcentration(5) = 100000-dConcentration(4);

dConcentration = dConcentration';

end


