% MolterNetwork.m
% A program that models the hippocampal network proposed by Colin Molter et al
% for space learning in rats.

% By Zachary Glass
% for BioNB222 Computational Section
% Professor Bruce Land
% Tuesdays 10:30-12:00, Mudd 325

% Coding Record
% Begun: April 10, 2008
% Edited: April 14, 2008
%         April 21, 2008
%         April 27, 2008
%         April 28, 2008
%         May   01, 2008

clc
hold off
% =========================================================================
% =========================================================================

% Parameter Definitions
global Alpha R Gamma n dimensions time OmegaNot OmegaS Beta_EC w33 fi Beta_CA3
global plotColorArray
OmegaNot = 1; 	% OmegaNot represents the underlying theta rhythm
OmegaS = 0.6;	% OmegaS represents the phase precession due to increase in firing
                %   while in a place cell's receptive field
Alpha = 1;      % Alpha is a constant that couples the two local field potentials
Beta_EC = 1.5;	% Beta_EC 
Beta_CA3 = 1.15;% Beta_Ca3
%time = 25;      % time is the number of timesteps to run
%n = 3;          % n is the number of inputs, EC cells, and CA3 cells
Knot_EC  = 1.5;	% Knot_EC is the influence of EC LFP on on EC cells
Knot_CA3 = 0.2; % Knot_CA3 is the influence of CA3 LFP on CA3 cells
delta_EC  = 0;	% threshold for output of EC cells
delta_CA3 = 0.9;% threshold for output of CA3 cells
K_EC = 3;       % K_EC is the influence of EC output on CA3 cells
k33 = 0;        % k33 regulates how recurrent connections influence CA3 units
fi = -6;        % fi is an in inhibition factor
R = 2;          % R is a saturation value for the Hebbian plasticity
Gamma = 0.01;   % Gamma is also reflects Hebbian plasticity
Tau = 1;        % Tau is a Hebbian asymmetric time constant
% randomType = 0; % randomType lets the user determine the kind of SPW input
% dimensions = 0; % dimensions lets you determine if the model is 1D or 2D

% =========================================================================
% LEARNING PHASE
% =========================================================================

% INPUTS

% Here we define the inputs into the system.
% There are n cells in the EC, so there must be n inputs.
% A different current can be input into each cell at each time, so
% the input matrix has a row for each cell and a column for each time step

% A function molterInputCalculator has been written to automatically
% generate the input matrix.  It assumes a computational time of 10 seconds
% per node of the model.  Currently, it can only handle a 1D model, but
% will eventually include 2 dimensional models

n = input('Input the number of nodes. ');
dimensions = input(' Is the model 1D (1) or 2D (2)? ');
time = 10*n;
I = molterInputCalculator;
I = 1.5 * I;    % Each input fires for 10 seconds, with an amplitude of 1.5

% *********************************************************************

% LOCAL FIELD POTENTIALS

% The brain has a constant rhythm.  In this model, it is represented by the LFP.
% There is a different LFP for each level of cells, but they are linked.

% Phinot_EC  is the LFP in the EC  level
% Phinot_CA3 is the LFP in the CA3 level

Phinot_EC  = OmegaNot * ones(1,time); % Initialize Phinot_EC
Phinot_CA3 = zeros(1,time);

% Phinot_CA3 is iterative, so we must define it at each time step.  First step must be
% defined alone, then the rest can be iterated

Phinot_CA3(1) = OmegaNot;
for ii = 2:time
	Phinot_CA3(ii) = OmegaNot - ((Alpha * cos(Phinot_EC(ii-1))) .* sin(Phinot_CA3(ii-1)));
end

%  The local field potentials should phase-lock.  We will plot their phase
%  to make certain that they do
figure(1)
clf(1)
title ('Local Field Potentials')
plot(1:1:time,Phinot_EC,'r*-',1:1:time,Phinot_CA3,'bx-')
legend('EC LFP','CA3 LFP')

% *********************************************************************

% NATIVE ANGULAR FREQUENCY

% Each unit has a native angular frequency at which it fires.  For the CA3 layer, it is
% assumed to be constant at OmegaNot.  For the EC layer, it varies as a place cell
% remains stimulated

% CA3 layer can be initialized easily

