% This MATLAB program is an example of code used in "Interhemispheric
% Dominance Switching in a Neural Network Model for Birdsong", byt
% D. Galvis, W. Wu, R. L. Hyson, F. Johnson, and R. BErtram,
% Journal of Neurophysiology, 120:1186-1197, 2018. 

% The code runs the two-hemisphere zebra finch model with a randomly
% selected motif (based on means that are selected below). It requires
% the Spiker and Burster Classes, which should be in the same directory.

% This is a simplified version of the model that can be used to understand
% the implementation of the model.  There is no medial/lateral HVC, PAm does
% not have a connection to HVC (for atypical transition), Uva does not
% inhibit HVC, PAm and RAm are not mutually inhibitory, 
% and Uva is not mutually inhibitory. None of these things are
% necessary for the model when there are no perturbations to the model
% (such as med/lat ablations, unilateral stimulations, cooling etc).

% This model runs much faster in C++, but MATLAB code is easier to parse.

clear;
rng('shuffle');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%The following parameters need to be %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%specified prior to running the model%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

total_attempts = 2; %Runs of the model
avg_syl_num = 6;  %avg syllable number
avg_gap_length = 10;%avg number of links in HVC gap chains
avg_syl_length = 37;%avg number of links in HVC syl chains
pam_length = 70; % number of links in pam chain
avg_switch_num = 1; %avg number of dom switches per syllable
dom_switcher = 4; %HVC ensemble n - dom_switcher activates the recurrent 
                  %network to initiate a intra-syl dominance switch


dom_frac = 0.90; %the strength of dominant connections from HVC to RA
                 %non-dominant connections take 1-dom_frac
           
% These are all in units of timestep . In this case, each timestep is 1
% ms. Therefore, each burst lasts 6 ms and has a refractory period of 12
% ms etc.

% The conversion of timestep to time is chosen by the user. It must be
% chosen so that each of the following parameters takes an integer value.

total_time = 1200; %Time per Run in timesteps
delays = [3;3];
refractory_period = 12;
burst_length = 6;


%End things to choose
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for attempts = 1:total_attempts

    % Determine number of syllables and lengths for all gaps and syllables for
    % this run
    % Switch number is the number of dominance switches (for each syllable)
    % Seg number is the number of segments per syllable (switch_num+1)
    %The number of ensembles per segments is randomly chosen from the avg
    %syllable length divided by the number of segments in that syllable


    %dom_hem_gap and dom_hem_syl determine the dominant hemisphere for each
    %gap and the first segment of each syllable (the rest of the segments have
    %dominance determined by number of switches).
    %The first column has the dom hem, the second column has the ndom hem
    clear syl_num gap_length syl_length;
    clear dom_hem_gap dom_hem_syl;
    %number of syllables
    syl_num = max(1,poissrnd(avg_syl_num));

    %Allocate zero arrays
    gap_length = zeros(syl_num,1);
    dom_hem_gap = zeros(syl_num,2);
    seg_num = zeros(syl_num,1);
    dom_hem_syl = zeros(syl_num,1);


    for j = 1:syl_num
        gap_length(j) = max(1,poissrnd(avg_gap_length));

        dom_hem_gap(j,1) = (binornd(1,0.5)+1);
        if(dom_hem_gap(j,1)==1);dom_hem_gap(j,2)=2;
        else;dom_hem_gap(j,2)=1;end

        seg_num(j) = poissrnd(avg_switch_num)+1;

        dom_hem_syl(j,1) = (binornd(1,0.5)+1);
        if(dom_hem_syl(j,1)==1);dom_hem_syl(j,2)=2;
        else;dom_hem_syl(j,2)=1;end    
    end


    syl_length = zeros(max(seg_num),syl_num);
    for j = 1:syl_num
        for l = 1:seg_num(j)
            syl_length(j,l) = max(dom_switcher,round(poissrnd...
                               (avg_syl_length)/seg_num(j)));

        end
    end
    % End syl/gap/motif lengths

    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%All bursting and tonic neural ensembles are defined in the following code
