clear all
close all

rng(1000); % seed 1000

%==========================================================================
%                              Choose and setup model
%==========================================================================

% For Figure 5, linear COI model 
Q1=[ [ -4   4/1e5             0 ]; ...    
     [  4   -(4/1e5+3)    3/1e5 ];
     [  0   3            -3/1e5 ];];

% For Figure S3, triangular COI model 
% Q1=[ [ -6   3/1e5      3/1e10          ]; ...    
%      [  3  -(3+3/1e5)  3/1e5           ]; ... 
%      [  3   3         -(3/1e5+3/1e10)  ];];
 
 
Q2=Q1; % identical channels

% Define labels
l=char({'C'; 'O'; 'I';}); % the labels for "closed" "open" and "inactivated"
sz=size(l);
nlabels=sz(1);
for i=1:nlabels 
   eval(['lab.' l(i,:) '=' num2str(i) ';']) % this will execute e.g., "C=1" and "O=2"
end

nstCh1= length(Q1); % number of states of channel 1
nstCh2= length(Q2); % number of states of channel 2

% Current generated by individual states (here 0 or -1)
% Different states may generate the same current, hence these states are
% indistinguishable (aggregated to the same value)
agr1=[0 -1 0]; % 0 means no current, -1 means unitary current is -1 (pA)
agr2=[0 -1 0]; % 0 means no current, -1 means unitary current is -1

% Aggregate according to shut/open property: being open 1 or not open 0
agr1o=[0 1 0]; % 0 means not open (shut), 1 means open
agr2o=[0 1 0];  % 0 means not open (shut), 1 means open

% Initial probabilities (initial conditions)
p01=[1 0 0]'; % all channels closed
p02=[1 0 0]'; % all channels closed

%------------------------

% Simulation parameters for the 2 channel models are assembled into structure "modelParam"
modelParam.matrixCh1=Q1;
modelParam.matrixCh2=Q2;
modelParam.aggrCh1=agr1;
modelParam.aggrCh2=agr2;
modelParam.OaggrCh1=agr1o;
modelParam.OaggrCh2=agr2o;
modelParam.P0Ch1=p01;
modelParam.P0Ch2=p02;
modelParam.nstatesCh1=nstCh1;
modelParam.nstatesCh2=nstCh2;
modelParam.labels=l; % should here be the same labels for both channels

%------------------------

% Further simulation parameters are assembled into structure "simParams"
simParam.Ns=3000; % number of samples = time steps; default 3000 * dt = 3 ms
simParam.Nr=1000; % number of sweeps
simParam.dt=0.001; % time step
simParam.t_selected=1; % time in ms at which to extract certain channel and current observables
simParam.pLevelFisher=0.01; % significance level for Fisher's exact test
simParam.EntropyThr=-0.01; % significance threshold for the entropy difference

%------------------------

% Parameters for the interactions between the channels
% ----------------------------------------------------

% Energy barriers

% No change in free energy of the barriers
bValue=[]; % changes in multiples of kT
bS1=[]; % compound state on one side of the barrier
bS2=[]; % compound state on the other side of the barrier

% % The following code would increse the CO-OO barriers by 2 kT
% bValue=[2]; % changes in multiples of kT
% bS1=[[lab.C lab.O];]; % compound state on one side of the barrier
% bS2=[[lab.O lab.O];]; % compound state on the other side of the barrier

% Assemble barriers energy change into structure "barrier"

barrier.value=bValue;
barrier.b1=bS1;
barrier.b2=bS2;

%------------------------

% Energies of the compound states

% The following code would stand for "no change in the energies of composite states"
% eValue=[]; % changes in multiples of kT
% eS1=[]; % state of 1st channel
% eS2=[]; % state of 2nd channel

eValue=2; % changes in multiples of kT
eS1=[lab.C]; % state of 1st channel
eS2=[lab.O]; % state of 2nd channel

energy.value=eValue;
energy.state1=eS1;
energy.state2=eS2;

% Assemble composite states energy change into structure "energy"

%------------------------

% Assemble plot options into structure "plotOptions"

plotOptions.plotFlag=1; % 0 to not plot, 1 to plot
plotOptions.saveFig=0; % 1 to save figure as *.fig (if plotted), 0 to not save
plotOptions.saveTif=0; % 1 to save figure as TIF bitmap (if plotted), 0 to not save
plotOptions.delFig=0; % 0 to leave figure at the end of the function call below, 1 to delete (useful if the function is called repeatedly by a loop)

% Special plot settings for COI linear model
plotOptions.po_scale=[0 0.5]; % y-range for the open probabilities of the 2 channels
plotOptions.i_scale=[-1 0]; % y-range for the ensemble average current
plotOptions.p_scale=[0 0.55]; % y-range for the plot of p1, p2, p1bar, p2bar
plotOptions.so_scale=[0 0.5]; % x-range for the histogram of latencies between successive openings
plotOptions.sc_scale=[0 0.5]; % x-range for the histogram of latencies between successive closings
plotOptions.ds_scale=[-0.2 0.01]; % y-range for the entropy difference scale
plotOptions.p1_scale=[0 0.55]; % x-range for p1 in the p2 vs. p1 plot
plotOptions.p2_scale=[0 0.27]; % y-range for p2 in the p2 vs. p1 plot

% Special plot settings for COI triangular model 
% plotOptions.po_scale=[0 0.3]; % y-range for the open probabilities of the 2 channels
% plotOptions.i_scale=[-0.5 0]; % y-range for the ensemble average current
% plotOptions.p_scale=[0 0.4]; % y-range for the plot of p1, p2, p1bar, p2bar
% plotOptions.so_scale=[0 0.5]; % x-range for the histogram of latencies between successive openings
% plotOptions.sc_scale=[0 0.5]; % x-range for the histogram of latencies between successive closings
% plotOptions.ds_scale=[-0.2 0.01]; % y-range for the entropy difference scale
% plotOptions.p1_scale=[0 0.4]; % x-range for p1 in the p2 vs. p1 plot
% plotOptions.p2_scale=[0 0.1]; % y-range for p2 in the p2 vs. p1 plot

%return

% Call the function using all the parameters defined above
% For the description of N,I,F, see the function SymmInteractionAnalysisFigure

[N,I,F]=SymmInteractionAnalysisFigure(modelParam,simParam,barrier,energy,plotOptions);

