%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is just an example/guide for you to refer to when submitting manuscripts to Frontiers, it is not mandatory to use Frontiers .cls files nor frontiers.tex  %
% This will only generate the Manuscript, the final article will be typeset by Frontiers after acceptance.                                                 %
%                                                                                                                                                         %
% When submitting your files, remember to upload this *tex file, the pdf generated with it, the *bib file (if bibliography is not within the *tex) and all the figures.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Manuscript frontiers
% Quan Zou, GWU, DC
% Jun 04, 2015
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Version 3.1 Generated 2015/22/05 %%%
%%% You will need to have the following packages installed: datetime, fmtcount, etoolbox, fcprefix, which are normally inlcuded in WinEdt. %%%
%%% In http://www.ctan.org/ you can find the packages and how to install them, if necessary. %%%

\documentclass{frontiersSCNS} % for Science, Engineering and Humanities and Social Sciences articles
%\documentclass{frontiersHLTH} % for Health articles
%\documentclass{frontiersFPHY} % for Physics and Applied Mathematics and Statistics articles

%\setcitestyle{square}
\usepackage{url,hyperref,lineno,microtype}
\usepackage[onehalfspacing]{setspace}
\usepackage{amsmath}				% amsmath package, useful for mathematical formulas
\usepackage{amssymb}				% amssymb package, useful for mathematical symbols
\usepackage[normalem]{ulem}
\usepackage{xcolor}
\usepackage{color}
% please place your own definitions here and don't use \def but
% \newcommand{}{}
\newcommand\redsout{\bgroup\markoverwith{\textcolor{red}{\rule[0.5ex]{2pt}{0.4pt}}}\ULon}
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\bea}{\begin{eqnarray}}
\newcommand{\eea}{\end{eqnarray}}
%
\newcommand{\wbar}{${\overline \omega}$}
\newcommand{\wt}{${\overline \omega}_{\mathrm {train}}$}
\newcommand{\wc}{${\overline \omega}_{\mathrm {recon}}$}
\newcommand{\wrr}{${\overline \omega}_{\mathrm {recurr}}$}
\newcommand{\wji}{${\omega}_{\mathrm {ji}}$}
\newcommand{\wffji}{${\omega}^{\mathrm {ff}}_{\mathrm {ji}}$}
\newcommand{\wrcji}{${\omega}^{\mathrm {rc}}_{\mathrm {ji}}$}
%
\newcommand{\rmax}{$r_{\mathrm {max}}$}
\newcommand{\sff}{$s_{\mathrm {ff}}$}
\newcommand{\src}{$s_{\mathrm {rc}}$}
%
\newcommand{\wobsrv}{${\omega}^{\mathrm {observe}}_{\mathrm {ji}}$}
\newcommand{\wpredt}{${\omega}^{\mathrm {predict}}_{\mathrm {ji}}$}
\newcommand{\wbarji}{${\overline \omega}_{\mathrm {ji}}$}
%
\newcommand{\wltp}{$\omega_{\mathrm {LTP}}$}
\newcommand{\wltd}{$\omega_{\mathrm {LTD}}$}
\newcommand{\tltp}{T$_{\mathrm {LTP}}$}
\newcommand{\tltd}{T$_{\mathrm {LTD}}$}
\newcommand{\fltp}{F$_{\mathrm {LTP}}(t)$}
\newcommand{\fltd}{F$_{\mathrm {LTD}}(t)$}
\newcommand{\ap}{A$_{\mathrm p}$}
\newcommand{\aq}{A$_{\mathrm q}$}
\newcommand{\taup}{$\tau_{\mathrm p}$}
\newcommand{\tauq}{$\tau_{\mathrm q}$}
\newcommand{\tauc}{$\tau_{\mathrm {corr}}$}
\newcommand{\et}{\it et al.}
\newcommand{\ie}{\it i.e.,}
\newcommand{\rms}{\it r.m.s.}
\newcommand{\trms}{$\tau_{\mathrm {rms}}$}
\newcommand{\erms}{$r.m.s_{\mathrm {error}}$}
\newcommand{\trmst}{$\tau^{\mathrm {train}}_{\mathrm {rms}}$}
\newcommand{\trmsc}{$\tau^{\mathrm {recon}}_{\mathrm {rms}}$}
\newcommand{\Stdp}{Spike-timing dependent plasticity}
\newcommand{\stdp}{spike-timing dependent plasticity}
\linenumbers

% Leave a blank line between paragraphs instead of using \\


\def\keyFont{\fontsize{8}{11}\helveticabold }
\def\firstAuthorLast{Zou {et~al.}} %use et al only if is more than 1 author
\def\Authors{Quan Zou\,$^{1,*}$, Philip H. Goodman\,$^{2}$ and Frederick C. Harris, Jr.\,$^3$}
% Affiliations should be keyed to the author's name with superscript numbers and be listed as follows: Laboratory, Institute, Department, Organization, City, State abbreviation (USA, Canada, Australia), and Country (without detailed address information such as city zip codes or street names).
% If one of the authors has a change of address, list the new address below the correspondence details using a superscript symbol and use the same symbol to indicate the author in the author list.
\def\Address{$^{1}$Department of Statistics, The George Washington University, Washington, DC 20052, USA \\
$^{2}$School of Medicine and Program in Biomedical Engineering, University of Nevada, Reno, NV 89557, USA\\
$^{3}$Department of Computer Science and Engineering, University of Nevada, Reno, NV 89557, USA}
% The Corresponding Author should be marked with an asterisk
% Provide the exact contact address (this time including street name and city zip code) and email of the corresponding author
\def\corrAuthor{Quan Zou}
\def\corrAddress{my current address}
\def\corrEmail{qzou@gwu.edu}




\begin{document}
\onecolumn
\firstpage{1}

\title[memory encoding and retrieval]{The role of spatiotemporal correlations in the encoding and 
retrieval of synaptic patterns by STDP in recurrent spiking networks} 

\author[\firstAuthorLast ]{\Authors} %This field will be automatically populated
\address{} %This field will be automatically populated
\correspondance{} %This field will be automatically populated

\extraAuth{}% If there are more than 1 corresponding author, comment this line and uncomment the next one.
%\extraAuth{corresponding Author2 \\ Laboratory X2, Institute X2, Department X2, Organization X2, Street X2, City X2 , State XX2 (only USA, Canada and Australia), Zip Code2, X2 Country X2, email2@uni2.edu}