%Look at Figures 2 and 3 of Galvis et al. 2018. Each population of
%ensembles is identified there.

    % Define Bursters
    clear HVC_gap;
    for j = 1:syl_num
        for k = 1:gap_length(j)
            HVC_gap(dom_hem_gap(j,1),j,k) = Burster(delays,burst_length,...
                                                    refractory_period,...
                                                    dom_frac);

            HVC_gap(dom_hem_gap(j,2),j,k) = Burster(delays,burst_length,...
                                                    refractory_period,...
                                                    1-dom_frac);                                                                                    
        end
    end
    clear HVC_syl;
    for j = 1:syl_num
        for l = 1:seg_num(j)
            for k = 1:syl_length(j,l)
                
                %If statement switches dominance with segment
                if(mod(l,2)==1)    
                HVC_syl(dom_hem_syl(j,1),j,l,k) = Burster(delays,burst_length,...
                                                        refractory_period,...
                                                        dom_frac);
                HVC_syl(dom_hem_syl(j,2),j,l,k) = Burster(delays,burst_length,...
                                                        refractory_period,...
                                                         1-dom_frac); 
                else
                HVC_syl(dom_hem_syl(j,1),j,l,k) = Burster(delays,burst_length,...
                                                        refractory_period,...
                                                        1-dom_frac);

                HVC_syl(dom_hem_syl(j,2),j,l,k) = Burster(delays,burst_length,...
                                                        refractory_period,...
                                                         dom_frac);       

                end 
                
            end 
       end
    end
    clear PAm;
    for i = 1:pam_length
        PAm(i) = Burster(delays,burst_length,refractory_period,1.0);
    end
    
    %Define Tonic Ensembles
    clear Uva_syl Uva_gap RAm RA_dR;
    for j = 1:syl_num
        Uva_syl(j) = Spiker(delays);
        Uva_gap(j) = Spiker(delays);
        for l = 1:seg_num(j)
            for k = 1:syl_length(j,l)
                RAm(j,l,k) = Spiker(delays);
                for i = 1:2
                    RA_dR(i,j,l,k) = Spiker(delays);
                end
            end
        end
    end
    clear RA_dPI DM_syl DM_gap RA_dPE RA_dD;
    for i = 1:2
        RA_dPI(i) = Spiker(delays);
        for j = 1:syl_num
            DM_syl(i,j) = Spiker(delays);
            DM_gap(i,j) = Spiker(delays);
            RA_dPE(i,j) = Spiker(delays);
            RA_dD(i,j)  = Spiker(delays);
        end
    end
    clear DM_intra RA_intra Uva_intra;
   for j = 1:syl_num
       for l = 1:seg_num(j)
            for i = 1:2
                DM_intra(i,j,l) = Spiker(delays);
                RA_intra(i,j,l) = Spiker(delays);
            end
            Uva_intra(j,l) = Spiker(delays);
        end
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Run the simulation!!!! This has 3 broad steps
%1) Identify the synaptic weights for all ensembles
%2) Run the Step function to update all ensembles
%3) Run the Syn_Step function to update all synapses for the next step



clear RA_plot RAm_plot PAm_plot;
% Allocates arrays to hold the solutions for plotting
HVC_plot = zeros(total_time,2,syl_num,max(seg_num),max(max(syl_length)));
RA_plot = zeros(total_time,2,syl_num,max(seg_num),max(max(syl_length)));
RAm_plot = zeros(total_time,syl_num,max(seg_num),max(max(syl_length)));
PAm_plot = zeros(total_time,pam_length);
% End allocation



   for times = 0:total_time
        if(mod(times,10)==0);disp(num2str(times));end
        
%Synaptic Weight Variables Allocated Here
PAm_exc = zeros(pam_length,1);
PAm_inh = zeros(pam_length,1);
Uva_syl_exc = zeros(syl_num,1);
Uva_syl_inh = zeros(syl_num,1);
Uva_gap_exc = zeros(syl_num,1);
Uva_gap_inh = zeros(syl_num,1);
RAm_exc = zeros(syl_num,max(seg_num),max(max(syl_length)));
RAm_inh = zeros(syl_num,max(seg_num),max(max(syl_length)));
RA_dR_exc = zeros(2,syl_num,max(seg_num),max(max(syl_length)));
RA_dR_inh = zeros(2,syl_num,max(seg_num),max(max(syl_length)));
HVC_syl_exc = zeros(2,syl_num,max(seg_num),max(max(syl_length)));
HVC_syl_inh = zeros(2,syl_num,max(seg_num),max(max(syl_length)));
HVC_gap_exc = zeros(2,syl_num,max(gap_length));
HVC_gap_inh = zeros(2,syl_num,max(gap_length));
RA_dPI_exc = zeros(2);
RA_dPI_inh = zeros(2);
DM_syl_exc = zeros(2,syl_num);
DM_syl_inh = zeros(2,syl_num);
DM_gap_exc = zeros(2,syl_num);
DM_gap_inh = zeros(2,syl_num);
RA_dPE_exc = zeros(2,syl_num);
RA_dPE_inh = zeros(2,syl_num);
RA_dD_exc = zeros(2,syl_num);
RA_dD_inh = zeros(2,syl_num);
DM_intra_exc = zeros(2,syl_num,max(seg_num));
DM_intra_inh = zeros(2,syl_num,max(seg_num));
RA_intra_exc = zeros(2,syl_num,max(seg_num));
RA_intra_inh = zeros(2,syl_num,max(seg_num));
Uva_intra_exc = zeros(syl_num,max(seg_num));
Uva_intra_inh = zeros(syl_num,max(seg_num));        
        
        
        
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%1) Identify the synaptic weights for all ensembles

