Optimal Design for Hetero-Associative Memory: Hippocampal CA1 Phase Response Curve and Spike-Timing-Dependent Plasticity

Recently reported experimental findings suggest that the hippocampal CA1 network stores spatio-temporal spike patterns and retrieves temporally reversed and spread-out patterns. In this paper, we explore the idea that the properties of the neural interactions and the synaptic plasticity rule in the CA1 network enable it to function as a hetero-associative memory recalling such reversed and spread-out spike patterns. In line with Lengyel’s speculation (Lengyel et al., 2005), we firstly derive optimally designed spike-timing-dependent plasticity (STDP) rules that are matched to neural interactions formalized in terms of phase response curves (PRCs) for performing the hetero-associative memory function. By maximizing object functions formulated in terms of mutual information for evaluating memory retrieval performance, we search for STDP window functions that are optimal for retrieval of normal and doubly spread-out patterns under the constraint that the PRCs are those of CA1 pyramidal neurons. The system, which can retrieve normal and doubly spread-out patterns, can also retrieve reversed patterns with the same quality. Finally, we demonstrate that purposely designed STDP window functions qualitatively conform to typical ones found in CA1 pyramidal neurons.


Introduction
It has been reported that characteristic ensemble spiking patterns are consistently repeated in the hippocampal CA1 region during waking and sleep periods [1][2][3][4][5][6][7][8][9]. Louie and Wilson (2001) reported that in rats, spike patterns produced during rapid eye movement (REM) episodes are very similar to those observed while the animals are running [8]. There are cases in which the timescale of these reactivation patterns during REM episodes is on average twice as long as that of the running periods. Foster and Wilson (2006) reported that the spike patterns observed during running periods are reproduced in a temporally reversed order during rest periods [9]. These experiments suggest that the CA1 network stores spatio-temporal spike patterns and retrieves reversed and spread-out patterns.
These experimental results raise a big question as to whether the hippocampal network including the CA1 and CA3 regions has an optimal structure for storing and retrieving such spike patterns. Lengyel et al. (2005) [10] developed a normative theory for autoassociative memory networks that specifies optimal pairs of the synaptic plasticity rule for embedding memories and the form of neural interactions for auto-associative memory retrieval. Under the speculation that a phase response curve (PRC) is appropriate to formulate the neural interactions if memories are embedded by spike-timing-dependent plasticity (STDP), they derived pairs of STDP window functions and PRCs optimally functioning as an auto-associative memory. They showed that the features of the PRCs of hippocampal CA3 pyramidal neurons qualitatively conform to ones theoretically derived from typical STDP window functions. However, they asked the question only as it relates to restoring phase patterns to the original stored state through mutual recurrent interactions, not retrieval of reversed and spread-out patterns. Moreover, although the possible existence of STDP at recurrent synapses between CA3 pyramidal neurons has been suggested [11,12], as far as we know, there are no reports on capturing the entire shape of the STDP window function.
Here, we focus on the CA1 network in which STDP has been reported [13][14][15][16][17]. We explore the idea that the properties of the neural interactions and the synaptic plasticity rule support the function of hetero-associative memory in which spike patterns are embedded in synapses and reversed and spread-out patterns are retrieved. In line with Lengyel's speculation, we search for optimal pairs of STDP window functions and PRCs. Whereas Lengyel et al. used a top-down approach, treating the auto-associative memory retrieval as optimal probabilistic inference and inferring the retrieval dynamics that are normatively matched to the typical STDP window functions, we take a synthetic approach of optimal design for hetero-associative memory under the physical limitations of the neural implementation. Figure 1 illustrates our approach, consisting of bottom-up and top-down steps. In the bottom-up steps, under the assumption of regular firing and weak coupling, we firstly formulate a hetero-associative memory network recalling not only the normal spike patterns, but also the reversed and doubly spread-out patterns as a phase oscillator model consisting of an STDP window function and a PRC. This network model associates pre-and postsynaptic phase patterns. For example, when presented with a stored presynaptic phase pattern that is temporally reversed or spread out, the postsynaptic neurons can recall the associated phase pattern that is temporally reversed or spread out (for more detail, see Table 1). Secondly, we analytically derive the mutual information between a stored phase pattern and a network output, and use it to evaluate memory retrieval performance. In the top-down steps, by maximizing the objective function given by the mutual information, we search for a set of optimal STDP window functions under the constraint of PRCs recorded in vitro from hippocampal CA1 pyramidal neurons [18,19].
The theoretically derived STDP window functions are compared with those reported for CA1 pyramidal neurons. The typical STDP window functions observed in the CA1 region are classified into two types [13,14]: symmetric [15,16] and asymmetric [17] plasticity rules. We show that both of these rules are included in a theoretically derived set of optimal STDP window functions and they allow network models with them to work as an associative memory.

