
function geneticnetwork

theta=1;

%plot a genetic network approximated by step and sigmoid functions
tspan = [0 2.5];
x0 = [0 0];

[t,x] = ode45(@sig_system, tspan, x0)
lineplot(t,x);
paraplot(t,x);

[t,x] = ode15s(@step_system, tspan, x0)
lineplot(t,x);
paraplot(t,x);

%the sigmoid function
function y=sigmoid(x,m)
    y=x^m/(theta+x^m);
end

%the step function
function y=step(x,sgn)
    if sgn*x<sgn*theta
        y=0;
    else
        y=1;
    end
end

%define a system of sigmoid equations
% g0'=g1
% g1'=~g0
function dx = sig_system(t, x)
dx = zeros(2, 1);
dx(1) = 2*sigmoid(x(2),11)-(1*x(1));
dx(2) = 2*(1-sigmoid(x(1),11))-(1*x(2));
end

%define a system of step equations
% g0'=g1
% g1'=~g0
function dx = step_system(t, x)
dx = zeros(2, 1);
dx(1) = 2*step(x(2),1)-(1*x(1));
dx(2) = 2*step(x(1),-1)-(1*x(2));
end

% find state changes
function [ys,states]=statechanges(t,x,theta)
    ys=[0];
    states=[];
    state=[0,0];
    ds=[0,0];
    [rows,cols]=size(x);
    for i = 2:rows
        for j=1:cols
            if x(i,j)>theta && x(i-1,j)<theta 
                ds(j)=1;
            end
            if x(i,j)<theta && x(i-1,j)>theta
                ds(j)=-1;
            end
        end
        if ds(1)~=0 || ds(2)~=0
            ys=cat(1,ys,[(t(i-1)+t(i))/2]);
            states=cat(1,states,state);
            state=state+ds;
            ds=[0,0];
        end
    end
    states=cat(1,states,state);
end

%plot concentrations of gene products by time
function lineplot(t,x)
    span=[min(min(x)),max(max(x))];
    [ys,states]=statechanges(t,x,theta);
    
    h=figure;
    hold on; axis([[0,max(t)],span(1)-.1,span(2)+.1])
    hd=plot(t,x);
    hc(1)=line([min(t),max(t)],[theta,theta],'color',[0,0,0]);
    for i=2:length(ys)
        hc(2*i)=line([ys(i),ys(i)],span,'color',[1,0,0]);
        hc(1+2*i) = text((ys(i-1)+ys(i))/2, span(2)/4, sprintf('[%d,%d]',states(i-1,:)), 'horizontalalignment', 'center');
    end
try
    hc(1+2*length(states)) = text((max(t)+ys(i))/2, span(2)/4, sprintf('[%d,%d]',states(length(states),:)), 'horizontalalignment', 'center');
catch
end
    xlabel('Time');
    ylabel('Concentration');
    legend('gene product 0','gene product 1','theta g0,g1');
end

%parametric plot of concentrations (state space)
function paraplot(t,x)
    span1=[min(min(x(:,1)),0),max(max(x(:,1)),2*theta)];
    span2=[min(min(x(:,2)),0),max(max(x(:,2)),2*theta)];
    h=figure;
    hold on; axis([span1,span2]);
    hd=plot(x(:,1),x(:,2));

    hc(1)=line([theta,theta],span2,'color',[0,0,0]);
    hc(2)=line(span1,[theta,theta],'color',[0,0,0]);
    hc(3) = text(theta*.5,theta*.5,'[0,0]', 'horizontalalignment', 'center');
    hc(3) = text(theta*.5,theta*1.5,'[0,1]', 'horizontalalignment', 'center');
    hc(3) = text(theta*1.5,theta*.5,'[1,0]', 'horizontalalignment', 'center');
    hc(3) = text(theta*1.5,theta*1.5,'[1,1]', 'horizontalalignment', 'center');

    xlabel('gene 0');
    ylabel('gene 1');
    legend('concentrations','thresholds');
end
    
end