\maketitle

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% The sections below are for reference only.
%%%
%%% For Original Research Articles, Clinical Trial Articles, and Technology Reports the section headings should be those appropriate for your field and the research itself. It is recommended to organize your manuscript in the
%%% following sections or their equivalents for your field:
%%% Abstract, Introduction, Material and Methods, Results, and Discussion.
%%% Please note that the Material and Methods section can be placed in any of the following ways: before Results, before Discussion or after Discussion.
%%%
%%%For information about Clinical Trial Registration, please go to http://www.frontiersin.org/about/AuthorGuidelines#ClinicalTrialRegistration
%%%
%%% For Clinical Case Studies the following sections are mandatory: Abstract, Introduction, Background, Discussion, and Concluding Remarks.
%%%
%%% For all other article types there are no mandatory sections.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{abstract}

%%% Leave the Abstract empty if your article falls under any of the following categories: Editorial Book Review, Commentary, Field Grand Challenge, Opinion or specialty Grand Challenge.
\section{}
%As a primary goal, the abstract should render the general significance and conceptual advance of the work clearly accessible to a broad readership. References should not be cited in the abstract.
% why combine correlation with STDP. 
{\Stdp} (STDP) is considered  % 7
to be an important synaptic mechanism to encode information within and between cerebral  %13
cortical networks. It remains unclear, however, how temporal and spatial % 10
correlations of signals and ongoing background spiking activity interact with STDP % 11
mechanisms to give rise to long-term memory. Because cortical networks mostly  % 11
contain local recurrent connections, it is % 6
critical to explore STDP-based encoding and decoding of synaptic patterns in such architectures.  % 10
In order to investigate their dynamics, we applied a phenomenological STDP % 9
model within recurrent networks, and between recurrent and % 8
feed-forward networks. We find that emergence of transient spatiotemporal correlations  % 10
of ongoing activity lead to the storage of self-organized, simple synaptic patterns in  % 13
the recurrent networks. These networks can be probed at a later time to reconstruct the % 15
encoded patterns by projecting them onto another network in a feed-forward, % 11
readout fashion with highly correlated spatiotemporal structure.  We hypothesize that  % 10
transient spatiotemporal correlations among networks can serve as a biologically % 10
plausible mechanism of memory storage and retrieval based on STDP. % 10


\tiny
 \keyFont{ \section{Keywords:} {\stdp}, {spatiotemporal correlation}, {recurrent network}, {synaptic pattern},
{memory encoding and retrieval}} %All article types: you may provide up to 8 keywords; at least 5 are mandatory.
\end{abstract}

\section{Introduction}\label{intro}

% For Original Research Articles, Clinical Trial Articles, and Technology Reports the introduction should be succinct, with no subheadings.
%
% For Clinical Case Studies the Introduction should include symptoms at presentation, physical exams and lab results.
%
% -------------------------------------------------------------------------- %
%                                                              Introduction  %
% -------------------------------------------------------------------------- %
% general introduction of STDP
The theoretical dependence of excitatory synaptic efficacy upon the relative timing 
of pre-synaptic and post-synaptic action potentials, also referred to as {\stdp} (STDP), 
has been experimentally confirmed (\cite{Mar97,Bgq98,Bel97,Mag97}). Potentiation is typically 
observed when the excitatory post-synaptic potential (EPSP) precedes the post-synaptic spike, whereas 
the reversed temporal order induces depression of the synaptic transmission. A great number of 
computational models have been proposed to account for the dynamics of STDP (see review in~\cite{Mor07}). 
They can be categorized into two basic types: ``hard-bound'' models, in which synaptic weights 
are maintained between minimal and maximal values (\cite{Son00}), 
and ``soft-bound'' models, in which synaptic weights progressively saturate to equilibrium (\cite{Ksl00, 
Ros00, Bdl05, Mor07, Zou07}). These models have been used, for example, to explain cortical remapping 
during development (\cite{Son01}) and cross-model coordinate transformation (\cite{Dav06}). 

% STDP in Recurrent network
Recurrent networks, comprised of ``internal'' links, connecting neurons to one another in an 
undirected manner, may generate their own internal dynamic representations when 
subjected to external stimuli. To perform computation in recurrent networks, STDP can be 
used as the mechanism modifying connection strengths across all synapses.  \cite{Kit02}, 
demonstrated that STDP, when applied to a self-organized recurrent network, can result 
in asynchronous activity at 20$\sim$30 Hz, during delay periods of working memory task. 
It also has been claimed that STDP can initiate resonant or rhythmic activity dominated 
by high spatiotemporal correlations in recurrent networks (\cite{Dau02, Yos07}).
Although this correlated activity has been experimentally observed in many brain regions 
during memory-related tasks, including hippocampal CA3 (\cite{Oke93, Trb89}), it remains unknown 
how correlations regulate and organize the memory within recurrent networks at synaptic level.  

% what we do here
In this paper, we first investigate the mechanism whereby STDP can lead to reliable storage 
and retrieval of simple synaptic patterns within a generic self-organized recurrent network. We then 
show how spatiotemporal correlated activity influences simple pattern formation and completion. 
We demonstrate that the storage process does not {\it require}, but is {\it accelerated} 
by spatially correlated firing among both layers of cells. However, we find that 
at the ``recall'' stage, transferring the stored recurrent pattern into the feedforward 
synapses of a reconstruction layer {\it does} require the presence of spatially correlated 
background activity. We then explore the relationship of learning curves and the time constants 
of correlated activity.  Finally, we analyze the memory capacity of a recurrent network under 
assumptions of differing network size and connection probability.  Our model is consistent with 
the reported literature on both ``hard-bound'' and ``soft-bound'' of STDP dynamics. These results 
suggest a possible generic mechanism of cortical memory organization and retrieval by STDP in 
self-organized recurrent networks.


%\begin{methods}
\section{Material \& Methods}
\subsection{Neuron models}

To investigate the synaptic dynamics, we simulated leaky integrate-and-fire neurons(LIF) 
as in~\cite{Bre07}, where each neuron was described by a single-compartment 
model with a time constant, {$\tau = 20\ ms$}, a leaky conductance of $5\ nS$, and leak reversal 
potential, $E_{L} = -60\ mV$. These parameters corresponded to a $100\ pF$ capacitance.  
The subthreshold membrane potential obeyed the following equation:

\be
	{\tau}\frac{dV}{dt}\ =\ -(V-E_{L})-{\tilde \omega}_{exc}\ (V-E_{exc})-{\tilde \omega}_{inh}\ (V-E_{inh}) %\nonumber 
	\label{eq01} ~.
\ee

Both excitatory and inhibitory synapses were simulated as conductance changes with instantaneous 
jump at maximal value and exponential decay according to the following equation:

\be
	{\tilde \omega}_{exc,inh}(t)\ =\ {\omega}_{exc,inh}\ \times\ e^{-t/{\tau_{exc,inh}}}
	\label{eq02} ~ , 