Working Hypothesis
Hippocampal CA1 network works as hetero-associative memory. On the basis of the anatomical structure and physiological properties of the hippocampus, many researchers have long hypothesized that auto-associative memory resides in the CA3 region because of its recurrent connections, and heteroassociative memory resides in the CA1 region because of its feedforward connections [20][21][22][23][24][25][26][27][28].
CA1 pyramidal neurons receive inputs from the entorhinal cortex and the CA3 region, and the temporal correspondence between the activity patterns of CA1 pyramidal neurons and these presynaptic activity patterns may result in hetero-association between them by modification of the synapses onto the CA1 pyramidal neurons [27]. In line with these considerations, we assume that the hippocampal CA1 network works as heteroassociative memory and introduce a feedforward network model. The fundamental requirements for hetero-associative memory are to recall an associated activity pattern of postsynaptic neurons upon presentation of a presynaptic activity pattern.
Lengyel's speculation. The STDP is an associative plasticity that adjusts synaptic efficacy depending on the relative timing of pre-and postsynaptic spikes. In the case of asymmetric STDP window functions, a synapse increases its efficacy if presynaptic spikes repetitively arrive within 5-40 msec before the postsynaptic spikes, whereas the same synapse decreases its efficacy if presynaptic spikes repetitively arrive with a similar time window after the postsynaptic spikes. On the other hand, the PRC reflects the sensitivity of oscillatory postsynaptic spike timing in response to presynaptic spike activation or a current perturbation mimicking presynaptic spike activation. The experimental protocols for measuring PRCs generate presynaptic spikes or inject perturbation currents at various timings relative to the last spike of repetitively firing postsynaptic neuron, and measure the inter-spike interval of the cycle containing the perturbation. The STDP window function and the PRC respectively indicate the effect of the timing of the presynaptic spikes relative to the postsynaptic ones on the synaptic efficacy and the timing of postsynaptic spikes, and they are based on the premise that neurons act as oscillators. In light of this, Lengyel et al. speculated that PRCs are an appropriate way to formulate neural interactions if memories are embedded by STDP [10]. In line with this speculation, we decided to search for optimal pairs of STDP window functions and PRCs.
Phase patterns of presynaptic neurons are associated with those of postsynaptic neurons in the hetero-associative memory. In the storage process, p pairs of pre-and postsynaptic phase patterns, g m and h m (m~1,2, Á Á Á ,p), are embedded by modifying synaptic weights in accordance with an STDP learning rule. In the retrieval process, when presented with a phase pattern of presynaptic neurons which resembles the m-th memory key pattern that is temporally reversed and/or stretched to jaj times its original timescale, y (&ag m ) (a~+1,+2), the postsynaptic neurons recall a phase pattern which resembles the associated memory output pattern that is temporally reversed and/or stretched to jaj times its original timescale, w (&ah m ). doi:10.1371/journal.pone.0077395.t001 Phase reduction of weakly coupled limit-cycle oscillators. The hippocampal CA1 region, as well as other regions involved in memory processing, exhibits stable oscillations of the local field potential (LFP) in several situations including awake and sleep states [29]. In such cases, the temporal order of neuronal spiking relative to the LFP oscillation are preserved [30] and correlated with the animal's location in space [31]. This evidence suggests that memories seem to be encoded in spike times relative to ongoing LFP oscillations. Here, we formulate the storage and retrieval processes as coupled limit-cycle oscillators. As described below, under the assumption of regular firing and weak coupling, we can reduce this oscillator system to a phase equation on the basis of PRC [32,33]. Thus, the use of PRCs is consistent with Lengyel's speculation described above. However, the temporal spike patterns of hippocampal CA1 pyramidal neurons in vivo show irregular and bursting activities that differ from the regular spike activity that we assume here [29]. The phase equation derived here is valid in the limit that each neuron generates a single spike once during each period of the collective oscillations.
Our assumption does not correspond to actual behaviors very well, but it is expected to capture the dominant factor of cooperative behavior in coupled oscillating systems [34]. Another big advantage of this analysis method is that it can be directly applied to real neurons by electrophysiologically measuring the PRCs [33].
Here, we use the PRCs of hippocampal CA1 pyramidal neurons recorded in vitro [18,19], and predict the behavior of a virtual hetero-associative memory network composed of pyramidal neurons.
Retrieval of doubly spread-out patterns under the same collective theta oscillations in the running periods. It has been reported that in rats, awake neural ensemble activities are reproduced during REM episodes associated with increases in LFP theta power, and that the timescale of reactivation patterns during the REM episodes is on average twice as long as that of running periods (see Fig. 5B in [8]). This result suggests that during REM episodes, doubly spread-out patterns are reactivated under the same collective theta oscillations as those during running periods.
Here, under the assumption that oscillatory dynamical properties of each neuron does not change in the CA1 network even though doubly spread-out patterns are reactivated, we use the same PRC when formulating the retrieval processes for normal and doubly spread-out patterns. This assumption simplifies the problem to find pairs of PRCs and STDP window functions optimally recalling normal and spread-out patterns.

