%infToy()
%
%Toy model for the simulation of the effect of some mutations on the
%inflammatory response of a cell.
%
%Simulation of a cell responding to a inflammatory signal IL-1B.
%It provides the solution of the original system as "Wild type" and also
%simuates the effect on the production of IL-6 and IL-10 after mutation on
%IRAK4 and IKK.
%
%ODE toy model for inflammation example, compare with Petri Nets version.
%Material for the review for Frontiers in Genetics.
%
%Guido Santos, 14th october 2015. Ver: 1.0
%

function infToy
tic
%Parameter values
par = [0.095494006372294021     %01
       0.095494006372294021     %02
       0.00095092672985701396   %03
       3.148117756260155        %04
       6.6237841144233673       %05
       0.14969105208857217      %06
       0.14969105208857217      %07
       83.992067558019102       %08 *10
       0.78238239329941495      %09
       0.011554688772929096     %10
       1.7167747622178087       %11
       540.25955106495292       %12
       1.5603360734858587       %13
       19.809822965869529       %14
       123.7572891589779        %15
       0.00070254758737242802   %16
       252.63310249152755       %17
       0.59996940751789452      %18
       0.074756548942597978     %19
       0.13584846690209032      %20
       1                        %21
       1                        %22
       1                        %23
       1                        %24
       0.14969105208857217      %25
       0.14969105208857217      %26
       83.992067558019102       %27
       0.78238239329941495      %28 
       0.14969105208857217      %29
       0.14969105208857217      %30
       83.992067558019102       %31
       0.78238239329941495      %32
       0.00070254758737242802   %33
       252.63310249152755       %34
       0.59996940751789452      %35
       0.074756548942597978     %36
       0.13584846690209032];    %37
par2 = par;

%Name of the variables
variables = {'IRAK1','IRAK1p','IKK','IKKa','NFkB','mIkBa','IkBa','mIL6','IL-6','MAPK','MAPKa','AP1','AP1a','mIL10','IL-10'};

[t1,y1] = ode45(@eqs_inf,[0 50],[1 0 1 0 0 0 1 0 0 1 0 1 0 0 0],[],par);

%Plot with the interleuklins
figure(1)
c = 0;
par2 = par;
par2(1) = par2(1)*0.1; %IRAK-
[t2,y2] = ode45(@eqs_inf,[0 50],[1 0 1 0 0 0 1 0 0 1 0 1 0 0 0],[],par2);

for n = [9 15]
    c = c + 1;
    subplot(2,2,c)
    plot(t1,y1(:,n))
    hold on
    plot(t2,y2(:,n),'r')
    hold off
    legend('Wild type','IRAK1^-','location','SouthEast')
    title(variables{n},'FontSize',14)
    xlabel('Time (hours)','FontSize',14)
    ylabel('Amount (arbitrary units)','FontSize',14)
    set(gca,'FontSize',12)
end

par2 = par;
par2(6) = par2(6)*0.1; %IKK-
[t2,y2] = ode45(@eqs_inf,[0 50],[1 0 1 0 0 0 1 0 0 1 0 1 0 0 0],[],par2);

for n = [9 15]
    c = c + 1;
    subplot(2,2,c)
    plot(t1,y1(:,n))
    hold on
    plot(t2,y2(:,n),'r')
    hold off
    legend('Wild type','IKK^-','location','SouthEast')
    title(variables{n},'FontSize',14)
    xlabel('Time (hours)','FontSize',14)
    ylabel('Amount (arbitrary units)','FontSize',14)
    set(gca,'FontSize',12)
end

c   = 0;
esp = [0.1 0.105 0.11 0.15 0.2 0.25 0.5 1 2 5 7 10];

for k = esp
    c       = c + 1;
    par2    = par;
    par2(1) = par2(1)*k; %IRAK-
	[~,y2]  = ode45(@eqs_inf,[0 50],[1 0 1 0 0 0 1 0 0 1 0 1 0 0 0],[],par2);
    k1(c)   = k;
    y91(c)  = max(y2(:,9));
    y151(c) = max(y2(:,15));
end

c = 0;

for k = esp
    c       = c + 1;
    par2    = par;
    par2(6) = par2(6)*k; %IKK-
    [~,y2]  = ode45(@eqs_inf,[0 50],[1 0 1 0 0 0 1 0 0 1 0 1 0 0 0],[],par2);
    k2(c)   = k;
    y92(c)  = max(y2(:,9));
    y152(c) = max(y2(:,15));
end