%NOTE for HVC syl ensembles S_far is weighted to account for the dominance.
%Look in the class Burster details

        %THIS ACTIVATES THE MODEL BY STARTING SYLA AT TIME = 0
        if(times == 0) 
            HVC_syl_exc(1,1,1,1) = 1;
            HVC_syl_exc(2,1,1,1) = 1;
        end
        
        
        %PAm is a chain so each is excited by the previous one
        for i = 2:pam_length
            PAm_exc(i) = PAm_exc(i) + PAm(i-1).S_close;      
        end
        %The first PAm ensemble can be excited by the last ensemble
        % in the last segment of any syllable
        for i = 1:2
            for j = 1:syl_num
                PAm_exc(1) = PAm_exc(1) + RA_dPE(i,j).S_far; 
            end
        end
        % PAm is inhibited by RA_dPI
        for i = 1:2
            for ii = 1:pam_length
                PAm_inh(ii) = PAm_inh(ii) + RA_dPI(i).S_far; 
            end
        end
        % RAm is inhibited by RA_dPE
        for j = 1:syl_num
           for l = 1:seg_num(j)
               for k = 1:syl_length(j,l)
                   for ii = 1:2
                       for jj = 1:syl_num
                            RAm_inh(j,l,k) = RAm_inh(j,l,k) + ...
                                                 RA_dPE(ii,jj).S_far; 
                       end
                   end
               end
           end
        end
        
        % Gap chains
         for i = 1:2
            for j = 1:syl_num
                for k = 2:gap_length(j)
                    HVC_gap_exc(i,j,k)=HVC_gap_exc(i,j,k) + ...
                                       HVC_gap(i,j,k-1).S_close;
                end
            end
         end
        %Syllable chains
        %HVCsyl ->RAdR ->RAm
        for j = 1:syl_num
            for l = 1:seg_num(j)
                for k = 1:syl_length(j,l)
                    for i = 1:2
                        if(k~=1)
                        HVC_syl_exc(i,j,l,k) = HVC_syl_exc(i,j,l,k) + ...
                                                HVC_syl(i,j,l,k-1).S_close;
                        end
                        RA_dR_exc(i,j,l,k) = RA_dR_exc(i,j,l,k) + ...
                                              HVC_syl(i,j,l,k).S_far;
                        RAm_exc(j,l,k) = RAm_exc(j,l,k) + ...
                                            RA_dR(i,j,l,k).S_far; 
                                            
                    end
                end
            end
        end
        %HVC_syl first ensemble in a syllable inhibits PAm through RA_dPI
        for i = 1:2
            for j = 1:syl_num
           RA_dPI_exc(i) = RA_dPI_exc(i) + HVC_syl(i,j,1,1).S_far; 
            end
        end
        
        
        
        % DM gets activated by RA during transition
        % RA gets activated by the last syl or gap ensemble during
        % transition
        % Uva gets activated by DM (in the case of DMgap -> Uvagap,
        % increase the index).
        for i = 1:2
            for j = 1:syl_num
                DM_syl_exc(i,j) = DM_syl_exc(i,j) + RA_dD(i,j).S_far;
                DM_gap_exc(i,j) = DM_gap_exc(i,j) + RA_dPE(i,j).S_far;
                
                RA_dPE_exc(i,j) = RA_dPE_exc(i,j) + HVC_syl(i,j,...
                               seg_num(j),syl_length(j,seg_num(j))).S_far;
                RA_dD_exc(i,j) = RA_dD_exc(i,j) + HVC_gap(i,j,...
                                                  gap_length(j)).S_far;
                                              
                Uva_syl_exc(j) = Uva_syl_exc(j) + DM_syl(i,j).S_far;
                Uva_gap_exc(mod(j,syl_num)+1) = Uva_gap_exc(mod(j,syl_num)+1)...
                                                     + DM_gap(i,j).S_far;
            end
        end
        % HVC gets excited by Uva
        for i = 1:2
            for j = 1:syl_num
                HVC_syl_exc(i,j,1,1) = HVC_syl_exc(i,j,1,1) + ...
                                        Uva_syl(j).S_far;
                HVC_gap_exc(i,j,1) = HVC_gap_exc(i,j,1) + ...
                                        Uva_gap(j).S_far;
            end
        end

        % Intra-syllable transitions (only if there is more than one
        % segment in the syllable
        for j = 1:syl_num
            if(seg_num(j)>1)
               for l = 1:(seg_num(j)-1)
                    for i = 1:2

                       RA_intra_exc(i,j,l) = RA_intra_exc(i,j,l) + ...
                                      HVC_syl(i,j,l,...
                                      syl_length(j,l)-dom_switcher+1).S_far;

                       DM_intra_exc(i,j,l) = DM_intra_exc(i,j,l) + ...
                                             RA_intra(i,j,l).S_far;

                       Uva_intra_exc(j,l) = Uva_intra_exc(j,l) + ...
                                            DM_intra(i,j,l).S_far;
                                        
                                        
                        HVC_syl_exc(i,j,l+1,1) =HVC_syl_exc(i,j,l+1,1) + ...
                                            Uva_intra(j,l).S_far;
                    end

               end
            end
        end


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %2) Run the Step function to update all ensembles




        for i = 1:pam_length
            PAm(i) = PAm(i).Step(times,PAm_exc(i),PAm_inh(i));
        end

        for j = 1:syl_num
            Uva_syl(j) = Uva_syl(j).Step(times,Uva_syl_exc(j),Uva_syl_inh(j));
            Uva_gap(j) = Uva_gap(j).Step(times,Uva_gap_exc(j),Uva_gap_inh(j));
            for l = 1:seg_num(j)
                for k = 1:syl_length(j,l)
                    RAm(j,l,k) = RAm(j,l,k).Step(times,RAm_exc(j,l,k),...
                                                         RAm_inh(j,l,k));
                    for i = 1:2
                        RA_dR(i,j,l,k) = RA_dR(i,j,l,k).Step(times,...
                                     RA_dR_exc(i,j,l,k),RA_dR_inh(i,j,l,k));
                        HVC_syl(i,j,l,k) = HVC_syl(i,j,l,k).Step(times,...
                                   HVC_syl_exc(i,j,l,k),HVC_syl_inh(i,j,l,k));
                    end
                end
            end
        end
        for i = 1:2
            for j = 1:syl_num
                for k = 1:gap_length(j)
                    HVC_gap(i,j,k) =  HVC_gap(i,j,k).Step(times, ...
                                    HVC_gap_exc(i,j,k),HVC_gap_inh(i,j,k));
                end
            end
        end
        for i = 1:2
            RA_dPI(i) = RA_dPI(i).Step(times,RA_dPI_exc(i),RA_dPI_inh(i));
                                   
            for j = 1:syl_num
                DM_syl(i,j) = DM_syl(i,j).Step(times,DM_syl_exc(i,j),...
                                                DM_syl_inh(i,j));
                DM_gap(i,j) = DM_gap(i,j).Step(times,DM_gap_exc(i,j),...
                                                DM_gap_inh(i,j));
                RA_dPE(i,j) = RA_dPE(i,j).Step(times,RA_dPE_exc(i,j),...
                                                RA_dPE_inh(i,j));
                RA_dD(i,j) = RA_dD(i,j).Step(times,RA_dD_exc(i,j),...
                                                RA_dD_inh(i,j));
            end
        end

       for j = 1:syl_num
           for l = 1:seg_num(j)
                for i = 1:2
                    DM_intra(i,j,l) = DM_intra(i,j,l).Step(times,...
                              DM_intra_exc(i,j,l),DM_intra_inh(i,j,l)) ;
                    RA_intra(i,j,l) = RA_intra(i,j,l).Step(times,...
                              RA_intra_exc(i,j,l),RA_intra_inh(i,j,l)) ;
                end
                Uva_intra(j,l) = Uva_intra(j,l).Step(times,...
                                   Uva_intra_exc(j,l),Uva_intra_inh(j,l));
            end
        end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%3) Run the Syn_Step function to update all synapses for the next step
              for i = 1:pam_length
            PAm(i) = PAm(i).Syn_Step(times);
        end

        for j = 1:syl_num
            Uva_syl(j) = Uva_syl(j).Syn_Step(times);
            Uva_gap(j) = Uva_gap(j).Syn_Step(times);
            for l = 1:seg_num(j)
                for k = 1:syl_length(j,l)
                    RAm(j,l,k) = RAm(j,l,k).Syn_Step(times);
                    for i = 1:2
                        RA_dR(i,j,l,k) = RA_dR(i,j,l,k).Syn_Step(times);
                        HVC_syl(i,j,l,k) = HVC_syl(i,j,l,k).Syn_Step(times);
                    end
                end
            end
        end
        for i = 1:2
            for j = 1:syl_num
                for k = 1:gap_length(j)
                    HVC_gap(i,j,k) =  HVC_gap(i,j,k).Syn_Step(times);
                end
            end
        end
        for i = 1:2
            RA_dPI(i) = RA_dPI(i).Syn_Step(times);
                                   
            for j = 1:syl_num
                DM_syl(i,j) = DM_syl(i,j).Syn_Step(times);
                DM_gap(i,j) = DM_gap(i,j).Syn_Step(times);
                RA_dPE(i,j) = RA_dPE(i,j).Syn_Step(times);
                RA_dD(i,j) = RA_dD(i,j).Syn_Step(times);
            end
        end

       for j = 1:syl_num
           for l = 1:seg_num(j)
                for i = 1:2
                    DM_intra(i,j,l) = DM_intra(i,j,l).Syn_Step(times) ;
                    RA_intra(i,j,l) = RA_intra(i,j,l).Syn_Step(times) ;
                end
                Uva_intra(j,l) = Uva_intra(j,l).Syn_Step(times);
            end
        end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% If you want save points to plot, do that here 

 for j = 1:syl_num
       for l = 1:seg_num(j)
           for k = 1:syl_length(j,l)
                RAm_plot(times+1,j,l,k) = RAm(j,l,k).V(RAm(j,l,k).V_now);
                for i = 1:2
                HVC_plot(times+1,i,j,l,k) = HVC_syl(i,j,l,k).V(HVC_syl(i,j,l,k).V_now);
                RA_plot(times+1,i,j,l,k) = RA_dR(i,j,l,k).V(RA_dR(i,j,l,k).V_now);
                end
           end
       end
 end