Minimum Model Functioning as a Temporal Heteroassociative Memory
In this study, we use the hetero-associative memory model shown in Fig. 2A. The model consists of N presynaptic and M postsynaptic oscillator neurons. For simplicity, these pre-and postsynaptic neurons have the same firing period. A presynaptic neuron j is connected to a postsynaptic neuron i through a synapse with an efficacy (weight) J ij . The theoretical derivations presented below assume all-to-all connectivity.
For the purpose of mathematical tractability and simplicity, we assume that the timescale of synaptic dynamics in the memory storage process is far different from that of network dynamics in the memory retrieval process. Under such an assumption of timescale separation, the storage process and retrieval process can be separated from one another.
Set out below is the outline of the hetero-associative memory functions that we wanted to study (Table 1). In the storage process, p pairs of pre-and postsynaptic phase patterns, g m and h m (m~1,2, Á Á Á ,p), are embedded by modifying synaptic weights in accordance with an STDP learning rule. g m and h m are called the m-th memory key pattern and memory output pattern, respectively. In the retrieval process, after assigning a phase pattern of the presynaptic neurons that resembles the memory key pattern, y (&g m ), the postsynaptic neurons recall a phase pattern that resembles the memory output pattern, w (&h m ). Furthermore, we treat more general cases involving the normal phase pattern retrieval described above. When presented with a phase pattern of presynaptic neurons resembling the memory key pattern that is temporally reversed and/or stretched to jaj times its original timescale, y (&ah m ) (a~+1,+2), the postsynaptic neurons recall a phase pattern that resembles a memory output pattern that is temporally reversed and/or stretched to jaj times its original timescale, w (&ah m ). Here, y and w are called retrieval key pattern and retrieval output pattern, respectively. Note that the synaptic weights are then fixed during the retrieval process. The case a~1 corresponds to normal phase pattern retrieval, the case a~2 corresponds to doubly spread-out pattern retrieval, and the cases a~{1 and {2 correspond to retrievals of reversed patterns. The following subsections describe the storage and retrieval processes in the network.
Synapse dynamics in the storage process. In the storage process, we treat g j as the phase of presynaptic neuron j (~1,2, Á Á Á ,N), and h i as that of postsynaptic neuron i (~1,2, Á Á Á ,M). J ij denotes the synaptic weights between presynaptic neuron j and postsynaptic neuron i. Memory storage occurs as a result of synaptic modification depending on the relative phase of the pre-and postsynaptic neurons. The amount of synaptic modification, DJ ij , is determined according to the following synaptic plasticity rule: where V(h) is the STDP window function. This rule is local in that the change to J ij depends only on the phases of these two neurons and not on those of other neurons. When storing more than one pair (m~1,2, Á Á Á ,p), we also make a simplifying assumption that synaptic plasticity is additive across the memories: where 1 N J 0 (w0) is the initial synaptic weight, which avoids negative values of J ij . The local and additive plasticity rule is similar to the one described in the previous study [10]. The difference from the previous study is the scaling of the synaptic weight with the number of presynaptic neurons, N. This scaling is necessary for derivation of the order parameter m m k,l , which measures the overlap between the m-th memory key pattern g m and the retrieval key pattern y, as we shall discuss later. Since the magnitude of V(h) is arbitrary, there is no loss of generality due to the scaling in Eqs. (1) and (2). Each element g m j of the m-th memory key pattern g m stored in the model is assigned to an independent random number in ½0,2p) with a uniform probability, P(g m j )~1 2p .
By the same token, each element h m i of the m-th memory output pattern h m is assigned to an independent random number in ½0,2p) with a uniform probability, P(h m i )~1 2p . Thus, these memory patterns are not correlated with each other.

