
% This code computes the fixed points of the system, sorts and filters
% them, computes their eigenvalues and stores them in realmccoy.

clear all
close all


At = 10E3; % Total concentration of kinase A
Bt = 10E3; % Total concentration of kinase B

maxtol = 2E-3; 

options=optimset('TolFun', maxtol, 'MaxFunEvals', 1000,'Display','off');   % Option to display output. Was 200


% defining the range of appropriate k values

k1_vector = logspace(-4, -2, 5);  
k2_vector = logspace(-3, 2, 11); 

% defining the range of appropriate initial conditions

IC_A_vector = logspace(0, 3.99, 10);
IC_B_vector = logspace(0, 3.99, 10);


master = 0;
beginning = 1;
eigcount = 1;
mess = 1;
bogus = 0;

duplicated = zeros(length(k1_vector), length(k2_vector));

for a = 1:length(k2_vector)     % vary k2
    
    k2 = k2_vector(a)
    
    for c = 1:length(k1_vector) % vary k1 
        
        k1 = k1_vector(c)
        
        k3 = k1/2;

        % testing if this combo of k1 and k2 is appropriate. If so, then continue

        
        if (k2/k1 >= 1E-3 & k2/k1 <= 1E6)
            count1 = 0;
            count2 = 0;
            count4 = 1;
            errortotal = 0;
            errortemp = 0;

            clear solutions1;
            clear solutions2;
            clear solutions3;
            clear solutions4;
            clear u;
            clear v;
            clear eigenval;
            clear k;
            clear D;

            % searching fsolve space for solutions
            for i = 1:length(IC_A_vector)         
                i
                for j = 1:length(IC_A_vector)
                    for k = 1:length(IC_B_vector)
                        for m = 1:length(IC_B_vector)
                            
                            x0 = [IC_A_vector(i);
                                  IC_A_vector(j);
                                  At - IC_A_vector(i) - IC_A_vector(j);
                                  IC_B_vector(k);
                                  IC_B_vector(m); 
                                  Bt - IC_B_vector(k) - IC_B_vector(m)];
                            

                            if ( x0(3) >= 0 & x0(6) >= 0 )

                                [x,fval] = fsolve(@goo_k3_v3, x0, options, k1, k2, k3, At, Bt);  % Call optimizer
                                
                                if (sum(goo_k3_v3(x, k1, k2, k3, At, Bt)) > 1)
                                    bogus = bogus + 1;
                                end
                                
                                count1 = count1 + 1;
                                
                                total_input(count1,:) = [x0(1) + x0(2) + x0(3), x0(4) + x0(5) + x0(6)];

                                solutions1(count1,1)  = x(1);               % A
                                solutions1(count1,2)  = x(2);               % Ai
                                solutions1(count1,3)  = x(3);   % Aii
                                solutions1(count1,4)  = x(4);               % B
                                solutions1(count1,5)  = x(5);               % Bi
                                solutions1(count1,6)  = x(6);   % Bii
                                solutions1(count1,7)  = fval(1);
                                solutions1(count1,8)  = fval(2);
                                solutions1(count1,9)  = fval(3);
                                solutions1(count1,10) = fval(4);
                                solutions1(count1,11) = fval(5);
                                solutions1(count1,12) = fval(6);
                                
                                if (abs(At - x(1) - x(2) - x(3)) > At * 0.001  | abs(Bt - x(4) - x(5) - x(6)) > Bt * 0.001)
                                    messed_up(mess,1) = x(1) + x(2) + x(3);
                                    messed_up(mess,2) = x(4) + x(5) + x(6);
                                    messed_up(mess, 3:4) = total_input(count1,:);
                                    mess = mess + 1;
                                end

                            end
                        end
                    end
                end
            end

            % testing for non-negative equilibrium points
            countr = 0;
            sol2 = 0;
            for i = 1:count1
                if (min(solutions1(i,1:6)) >= 0 & max(abs(solutions1(i,7:12))) <= maxtol)
                    count2 = count2 + 1;
                    solutions2(count2,:) = solutions1(i,:);
                    sol2 = 1;
                else
                    countr = countr + 1;
                    rejected(countr,:) = solutions1(i,:);
                end
            end

            if (sol2 == 1)
                % sorting all non-negative solutions
                solutions3 = sortrows(solutions2,1);


                % elimnating all redundant solutions, and testing for redundant answers
                % to ensure that I have thoroughly searched variable space?

                solutions4(1,:) = solutions3(1,:);

                for i = 2:count2
                    if (solutions3(i,1) + 1 > 1.4 * (solutions3(i-1,1) + 1) | solutions3(i,1) + 1 < 0.60 * (solutions3(i-1,1) + 1) | ...
                        solutions3(i,2) + 1 > 1.4 * (solutions3(i-1,2) + 1) | solutions3(i,2) + 1 < 0.60 * (solutions3(i-1,2) + 1) | ...
                        solutions3(i,3) + 1 > 1.4 * (solutions3(i-1,3) + 1) | solutions3(i,3) + 1 < 0.60 * (solutions3(i-1,3) + 1) | ...
                        solutions3(i,4) + 1 > 1.4 * (solutions3(i-1,4) + 1) | solutions3(i,4) + 1 < 0.60 * (solutions3(i-1,4) + 1) | ...
                        solutions3(i,5) + 1 > 1.4 * (solutions3(i-1,5) + 1) | solutions3(i,5) + 1 < 0.60 * (solutions3(i-1,5) + 1) | ...
                        solutions3(i,6) + 1 > 1.4 * (solutions3(i-1,6) + 1) | solutions3(i,6) + 1 < 0.60 * (solutions3(i-1,6) + 1))
                    
                        if (abs(solutions3(i,1) - solutions3(i-1,1)) > 1 |...
                            abs(solutions3(i,2) - solutions3(i-1,2)) > 1 | ...
                            abs(solutions3(i,3) - solutions3(i-1,3)) > 1 |  ...
                            abs(solutions3(i,4) - solutions3(i-1,4)) > 1 |  ...
                            abs(solutions3(i,5) - solutions3(i-1,5)) > 1 |  ...
                            abs(solutions3(i,6) - solutions3(i-1,6)) > 1)
                    
                            count4 = count4 + 1;
                            solutions4(count4,:) = solutions3(i,:);
                        else
                            duplicated(c, a) = 1;
                        end
                    else
                        duplicated(c, a) = 1;
                    end
                end
           

                points(c,a) = count4;

                if (count4 >= 1)

                    realmccoy(eigcount: eigcount + D(1) - 1, 1:3) = real(u);  % k1 and k2 values
                    realmccoy(eigcount: eigcount + D(1) - 1, 4) = real(v);    % number of fixed points
                    realmccoy(eigcount: eigcount + D(1) - 1, 5:10) = real(solutions4(:,1:6));  % fixed points
                    eigcount = eigcount + D(1);
                end
            end
        end
    end
end
       

% plotting k1 and k2 vs. number of points

D_points = size(points);

for i = 1:(D_points(1))
    for j = 1:(D_points(2))
        if (points(i,j) == 1)
            points_adj(i,j) = 1;
        elseif(points(i,j) == 2)
            points_adj(i,j) = 2;
        elseif(points(i,j) == 3)
            points_adj(i,j) = 3;
        elseif(points(i,j) > 3)
            points_adj(i,j) = 4;
        end
    end
end


figure(3)
Dose_resp = gca;
set(gcf, 'Color', [1 1 1])

contour(log10(k2_vector), log10(k1_vector), points_adj)
xlabel('k2')
ylabel('k1')
zlabel('Number of fixed points')
colorbar

print -f3 -r200 -depsc -tiff /Users/Dao/Desktop/Lyap_files/fixed_points_v9

