%%
% Matlab code in this directory computes the results published in:
%    Heeger DJ, Mackey WE, Oscillatory Recurrent Gated Neural Integrator
%    Circuits (ORGaNICs): A Unifying Theoretical Framework for Neural
%    Dynamics, Proceedings of the National Academy of Sciences, 2019.
%
% See also:
%    Heeger DJ, Mackey WE, ORGaNICs: A Canonical Neural Circuit Computation,
%    bioRxiv 506337; doi: https://doi.org/10.1101/506337, 2018.
% and:
%    Heeger DJ, Mackey WE, ORGaNICs: A Theory of Working Memory in Brains
%    and Machines, arXiv:1803.06288, 2018.

%
% The Matlab code in this directory was developed and written by David
% Heeger (david.heeger@nyu.edu).
% 
% Copyright (c) 2017-2019, David J. Heeger, New York University

%%
% To run the code, you must edit the first two lines of 'makeAllFigures.m'
% so that this directory and its subdirectories get added your matlab path.
% It should look using something like the following:

chdir '~/Matlab/organics'
addpath(genpath('drawPublishAxis'));

%% Flow of computation

% To recompute all of the figures in the published paper and in the SI
% Appendix, simply run the 'makeAllFigures.m' script.

% If you edit that script, you will see that it calls in sequence scripts
% that compute the results for each of the figures and supplementary
% figures:
%
% makeFigure2_3;
% makeFigure4_S5;
% makeFigure5;
% makeFigure6;
% makeFigure7;
% makeFigureS1_S2_S3;
% makeFigureS4;

% The main functions for computing ORGaNICs are:
%
% organics.m: single modulator term (Eq. 1 in the published paper) but
% without normalization: f(y)=0 in Eq. 2.
%
% organicsNorm.m: single modulator with normalization.
%
% organics2.m: multiple modulators (Eq. 6 in the published paper) but
% without normalization.
%
% organicsNorm2.m: multiple modulators with normalization.
%
% The figures in the PNAS paper show simulations with normalization and the
% figures in the bioRxiv paper show simulations without normalization.

% Each of the makeFigureX scripts set some parameters and then execute a
% corresponding script to run the simulation:
%
% 'makeFigure2_3.m' and 'makeFigureS4.m' each call 'odrNormScript.m' or
% optionally 'odrScript.m' (without normalization).
%
% 'makeFigure4_S5.m' calls 'sequentialNormScript.m' or optionally
% 'sequentialScript.m' (without normalization).
%
% 'makeFigure5.m' calls 'premotorNormScript1.m' or optionally
% 'premotorScript1.m' (without normalization).
%
% 'makeFigure6.m' calls 'spatialUpdatingNormScript.m' or optionally
% 'spatialUpdatingScript.m' (without normalization).
%
% 'makeFigure7.m' calls 'premotorNormScript2.m' or optionally
% 'premotorScript2.m' (without normalization).

% Each of the following calls the organics.m function to run the
% simulations without normalization:
%
% odrScript.m
% sequentialScript.m
% premotorScript1.m

% Each of the following calls the organicsNorm.m function:
%
% odrNormScript.m
% sequentialNormScript.m
% premotorNormScript1.m

% Each of the following calls the organics2.m function to run the
% simulations with multiple modulator terms (see Eq. 6 of the manuscript),
% but without normalization:
%
% spatialUpdatingScript.m
% premotorScript2.m

% Each of the following scripts the organicsNorm2.m function to run the
% simulations with multiple modulator terms and with normalization:
%
% spatialUpdatingNormScript.m
% premotorNormScript2.m

%% Variants

% This code includes several variants of ORGaNICs, and you can easily
% modify the code to run each of the variants.

% For example, open 'makeFigure2_3.m' and look at the last two lines:

%odrScript;
odrNormScript;

% Runnig the script as is creates figures that include recurrent
% normalization (Fig 2 and Fig 3b in the published paper).

% Flip the comments:

odrScript;
%odrNormScript;

% Then run the script again. Now it computes the responses without
% normalization (including Fig 3a; see also our bioRxiv paper cited above).

% Likewise, you can compute versions of each of the other figures either
% with or without normalization included by changing the comments near the
% end of each of the other 'makeFigureX' scripts:

% In addition, this code includes different recurrent circuit models of
% normalization. For example, open 'odrNormScript.m' and look at lines
% 173-175 (similarly, lines 290-292):