Omega_CA3 = OmegaNot * ones(n,time);
Omega_EH  = OmegaNot * ones(n,time);

%%% FOR NOW, the change in native frequency in EH will be artificial
% Omega_EH(1,1:10)  = (OmegaNot - OmegaS) : (2* OmegaS / 9) : (OmegaNot + OmegaS);
% Omega_EH(2,8:17)  = (OmegaNot - OmegaS) : (2* OmegaS / 9) : (OmegaNot + OmegaS);
% Omega_EH(3,15:24) = (OmegaNot - OmegaS) : (2* OmegaS / 9) : (OmegaNot + OmegaS);

% A function was written to calculate  this in a less artificial manner.
Omega_EH = molterOmegaECCalculator;

% *********************************************************************

% INDIVIDUAL UNIT INPUT, PHASE, AND OUTPUT

% Each individual unit will receive input, adjust its own phase based on this input, and 
% send along output based on its phase.  This will happen for each time step.  For sake of
% clarity, we will store the EC layer and the CA3 layer separately, and store Input, Phase, 
% and Output separately.  First, we will initalize all six, for data storage purposes.

Input_EC   = zeros(n,time);
Phase_EC   = zeros(n,time);
Output_EC  = zeros(n,time);

Input_CA3  = zeros(n,time);
Phase_CA3  = zeros(n,time);
Output_CA3 = zeros(n,time);

% The synaptic strengths in this model are plastic.  We will now define a
% weight vector w33, which is a square matrix repesenting the synaptic
% weight of each cell in the CA3 layer with each other CA3 cell.  Also,
% there is a matrix of various inputs, Gi that affects each cell at all
% times.  Furthermore, the synaptic weights will be updated using a simple
% derivative, callled wdot

w33 = ones(n,n);
wdot33 = ones(n,n);
Gi  = zeros(n,time);

% For the EC region, the inputs do not vary with time, but vary with each cell.
% This can be calculated using vectorized math, looping through all cells

for ii = 1:n
	Input_EC(ii,:) = (Knot_EC .* cos(Phinot_EC)) + I(ii,:);
end

% Next we can calculate the phase and output of the EC region, as an iteration of time.
% We first define the initals, and then loop.
% Notice that the output depends on a threshold value delta

Phase_EC(:,1) = Omega_EH(:,1);
Output_EC(:,1) = ( cos(Phase_EC(:,1)) >= delta_EC ) .* ((cos(Phase_EC(:,1)) - delta_EC) / (1 - delta_EC) );

for ii = 2:time
	Phase_EC(:,ii) = Omega_EH(:,ii) + (Beta_EC - Input_EC(:,ii)) .* sin(Phase_EC(:,ii-1));
	Output_EC(:,ii) = ( cos(Phase_EC(:,ii)) >= delta_EC ) .* ((cos(Phase_EC(:,ii)) - delta_EC) / (1 - delta_EC) );
end


% Next we need to calculate the Input, Phase, and Output of each of the CA3 cells.  All of these take
% iteration, so we must define the first time step, and then loop the
% remaining steps.  I will assume no CA3-CA3 inputs in the first time step

% $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
% $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
% ASSSUME THIS IS WRONG!!! BELOW IS "CORRECTED"
% % Input_CA3(:,1)  = Knot_CA3 * cos(Phinot_CA3(1)) - (K_EC .* (Output_EC(:,1))) + Gi(:,1);
% % Phase_CA3(:,1)  = Omega_CA3(:,1);
% % Output_CA3(:,1) = ( cos(Phase_CA3(:,1)) >= delta_CA3 ) .* ((cos(Phase_CA3(:,1)) - delta_CA3) / (1 - delta_CA3) );
% % 
% % % Now loop for time
% % 
% % for ii = 2:time
% %     Input_CA3(:,ii)  = Knot_CA3 * cos(Phinot_CA3(ii)) - K_EC .* (Output_EC(:,ii)) -...
% %         (k33 * (w33 * Output_CA3(:,ii-1))) - (fi * sum(Output_CA3(:,ii-1))) + Gi(:,ii);
% %     Phase_CA3(:,ii)  = Omega_CA3(:,ii) + ((Beta_CA3 - Input_CA3(:,ii)) .* sin(Phase_CA3(:,ii-1)));
% %     Output_CA3(:,ii) = ( cos(Phase_CA3(:,ii)) >= delta_CA3 ) .* ((cos(Phase_CA3(:,ii)) - delta_CA3) / (1 - delta_CA3) );
% %     % For now, I will calculate wdot33 using a separate function
% %     wdot33 = wdot33Calculator (Output_CA3(:,ii),Output_CA3(:,ii-Tau), wdot33);
% %     w33 = w33 + wdot33;
% %     %
% % end