figure(2)
subplot(2,1,1),semilogx(k1,y91,k2,y92)
ylabel('Maximum value (arbitrary units)','FontSize',14)
xlabel('Value of mutation','FontSize',14)
xlim([0.1 10])
title('IL-6','FontSize',14)
legend('IRAK1^-','IKK^-')
set(gca,'FontSize',12)
subplot(2,1,2),semilogx(k1,y151,k2,y152)
ylabel('Maximum value (arbitrary units)','FontSize',14)
xlabel('Value of mutation','FontSize',14')
xlim([0.1 10])
title('IL-10','FontSize',14)
legend('IRAK1^-','IKK^-')
set(gca,'FontSize',12)

f       = 0;
espacio = linspace(-2,2,51);

for k1 = espacio
    f = f + 1;
    c = 0;
    for k2 = espacio
        c               = c + 1;
        par2            = par;
        par2([8 27 31]) = par2([8 27 31])*10^k1; %IRAK-
        par2(10)        = par2(10)*10^k2; %IKK-
        [~,y2]          = ode45(@eqs_inf,[0 50],[1 0 1 0 0 0 1 0 0 1 0 1 0 0 0],[],par2);
        y91(f,c)        = max(y2(:,9));
        y151(f,c)       = max(y2(:,15));
    end
end

y91plot  = fliplr(y91')';
y151plot = fliplr(y151')';

figure(3)
etiq1 = {'10^{-1}','10^0','10^1'};
etiq2 = {'10^{1}','10^0','10^{-1}'};
ptos  = [1 (length(espacio)+1)/2 length(espacio)];

subplot(2,1,1)
imagesc(y91plot)
set(gca,'XTick', ptos,'XTickLabel',etiq1,'YTick', ptos,'YTickLabel',etiq2)
ylabel('IRAK1^-')
xlabel('IKK^-')
title('IL-6')
colorbar

subplot(2,1,2)
imagesc(y151plot)
set(gca,'XTick', ptos,'XTickLabel',etiq1,'YTick', ptos,'YTickLabel',etiq2)
ylabel('IRAK1^-')
xlabel('IKK^-')
title('IL-10')
colorbar

toc

end

function dx = eqs_inf(~,x,par)

%Variable names
IRAK1  = x(1);
IRAK1p = x(2);
IKK    = x(3);
IKKa   = x(4);
NFkB   = x(5);
mIkBa  = x(6);
IkBa   = x(7);
mIL6   = x(8);
IL6    = x(9); 
MAPK   = x(10);
MAPKa  = x(11); 
AP1    = x(12); 
AP1a   = x(13); 
mIL10  = x(14);
IL10   = x(15);
       
%Parameter names
kirak1syn     = par(1);
kyrak1deg     = par(2);
kirak1ph1     = par(3);
kirak1ph2     = par(4);
kirak1pdeg    = par(5);
kikksyn       = par(6);
kikkdeg       = par(7);
kikkact       = par(8);
kikkadeg      = par(9);
knfkbgain     = par(10);
knfkbloss     = par(11);
kmikbatransc  = par(12);
kmikbadeg     = par(13);
kikbatransl   = par(14);
kikbaloss     = par(15);
kmil6transc1  = par(16);
kmil6transc2  = par(17);
kmil6deg      = par(18);
kil6transl    = par(19);
kil6deg       = par(20);
K1            = par(21);
K2            = par(22);
K3            = par(23);
Ntot          = par(24);
mapksyn       = par(25);
mapkdeg       = par(26);
mapkact       = par(27);
mapkadeg      = par(28);
ap1syn        = par(29);
ap1deg        = par(30);
ap1act        = par(31);
ap1adeg       = par(32);
kmil10transc1 = par(33);
kmil10transc2 = par(34);
kmil10deg     = par(35);
kil10transl   = par(36);
kil10deg      = par(37);
Flagellin     = 0;    
IL1B          = 1;    

%Differential equations
dx     = zeros(15,1);
dx(1)  = kirak1syn - kyrak1deg*IRAK1 - kirak1ph1*IRAK1*Flagellin - kirak1ph2*IRAK1*IL1B;
dx(2)  = kirak1ph1*IRAK1*Flagellin + kirak1ph2*IRAK1*IL1B - kirak1pdeg*IRAK1p;
dx(3)  = kikksyn - kikkact*IKK*IRAK1p - kikkdeg*IKK;
dx(4)  = kikkact*IKK*IRAK1p - kikkadeg*IKKa;
dx(5)  = knfkbgain*IKKa*((Ntot - NFkB)/(K1 + IkBa + IKKa)) - knfkbloss*((IkBa*NFkB)/(K2+NFkB));
dx(6)  = kmikbatransc*NFkB - kmikbadeg*mIkBa;
dx(7)  = kikbatransl*mIkBa - kikbaloss*IKKa*IkBa*((Ntot - NFkB)/(K3 + IkBa + IKKa));
dx(8)  = kmil6transc1 + kmil6transc2*NFkB - kmil6deg*mIL6;
dx(9)  = kil6transl*mIL6 - kil6deg*IL6;
dx(10) = mapksyn - mapkact*MAPK*IRAK1p - mapkdeg*MAPK;
dx(11) = mapkact*MAPK*IRAK1p - mapkadeg*MAPKa;
dx(12) = ap1syn - ap1act*AP1*IRAK1p - ap1deg*AP1;
dx(13) = ap1act*AP1*IRAK1p - ap1adeg*AP1a;
dx(14) = kmil10transc1 + kmil10transc2*AP1a - kmil10deg*mIL10;
dx(15) = kil10transl*mIL10 - kil10deg*IL10;

end