\ee
where ${\tau_{exc}} = 5\ ms$ and $E_{exc} = 0\ mV$ for excitatory synapses, ${\tau_{inh}} = 10\ ms$, and 
$E_{inh} = -80\ mV$ for inhibitory synapses.  All synapses had the same synaptic delay as $1\ ms$.  
In most situations, the quantal conductance was $5\ nS$ for both feedforward and 
background excitatory connections, $1\ nS$ for recurrent excitatory connections, and 
$50\ nS$ for inhibitory connections, unless specified. The time step of the integration 
is $0.1\ ms$ in all simulations. {We also performed the simulation with shorter time step, 
and there is no qualitative changes in the results.

\subsection{Network model}\label{network}
The network was composed of three layers of neurons: a training layer ($T$), a recurrent layer ($R$) 
and a reconstruction layer ($C$), where each layer contained 100 neurons. Both $T$ and $C$ layers 
maintained all-to-all feedforward excitatory connections to the $R$ layer 
with {\wt} = {\wc} = $5\ nS$}, while R contained local excitatory connections 
at {\wrr} = $1\ nS$ with connection probability of ${\rho} = 100\%$, unless specified. 
In general, we denote the {\wbar} as the upper bound of synaptic weights subject 
to STDP, and $\omega$ as the modulated maximal conductance in Equ.~\ref{eq02}.
In addition, each neuron of $R$ also received a single excitatory background input $E$ (corresponding to 
100 one-to-one connections) and 25 inhibitory independent inputs (pool of 2500 neurons), 
where both types maintain fixed maximal conductances according to Equ.~\ref{eq02}.
We chose a 4:1 ratio of feedforward excitatory to inhibitory connection to mimic 
a balanced state of network activity (\cite{Vre96}).  In sum, each neuron in $R$ received synaptic current from four sources at the same time, as following:

\be
	I_{i}^{recurr}\ =\ I_{i}^{bg}\ +\ \sum_{j=0,j{\neq}i}^{N_{R}=\rho N} I_{i}^{recurr}\ +\ 
                           \sum_{j=0}^{N_{T/C}=N} I_{j}^{train/recon}\ +\ 
                           \sum_{j=0}^{N_{I}} I_{j}^{inh}
	\label{eq03} ~ , 
\ee
where $N=100$ and $N_{I}= 25$.  In the experimental scenario, either $T$ or $C$ was presented as an input (not 
both at the same time). The local recurrent excitatory connections in $R$, on one hand, and feedforward 
connections between $C$ and $R$, on the other hand, were subject to STDP rule during learning and recall stage, 
respectively.  The connection probability ${\rho}$ and size of $R$ were varied to estimate their 
impact on the pattern storage and reconstruction. 

The feedforward layers, $T$ and $C$, were simulated to reproduce an ongoing activity at $50\ Hz$, 
which was spatiotemporally correlated with the background activity in $R$.
The single input of excitatory background activity $I_i^{bg}$ (denoted as $E$) 
fired at a peak rate of $1000\ Hz$ to represent the large number of afferents on $R$. 
 (see Sec.~\ref{ongoing} for an explanation).
The 25 inhibitory inputs $I_j^{inh}$ were simulated as Poisson process at a fixed 10 $Hz$, whose 
role was to prevent $R$ from over-firing. 

During the training phase, feedforward connection strengths {\wffji} between $T$ and $R$ were 
set to the following template pattern: 