Input_CA3(:,1)  = Knot_CA3 * cos(Phinot_CA3(1)) + (K_EC .* (Output_EC(:,1))) + Gi(:,1);
Phase_CA3(:,1)  = Omega_CA3(:,1);
Output_CA3(:,1) = ( cos(Phase_CA3(:,1)) >= delta_CA3 ) .* ((cos(Phase_CA3(:,1)) - delta_CA3) / (1 - delta_CA3) );

% Now loop for time
choice = input('Use Molter''s equation (1) or Zach''s equation (2) for wdot33? ');

for ii = 2:time
    Input_CA3(:,ii)  = Knot_CA3 * cos(Phinot_CA3(ii)) + K_EC .* (Output_EC(:,ii)) +...
        (k33 * (w33 * Output_CA3(:,ii-1))) - (fi * sum(Output_CA3(:,ii-1))) + Gi(:,ii);
    Phase_CA3(:,ii)  = Omega_CA3(:,ii) + ((Beta_CA3 - Input_CA3(:,ii)) .* sin(Phase_CA3(:,ii-1)));
    Output_CA3(:,ii) = ( cos(Phase_CA3(:,ii)) >= delta_CA3 ) .* ((cos(Phase_CA3(:,ii)) - delta_CA3) / (1 - delta_CA3) );
    % For now, I will calculate wdot33 using a separate function
    wdot33 = wdot33Calculator (Output_CA3(:,ii),Output_CA3(:,ii-Tau), wdot33,choice);
    w33 = w33 + wdot33;
    
end

% $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
% $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
% Now, we will do some preliminary displays
w33

plotColorArray = ['cs-';'ro-'; 'bx-'; 'g+-'; 'yp-'; 'kd-'];
figure(2)
clf(2)
subplot(2,1,1)
hold on

for ii = 1:n
    plotCount = (mod(ii,6) + 1);
    plot((1:1:time),Output_EC(ii,:),plotColorArray(plotCount))
    hold on
end
title ('Output of the EH layer with time, Learning Phase')
xlabel('time,t')
ylabel('Output')

subplot(2,1,2)
hold on
for ii = 1:n
    plotCount = (mod(ii,6)+1);
    plot((1:1:time),Output_CA3(ii,:),plotColorArray(plotCount))
    hold on
end
title ('Output of the CA3 layer with time, Learning Phase')
xlabel('time,t')
ylabel('Output')

figure(3)
clf(3)
RastorPlot(Output_EC,0.1,3)
title('Rastor Plot of the EC layer, Learning Phase')
figure(4)
clf(4)
RastorPlot(Output_CA3,0.1,4)
title('Rastor Plot of the CA3 layer, Learning Phase')


% =========================================================================
% REACTIVATION PHASE
% =========================================================================

% ReactivationOutputs = MolterReactivationPhase;

% FIRST, redefine parameters that change from learning to reactivation
Knot_EC = 0;
Knot_CA3 = 0;
K33 = 2.5;

% allow delta_CA3 to be independently modified for reactivation
delta_CA3 = 0.9;



% Next, determine type of "random" input.  For now, the random inputs are
% actually manufactured.  this will be fixed later.

randomType = input('1 for sleep, 2 for rest ');
if randomType == 1
    Gi=(0.15 + (0.2 *randn(n,time)));
elseif randomType == 2
    Gi=(0.15 + (0.2 *randn(n,time)));
    Gi(n-1:n,:) = (0.5 + (0.2 *randn(2,time)));
end