for i = 1:pam_length
    PAm_plot(times+1,i) = PAm(i).V(PAm(i).V_now);
end



   end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plots
    figure(1);clf;hold all;
    title('RAm');
   for j = 1:syl_num
       for l = 1:seg_num(j)
           for k = 1:syl_length(j,l)
           plot(1:total_time,RAm_plot(1:total_time,j,l,k)-j,'k');
           end
       end
   end
   yticks([]);
   xlabel('timesteps');
   
   figure(2);clf;hold all;
   for i = 1:2
       subplot(1,2,i);hold all;
      if(i==1); title('RA Left');
      else; title('RA Right');end
       for j = 1:syl_num
           for l = 1:seg_num(j)
               for k = 1:syl_length(j,l)
               plot(1:total_time,RA_plot(1:total_time,i,j,l,k)-j,'k');
               end
           end
       end
       yticks([]);
   end
   xlabel('timesteps');
   
      figure(4);clf;hold all;
   for i = 1:2
       subplot(1,2,i);hold all;
      if(i==1); title('HVC Left');
      else; title('HVC Right');end
       for j = 1:syl_num
           for l = 1:seg_num(j)
               for k = 1:syl_length(j,l)
               plot(1:total_time,HVC_plot(1:total_time,i,j,l,k)-j,'k');
               end
           end
       end
       yticks([]);
   end
   xlabel('timesteps');
   
   figure(3);clf;hold all;
   for i = 1:pam_length
    plot(1:total_time,PAm_plot(1:total_time,i)-i,'k');
   end
   title('PAm');
    yticks([]);
   xlabel('timesteps');
   disp('Hit enter for the next run');
   pause;
   
end









