%Coded by Cyril Tang, Team EchiDNA 2014 for BIOMOD
%Based on work by Bai et al. See DOI: 10.1126/science.1182105.

clear all

canvas=figure(1);
set(canvas, 'Position', [10,10,1000,700]);

%Binding is associated with activation
%i.e. b favours a, but B favours A

%---------------------------------------------
n = 30; %number of protomers
Ea = 1; %Coupling of binding to activation %3
Ej = 4; %Coupling of neighbouring protomers %3
c = 2; %as a proportion of c0.5 (conc req for neutral bias): c/c(0.5)

omegaa = 10^4;
omegab = 10;
lambdaa = 0.5;
lambdab = 0;
%---------------------------------------------

%Graphics parameter initialisation
radius = n/(2*pi);
subplot(1,2,1)
for x = 1:n
    graphicsobjects(x,1) = rectangle('Position',[radius*cos(x*2*pi/n),radius*sin(x*2*pi/n),1,1],'Curvature',[1,1],'FaceColor','w');
    graphicsobjects(x,2) = rectangle('Position',[radius*cos(x*2*pi/n)+0.375,radius*sin(x*2*pi/n)+0.375,0.25,0.25],'Curvature',[1,1],'FaceColor','w');
end
axis square
axis off
    
kbT = 1;
history = 1;

%Columns are: 1) a state (active/inactive as 1/0)
%   2) b state (bound/unbound as 1/0)
%   3) A time as double
%   4) B time as double
state = zeros(n,4);

%Rate constants for active/inactive and bound/unbound
ka = zeros(1,n);
kb = zeros(1,n);

%Binding of signal molecule is independent of protomer conformation
bB = c*omegab;
Bb = @(z) omegab*exp((-log(c)+z*Ea)/kbT); %z = +1 if favourable change, -1 if unfavourable change

%Conformation change of protomer depends on neighbours
ij = @(z) omegaa*exp(-z/(2*kbT));

%Initialise variables with starting values for a ring of "0"s

for x = 1:n
    
    %Assign ka

    ka(x) = ij(Ea + 2*Ej);
           
    %Assign kb
    
    kb(x) = Bb(-1);

    state(x,3) = -log(rand)/ka(x);
    state(x,4) = -log(rand)/kb(x);
            
end

%-------------------------------------------------------
%Here is the meat of the program:
%-------------------------------------------------------
timestep = 0;
while timestep < 10%true
    %Find next timepoint to execute and recalculate self
    if min(state(:,3)) < min(state(:,4))
        [timestep, stepindex] = min(state(:,3));
        steptype = 3;
    else
        [timestep, stepindex] = min(state(:,4));
        steptype = 4;
    end

    %Execute the event
    state(stepindex,steptype-2) = not(state(stepindex,steptype-2));
    subplot(1,2,1)
    if state(stepindex,steptype-2) == 1
        if steptype == 3
            set(graphicsobjects(stepindex, steptype-2),'FaceColor','b');
        else
            set(graphicsobjects(stepindex, steptype-2),'FaceColor','r');
        end
    else
        set(graphicsobjects(stepindex, steptype-2),'FaceColor','w');
    end

    %Recalculate neighbours
    leftstep = stepindex - 1;
    if leftstep == 0
        leftstep = n;
    end
    rightstep = stepindex + 1;
    if rightstep > n
        rightstep = 1;
    end    
    
    neighbours = [leftstep, stepindex, rightstep]; 
    
    for x = neighbours
        if x == 1
             tempa = [state(n,1), state(x,1), state(x+1,1)];
        elseif x == n
             tempa = [state(x-1,1), state(x,1), state(1,1)];
        else
             tempa = [state(x-1,1), state(x,1), state(x+1,1)];
        end

        switch sum(tempa)
            case {0,3}
                if state(x,1) == state(x,2) %+Ea if changed
                    ka(x) = ij(Ea + 2*Ej);
                else %-Ea if changed
                    ka(x) = ij(-Ea + 2*Ej);
                end
            otherwise
                if isequal(tempa,[0,1,0]) || isequal(tempa,[1,0,1])
                    if state(x,1) == state(x,2) %+Ea if changed
                        ka(x) = ij(Ea - 2*Ej);
                    else %-Ea if changed
                        ka(x) = ij(-Ea - 2*Ej);
                    end
                else
                    if state(x,1) == state(x,2) %+Ea if changed
                        ka(x) = ij(Ea);
                    else %-Ea if changed
                        ka(x) = ij(-Ea);
                    end
                end

        end

        %Assign kb
        if state(x,2) == 0
            kb(x) = bB;
        elseif state(x,2) == 1
            if state(x,1) == 0
                kb(x) = Bb(1); %favourable
            else
                kb(x) = Bb(-1); %unfavourable
            end
        end

        %Generate An and Bn
        state(x,3) = -log(rand)/ka(x) + timestep;
        state(x,4) = -log(rand)/kb(x) + timestep;

    end
    
    record(history,1) = timestep;
    record(history,2) = sum(state(:,1))/n;
    record(history,3) = sum(state(:,2));
    history = history + 1;

    subplot(1,2,2)
    plot(record(:,1),record(:,2))
    title(strcat('Order Parameter: ', num2str(sum(record(:,2))/length(record(:,2))), ' at time ', num2str(timestep)));
    figure(1)
    pause(0.0001)
end