Optimal Design for Hetero-Associative Memory
Network dynamics in the retrieval process. We assume that the postsynaptic neural population ( Fig. 2A) consists of neural oscillators which share common features: the M postsynaptic neurons fire rhythmically with a period T~2p=v (v is the angular frequency), and N presynaptic neurons also fire with a period T. Let us consider a situation in which rhythmical firing of a postsynaptic neuron i is perturbed by a total synaptic input current G i (t) from N presynaptic neurons and an additive noise current es i (t): The term (G i (t)zes i (t))U represents the driving current, in which the vector U (~½1,0, Á Á Á ,0 T ) results in perturbing one degree of freedom of the neural oscillator. s i (t) is a onedimensional Langevin force satisfying Ss i (t)T~0, Ss i (t)s j (t 0 )T~2d ij d(t{t 0 ), and is the intensity of the Langevin force. X i is a high-dimensional state vector that represents the activity of the postsynaptic neuron i, namely, the membrane potential, calcium concentration, and conductances for voltagegated ion channels. F(X i ) is a vector field that represents the intrinsic dynamics of neuron i. We assume that the unperturbed neural oscillator dX i =dt~F(X i ) has a stable limit cycle solution: wherew w i is the phase of postsynaptic neuron i, and w i is the initial phase corresponding to the retrieval output pattern. The total postsynaptic current G i (t) is given by where r(x) denotes the waveform of the postsynaptic current, and J ij is the synaptic efficacy determining the magnitude of the synaptic current.ỹ y j is the phase of presynaptic neuron j, and y j is the initial phase corresponding to the retrieval key pattern. When e and G i (t) are sufficiently small, a high-dimensional system (3) can be reduced to a one-dimensional one expressing the  [18,19]. The abscissa represents the phase of a perturbation arrival, and the ordinate represents the phase shift of the postsynaptic spike in response to the perturbation current. (C) Typical STDP window functions observed in hippocampal CA1 pyramidal neurons. In the storage process, synaptic weights fJ ij g are determined in accordance with an STDP learning rule depending on the phase difference between the pre-and postsynaptic spikes. Left: Symmetric plasticity rule [16]. Right: Asymmetric plasticity rule [17]. doi:10.1371/journal.pone.0077395.g002 motion of the phasew w i in the limit cycle orbit: where Z PRC (w) is the PRC reflecting the sensitivity to the perturbation current [32,33]. This is called the Langevin phase equation (LPE) [35,36]. We apply the variable transformationw w i~v tzw i and averaging to Eq. (6). w i represents a slow-moving initial phase driven by a small perturbation (synaptic input) and noise. Accordingly, we can write the slow dynamics of the initial phase of the postsynaptic neuron, w i , as s~e C(w) is called the coupling function, and it is represented as a convolution of r(x) with Z PRC (x). If the time constant of r(x) is much shorter than the period T of the postsynaptic neuron, r(x) can be written as a delta function (r(x)~d(x)). Thus, the coupling function can be written using the PRC: C(w)~1 2p Z PRC (w). s is the intensity of the noise, rescaled with the power of the PRC [36]. The effect of the white noise term ss i (t) can be regarded as temporal fluctuations in the firing rate of each postsynaptic neuron around the mean v (i.e., natural frequency). Moreover, if the fluctuation strength of the measured PRC is constant with respect to the perturbation timing, the effect of ss i (t) can also be regarded as a fluctuation of the measured PRC [37].
Here, we expand the STDP window function and the coupling function into their Fourier series: where a k and b l are the Fourier coefficients of the STDP window function and the coupling function (PRC), respectively. V(h) and C(w), which are real-valued functions, satisfy a {k~a Ã k and b {l~b Ã l (the superscript * denotes the complex conjugate), respectively. The parameters k (~0,+1,+2, Á Á Á ) and l (~0,+1,+2, Á Á Á ) denote the wavenumbers. Note that the initial synaptic weight J 0 can be involved in the DC component a 0 without loss of generality. Here, we define A k and f k as the amplitude and the phase of a k , respectively (a k~Ak exp (if k )). B l and x l represent the amplitude and the phase of b l (b l~Bl exp (ix l )). They satisfy the following equations: The order parameter m m k,l , the overlap between the k-th frequency component of the m-th memory key pattern g m and the l-th frequency component of the retrieval key pattern y, is defined as Because g m and y do not vary over time, each overlap m m k,l takes a constant value. Note that all the postsynaptic neurons share the same order parameter m m k,l and that m m k,l represents the characteristic function of presynaptic phase disrtibution P(y j jag m j ) at each wavenumber: m m k,l~Ð 2p 0 dyP(yjag m ) exp (i(kg m {ly)). By using a k , b l , and m m k,l , the LPE (7) can be transformed into Because the postsynaptic neurons share the same m m k,l and are driven by independent noise, the M neurons can be considered to be statistically independent of each other and have the same statistical characteristics.

Equilibrium Distribution when Storing a Finite Number of Patterns
For mathematical simplicity and tractability, we consider the case when the number of stored paired patterns p is finite in the thermodynamic limit as N?? (i.e., pvvN). Given a retrieval key pattern similar to ag 1 , ah 1 is to be retrieved. As described above, the memory key patterns g m (m~1,2, Á Á Á ,p) are not correlated with each other. The same goes for the memory output patterns h m .
The retrieval key pattern is generated with the following von Mises probability density function (PDF): where ag 1 j corresponds to the mean of this PDF. c is a measure of the concentration, and it controls the similarity between the retrieval key pattern and ag 1 .
Under the above definition of retrieval key pattern generation, each overlap m m k,l (m~1,2, Á Á Á ,p) can be calculated as follows. The average overlap with m~1 between pairs of frequency components which satisfy k~al is where I l (c)~1 2p Optimal Design for Hetero-Associative Memory PLOS ONE | www.plosone.org function of the first kind. On the other hand, the average overlap with m~1 between the other components (k=al) is and the deviation is O(1= ffiffiffiffi ffi N p ). Moreover, for any k and l, the overlap between the retrieval key y and a memory key pattern other than the first one (m §2) is on average and the deviation is O(1= ffiffiffiffi ffi N p ). In the case of a finite p, the number of terms with m §2 in Eq. (13) is finite. Moreover, the Fourier coefficients a l and b l rapidly approach zero as l increases. Thus, contributions of the terms except those with m~1 and k~al to Eq. (13) can be neglected in the limit N??. The phase dynamics can be rewritten as Note that the term consisting of DC components a 0 and b 0 in Eq. (13) is a constant. We can safely neglect this term, because the constant term can be involved in the natural frequency v without loss of generality. Equation (18) shows that the statistical properties of the hetero-associative memory model in the case of finite loading (pvvN) and those in the simplest case of just one pair of patterns to be stored (p~1) are identical. From Eq. (18), we obtain the equilibrium phase distribution of each postsynaptic neuron: where Z NF is the normalizing factor.

Mutual Information per Neuron
Mutual information, which measures how related two random variables are, is nonnegative and takes 0 only if these variables are independent.
When recalling the a times spread-out memory pattern ah 1 , the mutual information of the retrieval output w relative to ah 1 per neuron, 1 M H(w; ah 1 ), is given by Here, 1 M H(w) is the entropy of w per neuron, which measures the uncertainty associated with w. 1 M H(wjah 1 ) is the conditional entropy of w given ah 1 , which quantifies the remaining uncertainty of w given that ah 1 is known. Because each postsynaptic neuron is statistically independent of each other and has the same statistical characteristics as described above, 1 M H(w) and 1 M H(wjah 1 ) can be simply written as follows: . This is because the system is symmetric with respect to sign inversion of the key and output patterns. Namely, the system, which can retrieve normal and doubly spread-out patterns, can also retrieve reversed patterns with the same quality.
In what follows, we search for pairs of PRCs and STDP window functions that are optimal for retrieving both normal and doubly spread-out patterns by jointly maximizing two object functions, 1 M H(w; h 1 ) and 1 M H(w; 2h 1 ). To solve this joint optimization problem, we employ a simple sum of these functions, Furthermore, for comparison, we also use an objective function, As mentioned above, the optimal system derived by maximization of the objective function is also optimal for retrieving reversed patterns.

Phase Response Curve of Hippocampal CA1 Pyramidal Neurons
In our previous work, we obtained PRCs from rat hippocampal CA1 pyramidal neurons by performing whole-cell patch-clamp recording in vitro [18,19]. Figure 2B shows the PRCs from that research. In the protocol for measuring PRCs, we injected DC depolarizing currents into the somata of rat CA1 pyramidal neurons to evoke periodic firing. Using the dynamic clamp, the mean inter-spike interval (ISI) was adjusted to the target period by tuning the DC depolarizing current. In measuring these PRCs, the firing frequencies of those neurons were tuned in the theta frequency range (4-14 Hz). Next, we evoked a one-shot rectangle perturbation superimposed on the DC depolarizing current by using various timings relative to the spike, and measured how the perturbation current disturbed the timing of the succeeding spike, i.e., phase response. The spike times randomly fluctuated due to intrinsic noise in the neurons. To extract the PRCs from the noisy data of the phase responses, we made use of a PRC measurement model formulated as an LPE [38] (the same as the one used in the current study) and applied a maximum a posteriori (MAP) estimation algorithm based on the measurement model to the noisy phase response data. The effectiveness of the measurement model and the reliability of the estimated PRCs were verified by demonstrating that the LPEs with the estimated PRCs could predict the stochastic behaviors of the same neurons, whose PRCs had been measured, when they were perturbed by various periodic stimulus currents [18]. A detailed explanation of the experimental conditions and the MAP estimation algorithm can be found in [36,38,39], while the reliability and quality of PRCs used here has been discussed in detail in [18,19].

Performance of Hetero-associative Memory Model with Typical Parameters
We carried out numerical simulations on the hetero-associative memory model with typical parameters, and we compared the numerical results with theoretical predictions. Here, we use the hetero-associative memory model endowed with a PRC (cell #1 in Fig. 2B [18]) and a typical STDP window function (left panel of Fig. 2C [16]) of CA1 pyramidal neurons. In the following numerical simulations, we embedded three pairs of random phase patterns, g m and h m (m~1,2,3) in the synaptic weight J ij . We used the following to evaluate the retrieval quality of this model: This measure is the overlap between the k-th frequency component of the m-th memory output pattern h m and the l-th frequency component of the retrieval output pattern w. At equilibrium, the overlap M m k,l can be theoretically obtained with the PDF (19): As shown in Eq. (27), M m k,l represents the characteristic function of the postsynaptic phase distribution P(wjah m ) at each wavenumber.
First, we verified the effect of intrinsic noise on the retrieval quality of this model. For simplicity, we gave it a retrieval key pattern identical to the normal memory key pattern (y~g 1 ), which corresponds to the special case of c??. Figure 3A plots the amplitude of the overlap jM 1 l,l j (l~1,2, Á Á Á ,5) at equilibrium as a function of the noise intensity s. In this figure, the LPE (13) was solved numerically by using the Euler method, and the values of jM 1 l,l j calculated with Eq. (26) were compared with theoretical predictions obtained by Eq. (27). As shown in Fig. 3A, the numerical results coincide with the theoretical values for all wavenumbers l. When the noise intensity is sufficiently small (sv0:05), the retrieval output pattern w has an appreciable overlap with h 1 , i.e., M 1 1,1 *O(1). As s increases, the overlap M 1 l,l approaches zero. Figure 3B shows an example of the PDF (19) and a histogram of the phase difference w i {h 1 i obtained by numerically solving LPE (13) at equilibrium. Here, s~0:03 and c??. In this figure, the PDF (19), which forms a unimodal distribution, is in good agreement with the histogram normalized by the bin width.
Next, we verified the effect of the degraded retrieval key patterns on the retrieval quality of this model. We generated the phase patterns by using the conditional PDF (14) given g 1 and various values of c, and we used these generated patterns as retrieval key patterns. Figures 3C and D respectively show amplitudes of the overlaps jm 1 l,l j (defined in Eq. (12)) and jM 1 l,l j (l~1,2, Á Á Á ,5) at equilibrium as a function of c. Here, s~0:03. The numerical results coincide with the theoretical ones for all l. When c is sufficiently large (cw3), the retrieval output pattern w has an appreciable overlap with h 1 , i.e., M 1 1,1 *O(1). As c decreases, the overlap m 1 l,l approaches zero faster than M 1 l,l converges to zero. Note that we got similar results to those above by using the other PRCs (different from cell #1 in Fig. 2B) and another STDP window function (right panel of Fig. 2C).

STDP Window Functions Optimally Matched to PRCs of Hippocampal CA1 Pyramidal Neurons
We searched for STDP window functions optimally matched to the PRCs of the five hippocampal CA1 pyramidal neurons shown in Fig. 2B. As described in the Methods section, we considered the two cases. One is that we maximize the objective function I total defined in Eq. (25) to search for STDP window functions that are optimal for retrieving normal patterns. The other is that we maximize the objective function I total defined in Eq. (24) to search for STDP window functions that are optimal for retrieving normal and doubly spread-out patterns. In both cases, we assigned the Fourier coefficients of the PRCs of the hippocampal CA1 pyramidal neurons to b l of each mutual information constituting the objective function I total , and under the constraint of the measured PRC, searched for a l , the Fourier coefficients of the STDP window functions to maximize the objective function I total . Note that a system which can optimally retrieve normal and doubly spread-out patterns can also retrieve reversed ones, because it is symmetric with respect to sign inversion of the key and output patterns.
Because the mutual information 1 M H(w; ah 1 ) in Eq. (21) monotonically increases as A l (the amplitude of a l ) increases, we imposed the following constraint condition on the power of STDP window function: Here, we truncated the Fourier series of the PRCs and STDP window functions after the fifth term. Referring to the value of the STDP power in a previous study [10], we set Const:~0:12. To solve this optimization problem, we used the FMINCON function in the Matlab Optimization Toolbox.
First, we searched for STDP window functions that are optimal for retrieving normal patterns. By maximizing the objective function in Eq. (25) under the power constraint (28), we obtained a connected set of optimal STDP window functions described as V(g)~2 Next, we searched for STDP window functions that are optimal for both retrieving normal and doubly spread-out patterns. By maximizing the objective function in Eq. (24) under the power constraint (28), we obtained a connected set of optimal STDP window functions in the same manner as described above. Figures 4A9-D9 show five sets of optimal STDP window functions; each set is optimally matched to each PRC of the five different CA1 pyramidal neurons shown in Fig. 2B. The four panels of Figs. 4A9-D9 plot examples of optimal STDP window functions with different phases j. All sets of functions except for cell #1 have the same form, and they are almost completely composed of fundamental and second frequency components. The shape of the optimal STDP window function obtained from the PRC of cell #1 is very similar to the others, even though the PRC of cell #1 contains a relatively large number of higher frequency components compared with the PRCs of the other cells (see Fig. 2B).
The results of Fig. 4 suggests that the optimal STDP window function depends heavily on the fundamental frequency components of the PRCs, and thus its shape is nearly invariant for all of the PRCs of the CA1 pyramidal neurons. In addition, the joint optimization for retrieving normal and doubly spread-out patterns only yielded the second frequency components of the STDP window functions, and thus, the second frequency components of the STDP window function play a key role in recalling doubly spread-out phase patterns.
It has been reported that there are two types of STDP window functions in hippocampal CA1 pyramidal neurons [13,14], i.e., symmetric (left panel of Fig. 2C [16]) and asymmetric (right panel of Fig. 2C [17]). Here, we compared physiologically measured window functions with purposely designed ones for memory recalls. We computed the Fourier series of symmetric and asymmetric STDP window functions in Fig. 2C and compared the fundamental and second frequency components of the STDP window functions in Fig. 2C with the frequency components of the ones in Figs. 4A9-D9. Figure 5A plots symmetric and asymmetric STDP window functions composed of only the fundamental and second frequency components of the ones in Fig. 2C. The purposely designed STDP window functions shown in Figs. 4A9 and B9 qualitatively conform to those of Fig. 5A. Figure 5B shows the rate of the fundamental and second frequency components for STDP window functions in Fig. 5A and the purposely designed ones in Figs. 4A9-D9. We compared the amplitudes between the two Fourier coefficients of each STDP window function: . As shown in Fig. 5B, the joint optimization for retrieving normal and doubly spread-out patterns yields equal amounts of fundamental frequency component and second frequency component (A 1~A2 ), and the amount of second frequency component in the symmetric and asymmetric STDP window functions in Fig. 2C is almost equal to that of the fundamental frequency components. Moreover, Figs. 4C9 and D9 show an inverted symmetric window function and inverted asymmetric one in contradistinction to Figs. 4A9 and B9. These window functions were found in regions outside the hippocampal CA1 area (see [40][41][42]).

Memory Retrieval in the Hetero-associative Memory Model
By using numerical simulations, we confirmed that the system with the STDP window functions in Figs. 4A9-D9 can function as intended. The synaptic weight J ij was determined using the STDP window function (cell #5 in Fig. 4A9) to store three pairs of random phase patterns, g m and h m (m~1,2,3), and the retrieval performance of the system with the determined synaptic weight and the measured PRC (cell #5 in Fig. 2B) was verified. In the following simulations, we used the retrieval key pattern generated with the conditional PDF (Eq. (14)) given ag 1 , and under this condition, we checked whether the system could recall a temporally reversed memory output pattern and/or one stretched to jaj times its original timescale, ah 1 . The overlap M m k,l between the k-th frequency component of h m and the l-th frequency component of w (defined in Eq. (26)) was used as a measure of retrieval performance.
The three panels of the left column in Fig. 6 show the time evolution of jM m k,l j in the cases of the normal, reversed, and spread-out pattern retrievals. jM 1 1,1 j~O(1) and the others are almost zero in the normal pattern retrieval, whereas jM 1 {1,1 j~O(1) and the others are almost zero in the reversed pattern retrieval, and jM 1 2,1 j~O(1) and the others are almost zero in the doubly spread-out pattern retrieval. The panels of the center and right columns in Fig. 6 show samples of memory output patterns and retrieval output patterns at equilibrium in the cases of the normal, reversed, and doubly spread-out pattern retrievals. We also confirmed that the system with the other STDP window functions in Figs. 4B9-D9 and the other PRCs (different from cell #5 in Fig. 2B) works just as well as these.

Summary of Results and Conclusions
By maximizing the objective functions given by the mutual information, we derived pairs of STDP window functions and PRCs optimally recalling normal, reversed, and doubly spread-out phase patterns in a hetero-associative memory model. We searched for a set of optimal STDP window functions using the measured PRCs from in vitro experiments on hippocampal CA1 pyramidal neurons.
The optimal STDP window function heavily depends on the fundamental frequency component of the PRCs, and thus its shape is almost invariant with respect to the PRCs of CA1 pyramidal neurons (see Fig. 4). Even though the PRC of cell #1 contains a relatively large number of higher frequency components compared with the PRCs of the other cells (see Fig. 2B), the shape of the optimal STDP window function obtained from the PRC of cell #1 is very similar to the others. A comparison of results in Figs. 4A-D and A9-D9 suggests that second frequency components of STDP window functions play a key role in recalling doubly spread-out phase patterns. If the memory key and output patterns are respectively assigned to an independent number in ½0,2p) with a uniform probability, in the limit of N??, the doubly spreadout patterns are orthogonal to the original patterns (i.e., Ð 2p 0 dg exp (i(g{2g))~0). Because of this orthogonality, this system can retrieve spread-out patterns if the STDP window function contains higher frequency components. As shown in Fig. 6B, the system, which can retrieve normal and doubly spread-out phase patterns, can also retrieve reversed patterns. This is because the symmetry with respect to sign inversion for key and output patterns is satisfied as it is in the conventional associative memory model [43,44]. Thus, the mutual information in retrieving normal and doubly spread-out phase patterns is equal to the information involved in retrieving the reversed patterns.
Furthermore, the mutual information is invariant with respect to the phases of the STDP window functions and PRCs. Thus, the set of optimal STDP window functions forms a connected set homeomorphic to a ring, examples of which have a good qualitative match to those reported in hippocampal CA1 pyramidal neurons, such as the symmetric [15,16] and asymmetric window functions [17]. Note that the original data from Wittenberg and Wang (2006) exhibit a phase delay between the pre-and postsynaptic spikes in the peak of the symmetric STDP window function [16]. On the other hand, we used a simplified STDP window function (left panel of Fig. 2C) that ignored the phase delay in the peak. This is because the phase shift of the STDP window function has no effect on the retrieval performance of the associative memory model, as stated above. As shown in Fig. 5B, the fundamental and second frequency components of STDP window functions reported for CA1 neurons have roughly the same scale, which coincides with those of the theoretically derived STDP window functions.
Thus, the results obtained here suggest that the properties of the neural interaction and the synaptic plasticity rule in the CA1 region support a hetero-associative memory function recalling normal, doubly spread-out, and reversed patterns.

Effect of STDP Multiplicity in Single Neurons in Recalling Memories
Optical imaging studies have suggested that the shape of STDP window function in the CA1 pyramidal neuron depends on the location on the dendrite [13,14]. A symmetric STDP window function was observed in the proximal-to-soma dendrite, whereas an asymmetric STDP window function was observed in the distalto-soma dendrite.
Here, we verify the effects of symmetric and asymmetric STDP window functions in single neurons on recall memory. We assume that synapses between a postsynaptic neuron i (~1,2, Á Á Á ,M) and a presynaptic neuron j (~1,2, Á Á Á ,N s ) obey the symmetric STDP rules, and others between a postsynaptic neuron i (~1, Á Á Á ,M) and a presynaptic neuron j (~N s z 1,N s z2, Á Á Á ,N) obey the asymmetric STDP rules. N s is the number of synapses obeying the symmetric STDP rules in a single neuron. For mathematical simplicity, the symmetric and asymmetric STDP window functions coexisting in single neurons are described as Here, the variable j is the phase difference between the symmetric and asymmetric STDP window functions. The typical value of j is p=2. Under this condition, we can rewrite Eq. (13) as follows:  Figure 6. Confirmation that the system with the STDP window functions in Figs. 4(A9-D9) can function as intended. The synaptic weight J ij was determined using the STDP window function (cell #5 in Fig. 4A9) to store three pairs of random phase patterns, g m and h m (m~1,2,3), and when presented with the retrieval key pattern generated with the conditional PDF (Eq. (14)) given ag 1 , the retrieval performance of the system with the determined synaptic weight and the measured PRC (cell #5 in Fig. 2B) was verified by using numerical simulations (M~N~1000, c~20, s~0:03). where m m k,l is the same as the order parameter defined in Eq. (12). Thus, the above model is essentially equivalent to the homogeneous STDP model (Eq. (13)) except for the existence of the coefficient C k and the phase r k . If the fundamental frequency components are dominant in the PRCs, this inhomogeneous model can also work as a hetero-associative memory. Figure S1 (supporting information) shows the results of numerical simulations when l~0:5. Because C k v1 if l=1 or l = 0, the retrieval quality becomes worse than that of the homogeneous model. Lengyel et al. (2005) tried to determine whether the properties of neural interactions and the synaptic plasticity rule in the CA3 region support an auto-associative memory function [10]. Figure  S2 (supporting information) illustrates the top-down approach they used. First, they developed a theory treating auto-associative memory retrieval as a kind of Bayesian inference and constructed a gradient ascent algorithm for the MAP estimation. Next, they reinterpreted this algorithm as phase oscillators consisting of PRCs and STDP window functions. Finally, they qualitatively compared the PRCs of hippocampal CA3 pyramidal neurons with ones theoretically derived from a typical STDP window function. As a result of their contrasting the top-down approach with Marr's trilevel hypothesis [45], their study can be considered as a bridge from the computational level to the algorithmic level. Furthermore, they tried to bridge the algorithmic and physical levels by reinterpreting the algorithm they derived as phase oscillators.

Difference between Our Approach and Lengyel's
On the other hand, contrasting our approach shown in Fig. 1 with Marr's tri-level hypothesis, our study can be considered as a bridge between the physical level and the algorithmic level. Note that the phase equation (e.g., Eq. (13)) reduced from the coupled oscillator system (Eq. (3)) corresponds to the algorithm of heteroassociative memory. Unfortunately, there is no clear relationship between our phase equation derived from the bottom up and Lengyel's gradient ascent algorithm derived from the top down at this moment. In the future, we will explore the correspondence between our results and Lengyel's. [10].
How and Where are Key Patterns for Recalling Reversed and Spread-out Patterns Created?
As shown in Fig. 2, a key pattern has to be input in order to recall its associated output pattern in a hetero-associative memory. As summarized in Table 1, to retrieve reversed and spread-out patterns, the associated key patterns also have to be reversed and spread-out. Thus, we must answer new questions as to where and how key patterns for recalling reversed and spread-out patterns are created. It is possible for a recurrent network such as the CA3 network to create and provide reversed and spread-out patterns. The CA3 region provides one of the dominant inputs to the CA1 region [25,27]. By applying mean field approximations, our theory of hetero-associative memory can be straightforwardly extended to analyses of auto-associative memory. The preliminary results suggest that the properties of sine-like PRCs of hippocampal CA3 pyramidal neurons recorded by Lengyel et al. (2005) [10] and the typical STDP rule can support an auto-associative memory function for recalling reversed and spread-out phase patterns. Thus, the preliminary results and the results of this paper indicate the possibility that a combination of the CA1 network and the CA3 network can consistently work to retrieve reversed and spread-out patterns. We will report on this issue in our next study.

Role of Reversed and Spread-out Pattern Retrievals
It has been considered that memories are first stored in the hippocampus and are gradually moved to the neocortex in a more permanent form of storage. Temporally spread-out pattern retrieval, in which the temporal order of the memory spike sequence is preserved and the timescale of retrieval pattern is about two times longer, may be important for the memory translation and system consolidation [8]. On the other hand, temporally reversed pattern retrieval is suggestive of evaluating event sequences in the manner of reinforcement learning models [46]. During waking periods, reversed pattern retrieval occurs in situ, allowing immediately preceding events to be evaluated in precise temporal relation to the current one, and so it may be an integral mechanism for learning about recent events [9]. Figure S1 Effect of coexisting symmetric and asymmetric STDP window functions in single neurons on the memory retrieval. In this numerical simulation, we used the PRC (cell #2 in Fig. 2B) and symmetric and asymmetric STDP window functions (cell #2 in Figs. 4A9 and 4B9) at the same rate in single neurons. M~N~1000, N s~5 00, c~20, s~0:01.