% Now, simulate CA3 phase for a second lenght of time
time2 = time;
Input_CA3_prime  = zeros(n,time2);
Phase_CA3_prime  = zeros(n,time2);
Output_CA3_prime = zeros(n,time2);
%FOR NOW, artificially define Output_EC_prime
Output_EC_prime = zeros(n,time2);
Output_EC_prime(:,1) = cos(OmegaNot);
for ii = 2:time2
    Output_EC_prime(:,ii) = cos( (OmegaNot + (Beta_EC * sin(acos(Output_EC_prime(:,ii-1))))));
end

% and copy/paste from above

% *************** this is copy/pasted

for ii = 2:time2
    Input_CA3_prime(:,ii)  = Knot_CA3 * cos(Phinot_CA3(ii)) + K_EC .* (Output_EC_prime(:,ii)) +...
        (k33 * (w33 * Output_CA3_prime(:,ii-1))) - (fi * sum(Output_CA3_prime(:,ii-1))) + Gi(:,ii);
    Phase_CA3_prime(:,ii)  = Omega_CA3(:,ii) + ((Beta_CA3 - Input_CA3_prime(:,ii)) .* sin(Phase_CA3_prime(:,ii-1)));
    Output_CA3_prime(:,ii) = ( cos(Phase_CA3_prime(:,ii)) >= delta_CA3 ) .* ((cos(Phase_CA3_prime(:,ii)) - delta_CA3) / (1 - delta_CA3) );
    % For now, I will calculate wdot33 using a separate function
    wdot33 = wdot33Calculator (Output_CA3_prime(:,ii),Output_CA3_prime(:,ii-Tau), wdot33,choice);
    w33 = w33 + wdot33;
    %
end
% ******************** End Copy/paste

% Now lets display the output of CA3 in this phase
% % figure(5)
% % clf(5)
% % hold on
% % for ii = 1:n
% %     plotCount = ( mod(ii,6)+1);
% %     plot((1:1:time2),Output_CA3_prime(ii,:),plotColorArray(plotCount))
% %     hold on
% % end
% % title ('Output of the CA3 layer with time')
% % xlabel('time,t')
% % ylabel('Output')

figure(6)
clf(6)
RastorPlot(Output_CA3_prime,0.1,6)
title('Rastor Plot of the CA3 layer, Reactivation Phase')

figure(5)
clf(5)
subplot(2,1,1)
hold on

for ii = 1:n
    plotCount = (mod(ii,6) + 1);
    plot((1:1:time),Output_EC_prime(ii,:),plotColorArray(plotCount))
    hold on
end
title ('Output of the EH layer with time, Reactivation Phase')
xlabel('time,t')
ylabel('Output')

subplot(2,1,2)
hold on
for ii = 1:n
    plotCount = (mod(ii,6)+1);
    plot((1:1:time),Output_CA3_prime(ii,:),plotColorArray(plotCount))
    hold on
end
title ('Output of the CA3 layer with time, Reactivation Phase')
xlabel('time,t')
ylabel('Output')

% LastDitch way to determine sequentials
LastDitch = SequenceFinder(Output_CA3_prime);
figure(8)
clf(8)
RastorPlot(LastDitch,0.1,8)
% % 
% % % FOR NOW this is where the program ends.  It currently does not include
% % % synapse plasticity, and it does not simulate the "asleep" phase of the
% % % model.
% % 
% % % For now, we will just graph outputs with time
% % 
% % % % % % % % % % % figure(1)
% % % % % % % % % % % 
% % % % % % % % % % % subplot(2,1,1)
% % % % % % % % % % % hold on
% % % % % % % % % % % for ii = 1:n
% % % % % % % % % % %     plot((1:1:time),Output_EC(ii,:))
% % % % % % % % % % %     hold on
% % % % % % % % % % % end
% % % % % % % % % % % title ('Output of the EH layer with time')
% % % % % % % % % % % xlabel('time,t')
% % % % % % % % % % % ylabel('Output')
% % % % % % % % % % % 
% % % % % % % % % % % subplot(2,1,2)
% % % % % % % % % % % hold on
% % % % % % % % % % % for ii = 1:n
% % % % % % % % % % %     plot((1:1:time),Output_CA3(ii,:))
% % % % % % % % % % %     hold on
% % % % % % % % % % % end
% % % % % % % % % % % title ('Output of the CA3 layer with time')
% % % % % % % % % % % xlabel('time,t')
% % % % % % % % % % % ylabel('Output')
