% tristable_log.m
% for alpha12 = alpha13, alpha21 = alpha23, alpha31 = alpha32
i=5;
alpha1 = 10.^([1:i]);%transcription plus translation rate of pLac construct.  need to have same # of each alpha
alpha2 = 10.^([1:i]); %transcription plus translation rate of pTet construct
alpha3 = 10.^([1:i]); %transcription plus translation rate of pBAD construct
beta1 = 2; %Cooperativity ie when there are more repressors repression goes up by more than the increase in repressors
beta2 = 2;
beta3 = 2;
beta4 = 2;
beta5 = 2;
beta6 = 2;
J = zeros(3,3);%creates empty mat for Jacobian
digits(10); %# of significant digits used in vpa relavant calculations.  default is 32
count=0;
isstable = zeros(numel(alpha1), numel(alpha2), numel(alpha3));
%a 1 will replace each indece which is stable (indeces do not
%necessarily correspond to alpha values)
 
% plotting parameters
mksize = 10;%size of marker dots in plot
x = 0;
for a1 = 1:length(alpha1)
    for a2 = 1:length(alpha2)
        count%used to debug (see if program is still running/running the correct number of times)
        for a3 = 1:length(alpha3)
            count = count +1;
            eq = 0; %reset number of stable states
            Jisnan = 0;
            syms u v w; %creates the variables
            f = vpa(alpha1(a1)/(1+v^beta1)+alpha1(a1)/(1+w^beta2)-u);
            g = vpa(alpha2(a2)/(1+u^beta3)+alpha2(a2)/(1+w^beta4)-v);
            h = vpa(alpha3(a3)/(1+u^beta5)+alpha3(a3)/(1+v^beta6)-w);
            [x,y,z] = solve(f,g,h,u,v,w); %solves eqns f, g and h at the roots for u, v and w using vpa standard # of digits
            sizex = size(x);
            for ss = 1:sizex(1)%takes first dimension of size x (always 21x1)
                uu = double(x(ss)); %must use double type otherwise isreal() always returns 0
                vv = double(y(ss));
                ww = double(z(ss));
                if isreal(uu) && isreal(vv) && isreal(ww)
                    %compute jacobian
                    J(1,1) = -1;%du/dudt of alpha1/(v^beta1 +1)+alpha2/(z^beta2+1) -u
                    J(1,2) = -alpha1(a1)*beta1*vv^(beta1-1)/((1+vv^beta1)^2);
                    J(1,3) = -alpha1(a1)*beta2*ww^(beta2-1)/((1+ww^beta2)^2);
                    J(2,1) = -alpha2(a2)*beta3*uu^(beta3-1)/((1+uu^beta3)^2);
                    J(2,2) = -1;
                    J(2,3) = -alpha2(a2)*beta4*ww^(beta4-1)/((1+ww^beta4)^2);
                    J(3,1) = -alpha3(a3)*beta5*uu^(beta5-1)/((1+uu^beta5)^2);
                    J(3,2) = -alpha3(a3)*beta6*vv^(beta6-1)/((1+vv^beta6)^2);
                    J(3,3) = -1;
                    ev = eig(J);
                    if ev(1)<0 && ev(2)<0 && ev(3)<0 %if eigen value is negative then that set is stable, positive values are unstable
                        eq = eq + 1;
                    end
                end
            end
            if eq >= 3 %need three stable solutions for tristable switch
                isstable(a1,a2,a3) = 1;% puts 1 where switch is stable
                plot3(log10(alpha1(a1)),log10(alpha2(a2)),log10(alpha3(a3)),'b.','markersize',20)
                hold on
            %else
                %plot3(log10(alpha1(a1)),log10(alpha2(a2)),log10(alpha3(a3)),'ro','markersize',mksize-4)
                %hold on
            end
            %plot3(0:10, 0:10, 0:10, 'g-')%plots 1:1:1 lie for better readability in 3D
            xlim(log10(alpha1([1 end])))
            ylim(log10(alpha2([1 end])))
            zlim(log10(alpha3([1 end])))
            hold on
 
        end
    end
end
xlabel('alpha1 (log10)')%annotates plot
ylabel('alpha2 (log10)')
zlabel('alpha3 (log10)')
%title('Plots Region of Stable Alpha Values (production rate) in log10', 'fontsize', 14)
%legend('stable points');
csvwrite('U:\My Documents\iGEM\isstable102207.dat',isstable) %
%writes to csv file, matrix, row, col
count
hold off