[y, z, a, b, u] = organicsNorm(x, deltaT, tauA, tauB, tauY, tauU, Wyy, cy, Wzx, cz, Wax, Wbx, Wby, cb, sigma);
%[y, z, a, b, u, v] = organicsNormV(x, deltaT, tauA, tauB, tauY, tauU, tauV, Wyy, cy, Wzx, cz, Wax, Wbx, Wby, cb, sigma);
%[y, z, a, b, u] = organicsNormFlipA(x, deltaT, tauA, tauB, tauY, tauU, Wyy, cy, Wzx, cz, Wax, Wbx, Wby, cb, sigma);

% The function 'organicsNorm.m' was used to compute the figures in the
% published paper and corresponds to Eqs. S26-S29 in the SI Appendix
% (reproduced in 'organicsNormEquations.pdf' in this directory). The
% function 'organicsNormV.m' is an alternative expressed by Eqs. S44-S48 in
% the SI Appendix (reproduced in 'organicsNormVEquations.pdf' in this
% directory). The function 'organicsNormFlipA.m' is yet another alternative
% expressed by Eqs. S49-S52 in the SI Appendix (reproduced in
% 'organicsNormFlipAEquations.pdf' in this directory).

% Likewise for lines 156-158 and 301-303 'sequentialNormScript.m' which is
% called by 'makeFigure4_S5.m'. And likewise for lines 104-106 of
% 'premotorScript1.m' which is called by 'makeFigure5.m'

%% Parameter values

% You can vary the parameters of the simulations. For example, open
% 'makeFigure2_3.m' and compare it with 'makeFigureS4.m'. Line 15 of
% 'makeFigures4.m' specifies the non-zero oscillation frequencies. In
% addition, lines 18-19 of 'makeFigures2_3.m' specifies the amplitude of
% the input drive. Try changing the input drive amplitude to be 0.05 and
% compare the simulation results using 'organicsNorm.m' versus
% 'organicsNormV.m'.

% Likewise for each of the other figures.

% Open 'odrNormScript.m' and look at lines 8-16. This code specifies the
% time constants (the values of the tau values in the equations) and the
% semi-saturation constant (sigma in the equations). Line 129 specifies the
% temporal sampling rate (deltaT).

% Likewise for each of the other scripts:
% 'odrScript.m' (called by 'makeFigure2_3.m')
% 'odrNormScript.m' (called by 'makeFigure2_3.m')
% 'sequentialScript.m' (called by 'makeFigure4_S5.m')
% 'sequentialNormScript.m' (called by 'makeFigure4_S5.m')
% 'premotorScript1.m' (called by 'makeFigure5.m')
% 'premotorNormScript1.m' (called by 'makeFigure5.m')
% 'spatialUpdatingScript.m' (called by 'makeFigure6.m')
% 'spatialUpdatingNormScript.m' (called by 'makeFigure6.m')
% 'premotorScript2.m' (called by 'makeFigure7.m')
% 'premotorNormScript2.m' (called by 'makeFigure7.m')

% Some of the parameters are set to values such that the code only
% approximates the final figures in the published paper. For example, line
% 62 of 'premotorScript1.m' sets the temporal sample rate equal to 0.1 msec
% but a alue of 0.01 msec is needed to get accurate results. Of course,
% that takes 10x longer to run.

%% Stability

% We point out in the paper that, "The stability of the normalized
% responses did not depend on fine tuning any of the other synaptic weights
% in the circuit; perturbing those synaptic weights by random values within
% ±5% yielded virtually identical simulated responses and the responses
% were stable even when those synaptic weights were perturbed by random
% values ranging from 0.5x to 2x." These results are not shown in the
% figures but can be readily reproduced. Open 'organicsNorm.m' and look at
% lines 28-30:

p = 1 + p * (rand(6,N));
id = (rand(6,N)>0.5);
p(id) = 1./p(id);

% The value of p is an optional parameter to the function that can be used
% to perturb all the sypases other than the recurrent drive. Line 8 of this
% function sets % Perturb all the sypases other than the recurrent drivethe
% % default value for p. Change that line to be a non-zero number. For
% example, p=2 would multiple each synaptic weight by a randam value
% between 0.5 and 2. Likewise, p=10 would multiple each synaptic weight by
% a randam value between 0.1 and 10. Instead of editing line 8 of this
% funciton, you can instead call it with a value for p as the last argument
% in the function call.