\bea
	\omega^{ff}_{ji}  &=& \left\{ \begin{array}{cc}
	\overline \omega^{train} ~, &  {\mathrm if} {\mid{j-i}\mid}\bmod(100-d) \leq d \\
	0 nS ~,&  {\mathrm if} {\mid{j-i}\mid}\bmod(100-d) > d \\
	\end{array} \right. \label{eq04} ~. %\nonumber \\
\eea

Viewed as a square matrix of $N \times N$ binary connections (when $\rho = 1\%$), this template 
assumed the shape of a main diagonal of width $2d$, with periodic boundary conditions (creating 
two additional triangular corners of width $d$; see Fig.~\ref{recon_c08}). Typically, we set 
$d = 20$ for  $N = 100$. The upper bound value of synaptic weight was set to {\wt} = $5\ nS$.
Network development proceeded in two phases (Fig.~\ref{scheme}):

\begin{itemize}
\item During the training phase, the synaptic pattern of $T \xrightarrow[forward]{feed}R$ connections was clamped to {\wffji} according to Equ.~\ref{eq04}, while the recurrent $R \to R$ 
synapses were allowed to develop according to the STDP rule. This represented a learning process 
in which sensory stimulus was applied to $R$ to form a memory.

\item During the reconstruction phase, another layer of $C \xrightarrow[forward]{feed}R$ connections was put in place and subjected to the STDP learning rule, while $R$'s synaptic weights were maintained at 
their mature value created during the training phase. The results were the recovery of the same $T \xrightarrow[forward]{feed}R$ pattern used during training, {\ie} the $C \xrightarrow[forward]{feed} R$ connections converged to {\wffji}.
\end{itemize}

\subsection{Spatiotemporally correlated ongoing activity} \label{ongoing}

Ongoing activity is simply generated by correlating the background input $E$ of recurrent 
layer $R$ with either the training layer $T$ or the reconstruction layer $C$. The neurons of 
$T$ and $C$ generated spikes following an inhomogeneous Poisson random process, in which 
the spatio-periodic firing rate of neuron $i$ depends on the stimulus location $s$ (\cite{Son01}): 

\be
	r_i\ =\  r_{max}\ \times\ (e^{-(s-i)/2{\sigma}^2}\ +\ e^{-(s-i+100)/2{\sigma}^2}\
                 +\ e^{-(s-i-100)/2{\sigma}^2}) %\nonumber
	\label{eq05}  ~ ,
\ee
where $r_{max} = 50\ Hz$ and ${\sigma} = 10$. 
According to Equ.~\ref{eq05}, neurons exhibited maximal 
firing rate at location $s$ (the last two terms are added to implement periodic boundaries 
in 1D space). To generate time-varying firing rates, we update $r_i$ for each time 
interval, which was drawn from exponential distribution with mean {\tauc} = $20\ ms$, 
unless specified. Here, we define the decay time constant of stimuli {\tauc} 
as the temporal correlation constant. The background input of $R$ (denoted $E$ above) 
generated spikes through the same mechanism (${\sigma} = 10$, {\tauc} = $20\ ms$), 
except with $ r_{max} = 1000\ Hz$ to mimic massive afferent activity. 

Experiments show that neurons with higher instantaneous firing rate, are more likely to 
correlate with the spontaneous ongoing activity (\cite{Tso99}). In order to train and reconstruct 
the synaptic pattern at $R$ layer, the spatial correlation was introduced between the 
background input of $R$ layer and $T/C$ layers according to the following equation: 

\be
	s_{rc/ff}\ =\  s\ +\ (N + \sqrt{c} \times\ (1\ -\ N))\ \times\ g %\nonumber
	\label{eq06}  ~ ,
\ee
where $c\ (0 \le c \le 1) $ is the spatial correlation coefficient, 
$g \sim U[0,1]$ is standard uniform random number (\cite{Rud01,Zou07}). For every 
interval, a random location $s$ was first generated from a uniform distributions between 0 
and $N-$1 ($N=100$).  At the start of each interval, the stimulus location in feedforward 
$T/C$ layer, $s_{ff}$ and in recurrent $R$ layers, $s_{rc}$ were correlated at $c$ through 
Equ.~\ref{eq06}, then the firing rates of neurons in each layer were generated according 
Equ.~\ref{eq05} and held constantly until next interval.  In this way, the neurons which 
had the maximal instantaneous firing rates in both $R$ layer and $T/C$ layer, shared the 
spatiotemporally correlated ongoing activity (\cite{Tso99}). 

\subsection{STDP models}
The STDP mechanism we used in our simulation is similar to~\cite{Son00}.
\be
    {d\omega \over dt} \ = \ \Big [ F_{LTP}(t) \ - \ F_{LTD}(t) \ \Big]^{\omega_{LTP}}_{\omega_{LTD}}
    \label{eq07}  ~ ,
\ee
where synaptic weight is bounded between {\wltp} and {\wltd}. Therefore, the modulation 
of $\omega$ has ``hard bounds''.  In contrast, we also consider a model with ``soft 
bounds'', where the changes of $\omega$ depend on the value of the weight (\cite{Zou07}): 

\be
    {d\omega \over dt} \ = \ F_{LTP}(t)(\omega_{LTP}\ -\ \omega) \ - 
                           \ F_{LTD}(t)(\omega \ -\ \omega_{LTD}) \ \label{eq08}  ~ .
\ee
By default, {\wltp} and {\wltd} are set to {\wbar} and 0, respectively in 
both models.
 
We use functions {\fltp} and {\fltd} to describe the coincidence between the pre- and 
post-synaptic spikes for both hard and soft bounds models:
\bea
   F_{LTP} & = & \sum_{t_j<t_i\leq t} \ P(t-t_j) \ \delta(t-t_i) \nonumber \\
   F_{LTD} & = & \sum_{t\geq t_j>t_i} \ Q(t-t_i) \ \delta(t-t_j) 
   \label{eq09}  ~ ,
\eea
where $t_j$ and $t_i$ are the timing of pre- and post-synaptic spikes,
respectively, $\delta(t)$ is Dirac delta function representing a spike at time $t$, 
and all possible pairs of pre- and post-synaptic spikes are 
integrated by the summation.  $P(t) = A_p  \exp(-t/\tau_p)$ and  $Q(t) = A_q \exp[-|t|/\tau_q]$ 
define the time window of interaction
between spikes ({\taup} $=$ {\tauq} $= 20\ ms$).  For recurrent synapses, 
learning take {\ap} $=$ 0.005, and {\aq} $=$ 0.00525, which gives  
{\aq\tauq}$/${\ap \taup} $=$ 1.05 in the range of model from~\cite{Son00}. 
Between feedforward and recurrent synapses, {\ap} $=$ 0.005, and {\aq} $=$ 0.0058, which makes 
{\aq\tauq}$/${\ap\taup} $=$ 1.16. In addition, longer time window in either LTP or LTD are tested 
together with longer correlation time constant for the model of Song {\et} (2000). 

\subsection{Similarity of synaptic patterns}
The output synaptic patterns of $C$ layer were quantified by determining how similar the 
patterns were close to the initial synaptic patterns of $T$, which are constrained by 
Equ.~\ref{eq04}. Although synaptic pattern of $R$ developed in a self-organized way, 
interaction between decay time constant {\tauc} and STDP function
swaps the off-diagonal and diagonal connection in recurrent network $R$ (see Sec.~\ref{decay}). 
This dichotomy suggests that we could simply define the predicted recurrent pattern as the 
off-diagonal connection pattern for intial diagonal connection pattern in $T$:

\bea
	\omega^{rc}_{ji} &=& \left\{ \begin{array}{cc}
	0 nS ~, &  {\mathrm if} {\mid{j-i}\mid}\bmod(100-d) \leq d \\
	\overline \omega^{recurr} ~,&  {\mathrm if} {\mid{j-i}\mid}\bmod(100-d) > d \\
	\end{array} \right. \label{eq10} ~ . %\nonumber
\eea

To measure the magnitude of varying 
quantity, we take the {\it Root Mean Square} ($r.m.s$) between the observed synaptic weight 
and the estimated ones as following:

\be
   r.m.s_{error}\ =\ \sqrt{\frac{1}{M}\ \sum_{ji}(\omega^{observe}_{ji} / {\overline \omega}_{ji}\ -\ 
                     \omega^{predict}_{ji} / {\overline \omega}_{ji})^2} %\nonumber \\
   \label{eq11}  ~ ,
\ee
where {\wobsrv} represents the observed synaptic weight, {\wpredt} represents the initial 
synaptic weight being trained, and $M$ represents the total number of synapses.  {\wbarji} is the maximal 
synaptic weight, which is either {\wc} $= 5\ nS$ for feedforward synapses or {\wrr} $= 1\ ns$ 
for recurrent ones.  Since we take the percentage of the synaptic weights in equation above, 
the final {\rms} will vary between 0 and 1. 

All simulations and analyses were performed by NeoCortical Simulator (NCS) implemented with 
LINUX MPI on a Beowulf cluster (\cite{Bre07,Hog13}). NCS input code were generated by a Python tool 
kit BrainLab (\cite{Drw08}). 

%\end{methods}
% -------------------------------------------------------------------------- %
%                                                                    Results %
% -------------------------------------------------------------------------- %
\section{Results}

\subsection{Synaptic pattern training}
% describe the training
At training stage, we consider the synaptic pattern training in the recurrent 
network by taking the feedforward pattern of $d = 20$, with the maximal strengths 
at the periodic diagonal band.  Both feedforward and recurrent layers were 
spatiotemporally correlated at $c = 0.8$ and{\tauc} $= 20\ ms$. 
After 20k seconds with STDP applied at recurrent synapses, the initial 
uniformly distributed synaptic weights developed the maximal strengths at the off-diagonal 
band, complementary to the feedforward pattern as shown in  Fig.~\ref{train_c08}A (also see the Supplementary Video 1). 
The distribution of synaptic weights is bimodal in Fig.~\ref{train_c08}B, with most of the synaptic weights 
depressed (\cite{Son00}).  The average instantaneous firing rates of the recurrent layer were 
slightly decreased from $12.1\ Hz$ to $10.7\ Hz$ during the first and last 100 seconds, respectively 
during the training process, while the peak ISI stays unchanged (Fig.~\ref{train_c08}C). 

% talk about why off-diagonal
The mechanism of off-diagonal band formed at the recurrent layer could be explained 
with respect to the distribution of spike pair time interval $\Delta{t}$ in 
Fig.~\ref{corr_rms}A.  Each training neuron synapses 40 recurrent neurons according to 
Equ.~\ref{eq04}, and training layer generates periodic peak firing neurons for index $s\pm10$ 
under the correlation of 0.8, whenever a location $s$ is chosen.  Therefore, the recurrent 
neurons most likely will fire with correlation at least to their $\pm10$ recurrent local neighbor 
neurons, when the spiking training neuron was filtered by the fixed feedforward synaptic pattern. 
As a good indication of correlation, the time interval histogram of spike pairs in Fig.~\ref{corr_rms}A 
shows that locally connected synapses on the diagonal band were slightly shifted with $-5\ ms < \Delta{t} < 0\ ms$, 
which gave overall depression when convolving with the negative-favored 
STDP learning; On the contrary, the post-pre synaptic time interval in Fig.~\ref{corr_rms}A shows 
long-range connected synapses at off-diagonal band were slightly more correlated at 
$0\ ms < \Delta{t}< 5\ ms$. Because all connections have 1 $ms$ synaptic delay, 
there is a $1\ ms$ delay for the positive peak of spike timing correlogram in Fig.~\ref{corr_rms}A, dash line. 
The delay is more noticeable in Fig.~\ref{corr_rms}B, solid line. This implies that the coincidently 
spiking neurons will be depressed.  Although the peak of positive correlation is small in 
Fig.~\ref{corr_rms}A, dash line, the potentiation is sufficient enough to overcome the 
negative effect in the STDP learning rule and decorrelation due to synaptic delay, by carefully choosing
the ratio of {\aq\tauq}$/${\ap\taup} as 1.05.

% talk about why off-diagonal
To simplify the discussion, we identify a specific stimulus location as index 
40 for a given time, the training layer elicited maximal firing at recurrent neuron 
index $30\sim50$ through the fixed synaptic pattern.  These neurons will first fire spikes 
with correlation as in Fig.\ref{corr_rms}A (solid line) through local synapses, {\it i.e.}, the diagonal 
region of index $30\sim50$.  The spikes will then propagate to the rest of neuron through 
long-range synapses, {\it i.e.}, the off-diagonal region of $0\sim30$ and $50\sim100$, which will give 
more positive $\Delta{t}$ count in Fig.\ref{corr_rms}A (dash line). 

\subsection{Synaptic pattern reconstruction}
% describe reconstruction
At reconstruction stage, we consider the feedforward synaptic pattern between the reconstruction and 
recurrent layer by freezing the synaptic pattern formed at training stage.  Here, we still take the 
same spatiotemporal correlation ($c = 0.8$, {\tauc} $= 20\ ms$) 
as the training stage. The initial uniformly distributed feedforward 
synaptic weights recovered the maximal strengths at the diagonal band with 20k seconds 
plasticity applied as shown in Fig.~\ref{recon_c08}A (also see the Supplementary Video 1).  
The distribution showed strong bimodal distribution as early as 10k seconds in Fig.~\ref{recon_c08}B.  
The average instantaneous firing rate reduced from $7.9\ Hz$ to $5.8\ Hz$ of the first and last 100 seconds, 
respectively during reconstruction phase, with the peak ISI shifting to the left (Fig.~\ref{recon_c08}C). 

% explain how reconstruction can be done
The reason that the reconstruction layer can recover the pattern trained at recurrent layer 
before is due to the spatiotemporally correlated on-going activity.  If both reconstruction 
and recurrent layer are activated simultaneously at the same neuron index, the reconstruction 
layer will have the accurate spike-timing response from recurrent layer so that feedforward synapses 
can be recovered back to the original pattern.  Fig.~\ref{corr_rms}B. illustrates the spike pairs
distribution between reconstruction and recurrent layers.  
The synapses located at diagonal region had a large peak with $0\ ms <\Delta{t} < 5\ ms$ 
(Fig.~\ref{corr_rms}B, solid line), and was potentiated. Because of this large correlation peak, 
we could choose {\aq\tauq}$/${\ap \taup} up to 1.16, but still develop the potentiated diagonal pattern. 
This could be understood as neurons without local connection responded more to 
pre-synaptic neurons. While synapses located along the off-diagonal region had a trough at 
$0\ ms < \Delta{t} < 5\ ms$ (Fig.~\ref{corr_rms}B, dash line), which in turn depressed the synapses.  
This is due to the neurons with long-range connection firing spikes spontaneously and being reluctant 
to correlate with the pre-synaptic neurons.  Although the spike pair counts at $|\Delta{t}| > 25\ ms$ are 
unbalanced (Fig.~\ref{corr_rms}B), it has least impact on the overall effect when convolved with STDP 
modification function, because the largest modification happens at $|\Delta{t}| < 5\ ms$. 
The spatiotemporally correlated ongoing activity gave the reconstruction layer 
an exact reference of where the information is stored, and the asymmetric STDP rule 
shaped the synaptic pattern.  
As a control, we take a random recurrent pattern during reconstruction stage, and 
the resulting feedforward connection doesn't show any specific pattern, but rather synaptic 
competition as in~\cite{Son00} (figures not shown). This suggests that previously learned 
connectivity within $R$ is necessary to reconstruct the feedforward connectivity.

% root mean square calculation
To measure how well the synaptic patterns were trained and reconstructed, we calculated 
the root mean square ({\rms}) between the actual synaptic map and the synaptic template defined 
as Equ.~\ref{eq07}. The synaptic weights were sampled at every second, and fitted to exponential 
function as: $E_{rms} \times (1-e^{-t^2/\tau_{rms}})$, where $E_{rms}$ is the converged root 
mean square, and $\tau_{rms}$ is the time constant for the {\rms} to converge.  The trained synaptic 
map was index-wisely compared with the template map.  The recurrent {\rms} were almost saturated 
at $39.48\%$ after 20k seconds of training, while the reconstructed synaptic pattern had similar 
quality with $E_{rms} = 38.26\%$, although it took longer to saturate in Fig.~\ref{corr_rms}C. 

\subsection{Effect of spatial correlation on pattern formation}
% correlation dependency of training and reconstruction
We next investigated the effect of spatial correlation between feedforward layer 
$T/C$ and recurrent layer $R$ (see Sec.~\ref{network}).  
The {\rms} of the stable synaptic weights is calculated to estimate how 
the synaptic map was trained and reconstructed from the template pattern. The time constant of 
the {\rms}, {\trms} is also estimated from the exponential function to predict how fast the synaptic 
map reaches its stable state. 
Fig.~\ref{correl_ff_rc}A 
indicates that {\rms} decreases when the spatial correlation increases for both 
training and reconstruction. 
However, spatial correlation is NOT necessary to form the recurrent pattern during training: the 
quality of recurrent synaptic map does not change from  $c = 0.0$ to $c = 0.6$ 
(Fig.~\ref{correl_ff_rc}A, circle line and Fig.~\ref{correl_ff_rc}C, upper row). 
This suggests that the emerging weight 
structure in recurrent network does not depend on spatial correlation, and is ensured by 
the temporal correlation {\tauc} $= 20\ ms$ (see Sec.~\ref{decay} for effect of temporal correlation).
In contrast with the training, spatial correlation significantly 
changed the {\rms} during reconstruction (Fig.~\ref{correl_ff_rc}A, triangle line and 
Fig.~\ref{correl_ff_rc}C, bottom row). The pattern is not formed until $c = 0.6$. 
This suggests that STDP operates mainly via rate effects. 

The time constant of {\rms} during training {\trmst} ranges from 144 to 
220 seconds, and correlation monotonically decrease the time constant 
(Fig.~\ref{correl_ff_rc}B, circle line). By examining the scales of time units, 
we notice that the time constant for reconstruction is longer than the one of training by an 
order of magnitude. This result could be explained by the fact that fast learning usually 
is followed by more gradual improvements in memory (\cite{Aln07}).
However, the effect of the correlation on time constant of {\rms} during reconstruction {\trmsc} is 
rather complex: The synaptic map converged to stable value within 100 seconds without correlation 
introduced. {\trmsc} jumped up to about 7700 seconds 
when correlation was first introduced, then start to decrease when reconstructed synaptic pattern 
emerged at $c = 0.6$ (Fig.~\ref{correl_ff_rc}C, bottom row). 

\subsection{Influence of decay time constant of stimuli and STDP function} \label{decay}
% tau_corr = 100 ms
During both training and reconstruction of synaptic pattern, both recurrent and feedforward networks 
receive input stimuli generated by a short period of burst presented at random locations. The decay 
time constant of the stimuli, {\ie} the time lasting before an input burst moves to another 
random location, was chosen from an exponential distribution with a mean of 20 ms in previous simulation.
To investigate the influence of longer time constant of stimuli on the synaptic pattern formation, we 
increase the time constant of stimuli, {\tauc}, from 20 to 100 ms (Fig.~\ref{tau_ff_rc}A). The increase of 
{\tauc} strongly degrade the quality of the synaptic pattern formed in feedforward network 
(Fig.~\ref{tau_ff_rc}A, bottom middle), because of the disappearing of the potentiation peak in 
correlation (Fig.~\ref{tau_ff_rc}B, compare the dot dash line to solid line). 
However, longer {\tauc} preserved the pattern formed in recurrent network (Fig.~\ref{tau_ff_rc}A, 
top middle), although {\rms} has 3.84\% increase and it takes about 10 times longer to stabilize 
the pattern compared to {\tauc} $= 20\ ms$. 

% tau_corr = 100 ms, tau_p = 100 ms and tau_corr = 100 ms, and tau_q = 100 ms
To further test the influence of time constant of the STDP function, we matched {\taup} and {\tauq} with 
{\tauc} during both training and reconstruction processes, respectively.  We found that the formation of the 
pattern were improved, if {\taup} matches {\tauc} during training (Fig.~\ref{tau_ff_rc}A, top left) and 
{\tauq} matches {\tauc} during reconstruction (Fig.~\ref{tau_ff_rc}A bottom right). The potentiation 
peak in the correlation (Fig.~\ref{tau_ff_rc}B) is responsible for the improvement of quality of 
the pattern during reconstruction (figure not shown during training). 
Interestingly, the synaptic patterns were swapped, if {\tauq} matches {\tauc} during training 
(Fig.~\ref{tau_ff_rc}, top right) and {\taup} matches {\tauc} during reconstruction 
(Fig.~\ref{tau_ff_rc}, bottom left). 
Both diagonal and off-diagonal connection patterns in recurrent network during training 
(Fig.~\ref{tau_ff_rc}A, upper row) shows that strong synapses are bidirectional.
The relatively larger peak at positive intervals in 
Fig.~\ref{tau_ff_rc}C contributes to the dominance of the potentiation at the swapped region 
of the pattern during reconstruction (figure not shown during training). 


\section{Discussion}\label{discussion}
% summary of the work and compare to Davision's
We presented a novel model wherein simple synaptic patterns are stored into and retrieved 
from a recurrent network through STDP under conditions of correlated spatiotemporal 
activity. As signals, we used very general Poisson spike trains with a Gaussian 
distribution of firing rates (\cite{Son00}) , with correlation centered at the peak 
firing frequency.  Connectivity patterns in a recurrent cortical network have been shown to be 
nonrandom (\cite{Son05}). Here we limit the pattern to simple diagonal connectivity, 
{\ie} short-range synapses, suggest a role for STDP to selectively potentiate synaptic 
patterns (\cite{Zou07}), and promote resonant learning in recurrent 
networks (\cite{Dau02}). 

% biological plausiability. 
Our model could be applied generically to spiking networks that must transfer, store, then 
re-transfer simple synaptic patterns into another network.  
The paradigm might serve the underlying biological mechanism for the reward based 
system of~\cite{Kri13} or virtual neurorobotics (VNR) environment of~\cite{Gdm07}.
This is relevant in context of brain areas that have been shown to demonstrate 
both STDP and correlated spatiotemporal activity.  
The model provides a microcircuit-level explanation of how correlated activity temporally 
conjoins neurons into assemblies and selectively facilitates synaptic plasticity to support 
information storage and retrieval.  For example, this been observed in hippocampus, especially 
CA3 (\cite{Oke93, Trb89}), which is known for its prominent recurrent synapses.  Functionally, 
synaptic modification in hippocampus during task performance (such as navigation) is compatible 
with the dynamical organization and computational concept of our model.  Likewise, correlated neuronal 
activity in neocortex has been observed while learning sensory input (\cite{Lau94, Han01, Tat05, 
Pil08}), retrieving sensorimotor associations (\cite{Abl93, Vad95}), and making decisions (\cite{Shd01}). 
Cerebellar studies also show correlated spatiotemporal patterns (\cite{Mid08, Pil08}).  A series 
of combined electroencephalogram (EEG) and intracellular {\it {in vivo}} measurements suggested 
that the required spatiotemporal correlation can arise from synchronous oscillation generated 
by the correlated activity of a population of neurons, during both awake and sleep states (\cite{Con97, 
Des99, Buz04}).

% differentiate the time correlation with previous work
In this work, we define the temporal correlation constant, {\tauc}, 
as the mean of the update time 
intervals for inhomogeneous Poisson process, and study the relationship between {\tauc} and STDP 
learning time constant (either {\taup} or {\tauq}).  
This temporal correlation mechanism represents global discontinuity
for saccades in vision and sniff in olfaction (\cite{Uch06}), 
or brain oscillation (\cite{Con97,Des99, Buz04}). 
This contrasts with previous works that temporal sequences were presented continuously where 
afferents fire continuously with a constant population rate (\cite{Mas08,Now03}).

% time constant issue
How sensitive are the parameters of STDP to achieve learning and recall in our model? We first 
explored the effect of changing the background temporal correlation time constant, {\tauc}, 
on pattern formation during each of the phases.  Surprisingly, storage of the feedforward 
to recurrent synaptic pattern was minimally affected by increasing {\tauc} alone, in the physiological 
range of 20 to 100 $ms$ (Fig.~\ref{tau_ff_rc}A, top left).  
Increasing {\tauc} in concert with an increase of the positive learning time 
constant, {\taup} as high as 100 $ms$ enhances the efficiency of initial pattern storage of into the 
recurrent network (Fig.~\ref{tau_ff_rc}A, top middle). In either case, the synaptic patterns
have formed a long-range bidirectional connection.  Although such a large {\taup} has 
not yet been reported, it is believed that such 
a widened LTP window can result from NMDA channel activation or calcium diffusion (\cite{Bdl05}).  
On the other hand, increasing {\tauc} alone will impair the recalled pattern (as reported also by~\cite{Son01}); however, we found that this impairment is alleviated if the negative 
learning time constant, {\tauq}, is concurrently increased to 100 $ms$ (which also results in an 
inversion of the stored pattern of synaptic strengths in both recurrent and feedforward networks). 
Prolongation of {\tauq} on this scale has been reported in rats (\cite{Fdm00}) and electric 
fishes (\cite{Bel97}).
In this case, a short-range bidirectional connection has been formed in the recurrent network.
Nevertheless, Both long- and short-range bidirectional connections are potentiated (\cite{Son05}).

% soft-bound model and size and connectivity
We also evaluated the impact of imposing ``soft-bound'' STDP rules (\cite{Ksl00, Ros00, Bdl05, 
Mor07, Zou07}). 
We found that ``soft-bound'' STDP results in the formation of synaptic patterns at about the same 
rate as does the ``hard-bound'' model, but it produces a low contrast pattern, attributable to the unimodal 
distribution of the synaptic weights. Thus our spatiotemporal correlation--based computing paradigm of 
synaptic pattern storage and retrieval is independent from whether ``soft-'' or ``hard-bound'' 
rule is applied.  
The effects of the recurrent network size and connectivity are investigated by performing the pattern 
training and reconstruction as well.  Interestingly, all simple recurrent patterns, although at different levels of connectivity and size,  contain essential information to reconstruct the 
feedforward patterns initially trained.

% talk about the limit of the connectivity and size
The accuracy of pattern storage of our model was also evaluated using root mean square ($r.m.s$) error 
while co-varying the connectivity and size of the recurrent network.  In general, a sparsely connected 
but larger size recurrent network (for example, 20\% with 100 neurons) will have better performance 
on pattern storage and retrieval. We did not study recurrent networks larger than 100 neurons, due to the 
computational overhead involved.  Further work should also include the evaluation of clustering of  
complex synaptic patterns under different conditions of brain activity. 


\section*{Disclosure/Conflict-of-Interest Statement}
%Frontiers follows the recommendations by the International Committee of Medical Journal Editors (http://www.icmje.org/ethical_4conflicts.html) which require that all financial, commercial or other relationships that might be perceived by the academic community as representing a potential conflict of interest must be disclosed. If no such relationship exists, authors will be asked to declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. When disclosing the potential conflict of interest, the authors need to address the following points:
%•	Did you or your institution at any time receive payment or services from a third party for any aspect of the submitted work?
%•	Please declare financial relationships with entities that could be perceived to influence, or that give the appearance of potentially influencing, what you wrote in the submitted work.
%•	Please declare patents and copyrights, whether pending, issued, licensed and/or receiving royalties relevant to the work.
%•	Please state other relationships or activities that readers could perceive to have influenced, or that give the appearance of potentially influencing, what you wrote in the submitted work.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

%\section*{Author Contributions}
%When determining authorship the following criteria should be observed:
%•	Substantial contributions to the conception or design of the work; or the acquisition, analysis, or interpretation of data for the work; AND
%•	Drafting the work or revising it critically for important intellectual content; AND
%•	Final approval of the version to be published ; AND
%•	Agreement to be accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved.
%Contributors who meet fewer than all 4 of the above criteria for authorship should not be listed as authors, but they should be acknowledged. (http://www.icmje.org/roles_a.html)

%The statement about the authors and contributors can be up to several sentences long, describing the tasks of individual authors referred to by their initials and should be included at the end of the manuscript before the References section.


\section*{Acknowledgments}
In memory of Dr. Philip H. Goodman. 
%\textit{Funding\textcolon} Text Text Text Text Text Text  Text Text.

\section*{Supplemental Material}
The Supplementary Information, Supplementary Videos 1 and 2 for this article can be found at \href{http://www.cse.unr.edu/brain/}{http://www.cse.unr.edu/brain/}. 

\bibliographystyle{frontiersinSCNS_ENG_HUMS} % for Science, Engineering and Humanities and Social Sciences articles, for Humanities and Social Sciences articles please include page numbers in the in-text citations
%\bibliographystyle{frontiersinHLTH&FPHY} % for Health and Physics articles
\bibliography{references}

%%% Upload the *bib file along with the *tex file and PDF on submission if the bibliography is not in the main *tex file

\section*{Figures}

% scheme figure @ Intro
\begin{figure}
\begin{center}
\includegraphics[viewport=0 0 695 194, width=6in]{figs/scheme.eps}
\end{center}
\caption{
{\bf Training and reconstruction of synaptic patterns.}  (A) During training phase, 
The $T \to R$ synaptic pattern was fixed (diagonal matrix). The recurrent synaptic weights 
in $R$ were subject to STDP rules and, starting from a uniform initial state, developed into 
an off-diagonal synaptic pattern. (B) During reconstruction phase, the recurrent synaptic 
pattern inside R was fixed, while $C \to R$ synapses were subject to STDP rules. Starting 
from a uniform distribution, the $C \to R$ synaptic weights reconstructed the diagonal 
pattern that was used during training in (A). Both $T,R$ and $C,R$ received 
spatiotemporally correlated background stimuli $E$ in the 2 phases.}
\label{scheme}
\end{figure}

% training @ Results
\begin{figure}
\begin{center}
 \includegraphics[viewport=0 0 596 585, width=4in]{figs/train_rc_corr08_dev.eps}
\end{center}
\caption{
{\bf Training diagonal pattern at recurrent layer.} 
(A) Snapshot of recurrent synaptic map was plotted at different time by averaging 
each synaptic strength over 10 seconds. Each pixel represents the synaptic strength 
from pre-synaptic neuron $j$ to post-synaptic neuron $i$, which is normalized between 
0 and 1 (see first snapshot for weight scale).  
All the synaptic weights are initialized uniformly. The synapses were 
selectively potentiated at the off-diagonal band, and depressed for the rest area 
during the training stage. The black pixel indicates self-connection which was excluded (see Supplementary Video 1).
(B) Histogram of the synaptic weights plotted every 10 seconds 
over time. The bimodal distribution was formed at the end of training. 
(C) ISI distribution pooled from all the neurons of recurrent layer of 100 seconds duration 
shows no difference at start and end of the training, where both peak at 6$\sim$7 $ms$ by noting 
the time is in log scale. 
}
\label{train_c08}
\end{figure}


% reconstruction @ Results
\begin{figure}
\begin{center}
 \includegraphics[viewport=0 0 594 585, width=4in]{figs/recon_ff_corr08_dev.eps}
\end{center}
\caption{
{\bf Rebuilding diagonal pattern at reconstruction layer.}  
(A) Snapshot of feedforward synaptic map was plotted at different time by averaging 
each synaptic strength over 10 seconds.  All the synaptic weights are initialized 
uniformly.  The synapse map were exactly rebuilt as the periodic diagonal band which 
is trained at the recurrent layer before (see Supplementary Video 2). 
(B) Histogram of the synaptic weights plotted every 10 seconds over time. The strong 
bimodal distribution was formed at the end of training. 
(C) ISI distribution pooled from all the neurons of recurrent layer of 100 seconds duration 
shows the peak ISI was slightly shifted to the left after the reconstruction.
}
\label{recon_c08}
\end{figure}

% reconstruuction correlation @ Results
\begin{figure}
\begin{center}
 \includegraphics[viewport=0 0 356 741, width=2in]{figs/corr_rms.eps}
\end{center}
\caption{
{\bf Spike pairs histogram and variances of synaptic patterns.}  
All possible spike pairs $\Delta{t} = t_{post} - t_{pre}$ are binned for synapses both 
at diagonal band (solid line), and at off-diagonal band (dash line) with 1 $ms$ time step 
after the synapses converged.  The distribution were normalized to the total spike 
pair counts. Spike pair counts increase at small $\Delta{t}$ for diagonal synapses, 
and decrease at small $\Delta{t}$ for off-diagonal synapses.  However,  there are
 more counts for positive $|\Delta{t}| < 5\ ms$ at off-diagonal recurrent synapses 
(A), dot line and at diagonal feedforward synapses (B), solid line. In contrast, 
there are more counts for negative $|\Delta{t}| < 5\ ms$ at diagonal recurrent synapses 
(A), solid line and at off-diagonal feedforward synapses (B), dash line. 
The convolution of histogram with STDP function was evaluated 
to verify the weight drift (see text for details).
(C) Root 
mean square ({\rms}) of the varying magnitude between the estimated synaptic weights and 
the observed ones were measured every second (only plot with every 1000 seconds).  The {\rms} 
were fitted to exponential functions. The recurrent synaptic pattern (filled circle) 
has $E_{rms} = 39.48\%$, $\tau_{rms} = 144.4\ s$; the reconstruction synaptic pattern 
(filled triangle) has $E_{rms} = 38.26\%$, $\tau_{rms} = 900.6\ s$. 
}
\label{corr_rms}
\end{figure}

% vary the correlation @ Results
\begin{figure}
\begin{center}
 \includegraphics[viewport=0 0 572 551, width=4in]{figs/correl_ff_rc.eps}
\end{center}
\caption{
{\bf Effect of correlation on the {\erms} and {\trms}.}
(A) Effect of correlation on the converged {\rms} of training template pattern (filled circle) 
and reconstruction of the trained pattern (filled triangle). (B) Corresponding time constants 
were estimated by fitting {\rms} to exponential function, and then plotted as a function of 
correlation.  The time constant for reconstruction is longer than the one of training on an order 
of magnitude by noting the different scales. (C) Synaptic map of both training (top row) and 
reconstruction (bottom row) were plotted at 20k seconds for correlation being 0.0, 0.6 and 1.0. 
}
\label{correl_ff_rc}
\end{figure}


% tau decay @ Results
\begin{figure}
\begin{center}
 \includegraphics[viewport=0 0 603 595, width=4in]{figs/tau_ff_rc.eps}
\end{center}
\caption{
{\bf Effect of decay time constant of stimuli and STDP window function.}  
(A) Both recurrent (top) and feedforward network (bottom) were undergone stimuli with correlation 
time, $\tau_{corr}$, at 100 ms.  Feedforward network had a degradation of the reconstructed 
synaptic pattern by noting the scale of the conductance (bottom middle), while recurrent 
network was not affected (top middle).  The synaptic pattern can be preserved by matching 
$\tau_{corr}$ with $\tau_{p}$ during training (top left) and $\tau_{q}$ during reconstruction 
(bottom right), respectively.  However, the synaptic patterns were swapped if 
$\tau_{q} = \tau_{corr}$ during recurrent training(top right), and $\tau_{p} = \tau_{corr}$ 
during feedforward reconstruction(bottom left). The correlation of the stimuli was 0.8 for all 
the simulations. Histogram of time interval, $t_{post} - t_{pre}$, for $\tau_{corr} = \tau_{p} = 100\ ms$ (dash), 
$\tau_{corr} =  100\ ms$ (solid) and $\tau_{corr} = \tau_{q} = 100\ ms$ (dot dash) at diagonal region (B), 
and off-diagonal region (C) was sampled during reconstruction shown at A bottom.  There are peaks appearing 
at positive time interval for each region, which are responsible for the synaptic pattern 
formed with potentiated synaptic weights. 
}
\label{tau_ff_rc}
\end{figure}



%%% If you don't add the figures in the LaTeX files, please upload them when submitting the article.

%%% Frontiers will add the figures at the end of the provisional pdf automatically %%%

%%% The use of LaTeX coding to draw Diagrams/Figures/Structures should be avoided. They should be external callouts including graphics.

\end{document}
