text:
# Improved Wilson QCD simulations at light quark masses

## 1 SIMULATION PERFORMANCE

Since 1997, UKQCD has been conducting dynamical fermion simulations on a Cray T3E machine using typically $`128`$ processors ($`35`$ Gflops sustained). This has permitted production studies on lattices up to $`16^3\times 32`$ and $`m_\pi /m_\rho `$ down to around 0.60 . Study groups within the collaboration are now planning physics programmes for the advent of the QCDOC multi-Teraflops machine which is expected to be our main resource from 2003 onwards. Actions and algorithmic performance are key concerns.

### 1.1 HMC performance

The next figure shows estimates of the work (in Teraflops) required to obtain a ‘new’ configuration using three of the main methods currently available.

The data for the plot was mostly taken from that provided last year at LAT01 for the session on algorithm performance . The top curve corresponds to the UKQCD ‘clover-Wilson’ estimate:

$$\frac{\text{Gflops}}{\text{config}}=0.157\left(\frac{L}{a}\right)^{3.41}\left(\frac{T}{a}\right)^{1.14}\left(\frac{1}{am_\pi }\right)^{2.77}.$$
(1)
The accompanying curves correspond to Wilson fermions as measured by the SESAM collaboration  and the Asqtad MILC code (‘order $`a^2`$ tadpole’ improved staggered fermions using the R-algorithm ). The dashed curve corresponds to the Asqtad data but where the lattice spacing has been increased by $`50\%`$. The Asqtad scaling curve was taken from the performance data in . We performed an independent check by carrying out Asqtad simulations with 2 flavours on the same size lattice and machine used to produce the UKQCD data. For this comparison, a ‘new configuration’ amounts to 20 MD trajectories for the Wilson and clover HMC runs and 6 for the Asqtad runs. These numbers are based on using at least twice the measured integrated plaquette autocorrelation time in each case. The motivation for the extra comparison of Wilson at lattice spacing $`a`$ with Asqtad at spacing $`1.5a`$ is the observation that the improvement for the latter is much greater ($`𝒪(a^2\alpha _s`$)) than for Wilson at a given spacing $`a`$. Since the work scales like $`a^8`$ this contributes very significantly to the saving in work at fixed physical lattice size and $`m_\pi `$.

The Wilson/SESAM simulation  appears to be $`8`$ times faster than the improved Wilson/UKQCD one. Detailed comparison shows that a factor of $`3.4`$ comes from solver iterations (better preconditioning and solver algorithm and better underlying condition number); a factor of $`1.8`$ from larger MD step size (higher acceptance, physics); a factor of $`1.3`$ from the numerical overhead of the clover term. These naive comparisons do not take into account other factors such as relative scaling behaviour and accessibility and efficiency of measurement operators.

### 1.2 HMC stability

These HMC runs are expected to be susceptible to instabilities (occasional large $`\mathrm{\Delta }H`$ values and zero acceptance) when the fermion force term gets too large . We have observed this effect directly on the clover-Wilson run at $`\kappa =0.1358`$ where we find it necessary to use a step size of $`1/400`$. When the step size is this, or smaller, we find that $`64`$ bit for our field storage and matrix-vector manipulations is required so as to avoid a dramatic loss of acceptance due to rounding errors. This is so even though we always use full 64-bit arithmetic and sophisticated summing techniques in the global summing of our $`\mathrm{\Delta }H`$ (energy difference) calculations.

### 1.3 Run parameters

The main parameters of the new light quark run are shown in the following table. In physical units, $`L1.5fm`$ and $`m_\pi 430`$ Mev before extrapolation.

| $`(\beta ,\kappa )`$ | $`(5.2,.1358)`$ |
| --- | --- |
| trajectories/hour | $`3`$ |
| no. trajectories | $`2440`$ |
| no. configs | $`122`$ |
| $`\tau _{AC}^{\text{int}}`$(plaquette) | $`6.9(6)`$ |
| $`r_0/a`$ | $`5.3(5)`$ |
| $`a`$ (using $`r_0=0.49`$ fm) | $`0.0925(9)`$ fm |
| $`am_{PS}`$ | $`0.207(5)`$ |
| $`m_{PS}/m_V`$ | $`0.43(2)`$ |
| $`m_{PS}L`$ | $`3.3(1)`$ |

The integrated autocorrelation time for the plaquette follows the previously observed trend that it decreases with decreasing quark mass. For comparison, $`\tau _{AC}^{\text{int}}`$=$`16(3)`$ at $`\kappa =.1350`$ where $`m_{PS}/m_V0.7`$. Although not expected, this trend can be accommodated in simple models .

The ratio ‘$`m_\pi /m_\rho `$’ is below the threshold for the continuum decay $`\rho \pi \pi `$ but not for decay on this size of lattice. We are using maximum entropy methods  to look for evidence of decay of the lattice state at the relevant energy.

## 2 MEASUREMENTS

### 2.1 Hadron masses and finite size effects

Since $`m_\pi L=3.3(1)`$ (significantly below $`5`$) we anticipate finite size effects. For small and intermediate size boxes there is evidence  that hadron mass corrections follow a power law ($`L^3`$) rather than exponential ($`\mathrm{exp}(L/L_0)`$) size dependence. A simple phenomenological model for this might be

$$\frac{\delta m_h}{m_h}=\frac{c_h}{(m_\pi L)^3}.$$
(2)
The dimensionless parameter $`c_h`$ is presumably related to the form factor of the hadron $`h`$.

The figure below shows a compilation (from UKQCD and QCDSF ) of pseudoscalar and vector masses (fixed $`\beta =5.2`$ and $`\kappa _s=\kappa _v`$). The filled points are the raw masses extracted from multi-channel fits.

Aside from the point at very heavy quark mass ($`\kappa =0.1342`$), the $`m_{PS}`$ data (filled circles) are well fitted by lowest order chiral PT (linear in $`1/\kappa `$).

We have estimated $`c_\pi 6`$ and $`c_\rho 7`$ from the finite size (FS) effects reported by the JLQCD Collaboration, who used the same action (clover-Wilson)  on a variety of lattice sizes. We have checked that the corrections implied by these values are consistent with the FS effects which UKQCD has previously measured using a slightly different action . Using this model to ‘correct’ the measured values, we get the open points shown in the figure.

The FS-corrected pseudoscalar points are still well fitted by a straight line and yield $`\kappa _{crit}=0.135978(9)(55)`$. The first error is statistical and the second systematic - estimated by including the heaviest quark point in a quadratic fit. Taking the above FS corrections for the pseudoscalar and vector into account we estimate $`m_{PS}/m_V0.45`$ for the lightest quark simulation.

### 2.2 Static potential

The next figure shows the difference of the measured static potential from the universal bosonic string model for all UKQCD data sets at fixed $`\beta =5.2`$  including this latest run (stars) at $`\kappa =0.1358`$ and a QCDSF set at $`\kappa =0.1340`$ .

Also shown is the corresponding result from a modest run we have made using the MILC Collaboration’s code at $`a0.130`$ fm . Violations of rotational symmetry are marked for $`r<0.7r_00.35`$ fm for all these coarse lattices. Although the improved staggered fermion measurements (Asqtad) are at a coarser lattice spacing, they show smaller violations of rotational symmetry than the clover-Wilson data. This is presumably due to the significantly improved gluonic action. The Asqtad data (at $`m_{PS}/m_V0.50`$) show a greater effect of the running coupling in the Coulomb term at short distances than the corresponding clover-Wilson data (lower at short distances). Within the fixed $`\beta `$ clover-Wilson set there is no very clear trend with decreasing quark mass. The corresponding plot for clover-Wilson data at fixed lattice spacing  shows more clearly the effect of the running coupling.

## 3 CONCLUSIONS

We have analysed data from more than 2200 trajectories of a clover-improved Wilson simulation at $`m_q/m_s0.3`$. We have estimated finite volume corrections to meson masses. Light quark masses are accessible using clover-improved Wilson fermions and standard HMC but simulations are slow and can suffer from instabilities. Improved staggered fermions appear to offer a more efficient route to configurations with light dynamical quarks, but physics measurements are in general more complicated.

## Acknowledgements

I thank my UKQCD colleagues in the dynamical fermions project, particularly Derek Hepburn, Craig McNeile and Steve Miller. I thank Thomas Lippert and Steve Gottlieb for detailed insights into their algorithm performance figures, and the MILC Collaboration for the use of their simulation code.


format: text

text:
# Robustness of retrieval properties against imbalance between long-term potentiation and depression of spike-timing-dependent plasticity

## I Introduction

Recent experimental finding indicates that synaptic modification in cortical neurons depends on the precise temporal relationship between pre- and post-synaptic spikes Bi and Poo (1998); Markram et al. (1997); Zhang et al. (1998). In particular, pre-synaptic spikes that precede post-synaptic firing induce long-term potentiation (LTP) by no more than $`20`$ ms, while those that follow post-synaptic firing induce long-term depression (LTD), with a rapid transition (a few ms). The magnitude of synaptic modification decays exponentially with the time interval between pre- and post-synaptic spikes. This form of synaptic modification has been called spike-timing-dependent plasticity (STDP) Song et al. (2000) or temporally asymmetric Hebbian learning (TAH) Abbott and Song (1999); Rubin et al. (2001).

The functional role of STDP has been investigated by many authors. They showed that STDP is a mechanism for synaptic competition Abbott and Song (1999); Song et al. (2000); Levy et al. (2001); van Rossum et al. (2001); Rubin et al. (2001); Song and Abbott (2001) or a learning mechanism of sequential patterns Gerstner et al. (1996); Kempter et al. (1999); Kistler and van Hemmen (2000); Rao and Sejnowski (2000); Munro and Hernandez (2000); Yoshioka (2002); Matsumoto and Okada (2002). Asymmetric learning window depending on spike timing like STDP has been studied and shown to be appropriate learning rule for sequential patterns Herz et al. (1989); Gerstner et al. (1993); Abbott and Blum (1996). However, this asymmetric learning rule does not involve LTD. Some authors showed that the balance between LTP and LTD of STDP is significant for a stable learning Song et al. (2000); Munro and Hernandez (2000); Yoshioka (2002); Matsumoto and Okada (2002). In our previous work, we analytically showed that STDP has the same qualitative effect as the covariance rule when the spatio-temporal patterns are stored since the differences between spike times that induce LTP or LTD are capable of canceling out the effect of the firing rate information Matsumoto and Okada (2002). In the brain, a situation in which the balance is maintained precisely is inconceivable. The data points obtained by experiments are fluctuated in the different trials Bi and Poo (1998); Zhang et al. (1998). Therefore, it is meaningful to discuss more biological plausible situation to investigate the neuronal mechanism for sequential learning in the brain. Some authors numerically investigated the impact of the imbalance between LTP and LTD on the network properties Song et al. (2000); Munro and Hernandez (2000).

The aim of this paper is to analytically discuss the retrieval properties of spatio-temporal patterns in an associative memory model that incorporates the imbalance between LTP and LTD of STDP using a method of the statistical neurodynamics Matsumoto and Okada (2002); Okada (1996); Amari and Maginu (1988). According to our previous work, when the balance is not precisely maintained, it is impossible to cancel out the information of firing rate. Consequently, a cross-talk noise diverges. However, if the magnitudes of LTP and LTD are equivalent on average in a learning process, it may be possible to stably retrieve spatio-temporal patterns. Since the ratio of LTP and LTD is crucial, the magnitude of LTD changes while that of LTP is fixed. We found that the storage capacity takes a finite value even at large $`\delta `$ when the fluctuation of LTD is assumed to obey a Gaussian distribution with mean $`0`$ and variance $`\delta ^2`$. This implies that the balance between LTP and LTD of STDP need not to be maintained precisely, but must be maintained on average. This mechanism might work in the brain. Furthermore, we found that a basin of attraction becomes smaller as $`\delta `$ increases while an initial critical overlap remains unchanged.

## II Model

The model contains $`N`$ binary neurons with reciprocal connections. Each neuron has a binary state $`\{0,1\}`$. We define discrete time steps and the following rule for synchronous updating:

$`u_i(t)`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{N}{}}}J_{ij}x_j(t),`$ (1)
$`x_i(t+1)`$ $`=`$ $`F(u_i(t)\theta ),`$ (2)
$`F(u)`$ $`=`$ $`\{\begin{array}{cc}1.\hfill & u0\hfill \\ 0.\hfill & u<0,\hfill \end{array}`$ (5)
where $`x_i(t)`$ is the state of the $`i`$-th neuron at time $`t`$, $`u_i(t)`$ is the internal potential of that neuron, and $`\theta `$ is a uniform threshold. If the $`i`$-th neuron fires at $`t`$, its state is $`x_i(t)=1`$; otherwise, $`x_i(t)=0`$. $`J_{ij}`$ is the synaptic weight from the $`j`$-th neuron to the $`i`$-th neuron. Each element $`\xi _i^\mu `$ of the $`\mu `$-th memory pattern $`𝝃^\mu =(\xi _1^\mu ,\xi _2^\mu ,\mathrm{},\xi _N^\mu )`$ is generated independently by

$$\text{Prob}[\xi _i^\mu =1]=1\text{Prob}[\xi _i^\mu =0]=f.$$
(6)
The expectation of $`𝝃^\mu `$ is $`\text{E}[\xi _i^\mu ]=f`$, and thus $`f`$ is considered to be the mean firing rate of the memory pattern. The memory pattern is sparse when $`f0`$, and this coding scheme is called sparse coding.

The synaptic weight $`J_{ij}`$ follows the form of synaptic plasticity, which depends on the difference in spike times between the $`i`$-th (post-) and $`j`$-th (pre-) neurons. The time difference determines whether LTP or LTD is induced. This type of learning rule is called spike-timing-dependent plasticity (STDP). The biological experimental findings show that LTP or LTD is induced when the difference in the pre- and post-synaptic spike times falls within about $`20`$ ms Zhang et al. (1998). We define a single time step in equations (15) as $`20`$ ms, and durations within $`20`$ ms are ignored. The learning rule based on STDP conforms to this equation:

$$J_{ij}=\frac{1}{Nf(1f)}\underset{\mu =1}{\overset{p}{}}\left\{\xi _i^{\mu +1}\xi _j^\mu (1+ϵ_{ij}^{\mu 1})\xi _i^{\mu 1}\xi _j^\mu \right\}.$$
(7)
The number of memory patterns is $`p=\alpha N`$, where $`\alpha `$ is defined as a loading rate. LTP is induced when the $`j`$-th neuron fires one time step before the $`i`$-th neuron, $`\xi _i^{\mu +1}=\xi _j^\mu =1`$, while LTD is induced when the $`j`$-th neuron fires one time step after the $`i`$-th neuron, $`\xi _i^{\mu 1}=\xi _j^\mu =1`$. Since the ratio of LTP and LTD is crucial, the magnitude of LTD changes while the magnitude of LTP and the time duration are fixed. $`ϵ_{ij}^\mu `$ is generated independently and obeys a Gaussian distribution with mean $`ϵ`$ and variance $`\delta ^2`$: $`ϵ_{ij}^\mu 𝒩(ϵ,\delta ^2)`$. Fig.1 shows the time function of STDP in our model. When $`ϵ_{ij}^\mu =0`$, the balance between LTP and LTD is precisely maintained and then the model is equivalent to the previous model Matsumoto and Okada (2002). A sequence of $`p`$ memory patterns is stored by STDP: $`𝝃^1𝝃^2\mathrm{}𝝃^p𝝃^1\mathrm{}`$. In other words, $`𝝃^1`$ is retrieved at $`t=1`$, $`𝝃^2`$ is retrieved at $`t=2`$, and $`𝝃^1`$ is retrieved at $`t=p+1`$. There is a critical value $`\alpha _C`$ of the loading rate, so that the loading rate larger than $`\alpha _C`$ makes retrieval of the pattern sequence unstable. $`\alpha _C`$ is called a storage capacity.

## III Theory

In this section, we derive dynamical equations that describe the retrieval properties of the network. In this paper, we consider the thermodynamic limit: $`N\mathrm{}`$. The $`i`$-th neuronal internal potential $`u_i(t)`$ at time $`t`$ can be expressed as

$`u_i(t)`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{N}{}}}J_{ij}x_j(t)={\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\mu =1}{\overset{p}{}}}\left\{\xi _i^{\mu +1}\xi _j^\mu (1+ϵ_{ij}^{\mu 1})\xi _i^{\mu 1}\xi _j^\mu \right\}x_j(t)`$ (8)
$`=`$ $`{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\mu =1}{\overset{p}{}}}(\xi _i^{\mu +1}\xi _j^\mu \xi _i^{\mu 1}\xi _j^\mu )x_j(t){\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\mu =1}{\overset{p}{}}}ϵ_{ij}^{\mu 1}\xi _i^{\mu 1}\xi _j^\mu x_j(t).`$ (9)
Using the periodic boundary condition of $`\xi _i^{p+1}=\xi _i^1`$ and $`\xi _i^0=\xi _i^p`$, $`_{\mu =1}^p(\xi _i^{\mu +1}\xi _j^\mu \xi _i^{\mu 1}\xi _j^\mu )=_{\mu =1}^p(\overline{\xi }_i^{\mu +1}\overline{\xi }_j^\mu \overline{\xi }_i^{\mu 1}\overline{\xi }_j^\mu )`$ with $`\overline{\xi }_i^\mu =\xi _i^\mu f`$. Using this relationship, $`u_i(t)`$ is given by

$`u_i(t)`$ $`=`$ $`{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\mu =1}{\overset{p}{}}}(\overline{\xi }_i^{\mu +1}\overline{\xi }_i^{\mu 1})\overline{\xi }_j^\mu x_j(t)`$ (12)
$`{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\mu =1}{\overset{p}{}}}ϵ_{ij}^{\mu 1}\xi _i^{\mu 1}\xi _j^\mu x_j(t)`$
$`=`$ $`{\displaystyle \underset{\mu =1}{\overset{p}{}}}(\overline{\xi }_i^{\mu +1}\overline{\xi }_i^{\mu 1})m^\mu (t)`$
$`{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\mu =1}{\overset{p}{}}}ϵ_{ij}^{\mu 1}\xi _i^{\mu 1}\xi _j^\mu x_j(t)`$
$`=`$ $`(\overline{\xi }_i^{t+1}\overline{\xi }_i^{t1})m^t(t)+{\displaystyle \underset{\mu t}{\overset{p}{}}}(\overline{\xi }_i^{\mu +1}\overline{\xi }_i^{\mu 1})m^\mu (t)`$
$`{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\mu =1}{\overset{p}{}}}ϵ_{ij}^{\mu 1}\xi _i^{\mu 1}\xi _j^\mu x_j(t),`$
where $`m^\mu (t)`$ is an overlap between $`𝝃^\mu (t)`$ and $`𝒙(t)`$ and is given by

$$m^\mu (t)=\frac{1}{Nf(1f)}\underset{i=1}{\overset{N}{}}\overline{\xi }_i^\mu x_i(t).$$
(13)
The first term in equation (12) is a signal term for the retrieval of the target pattern $`𝝃^{t+1}`$. The second and third terms are a cross-talk noise term that represents contributions from non-target patterns other than $`𝝃^{t1}`$ and that prevents $`𝝃^{t+1}`$ from being retrieved. The third term is also a compensation term originated by a deviation from the balance between LTP and LTD of STDP. Since this term is order of $`N`$ with respect to $`N`$, $`O(N)`$, when $`ϵ0`$, it diverges in the thermodynamic $`N`$ limit: $`N\mathrm{}`$. This means that the stored limit cycle using the present learning rule (equation (7)) is unstable in the limit of $`N\mathrm{}`$ when $`ϵ0`$. Therefore, we will discuss the $`ϵ=0`$ case: $`ϵ_{ij}^\mu 𝒩(0,\delta ^2)`$.

We derive the dynamical equations using the method of the statistical neurodynamics Matsumoto and Okada (2002); Okada (1996); Amari and Maginu (1988). When it is possible to store a pattern sequence, a cross-talk noise term, that is, the second and third terms in equation (12) is assumed to obey a Gaussian distribution with the average $`0`$ and time-dependent variance $`\sigma ^2(t)`$ Okada (1996); Matsumoto and Okada (2002). We derive the recursive equations for $`m^t(t)`$ and $`\sigma ^2(t)`$ to investigate whether the memory pattern $`𝝃^t`$ is retrieved or not. Since $`m^t(t)`$ depends on $`\sigma ^2(t)`$, we derive $`\sigma ^2(t)`$. The dynamical equations are derived as

$`m^t(t)`$ $`=`$ $`{\displaystyle \frac{12f}{2}}\text{erf}(\varphi _0){\displaystyle \frac{1f}{2}}\text{erf}(\varphi _1)+{\displaystyle \frac{f}{2}}\text{erf}(\varphi _2),`$ (14)
$`\sigma ^2(t)`$ $`=`$ $`{\displaystyle \underset{a=0}{\overset{t}{}}}{}_{2(a+1)}{}^{}\text{C}_{(a+1)}^{}\alpha q(ta){\displaystyle \underset{b=1}{\overset{a}{}}}U^2(tb+1)`$ (15)
$`+{\displaystyle \frac{\alpha \delta ^2}{(1f)^2}}q(t),`$
$`U(t)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi }\sigma (t1)}}\{(12f+2f^2)e^{\varphi _0^2}`$ (16)
$`+f(1f)(e^{\varphi _1^2}+e^{\varphi _2^2})\},`$
$`q(t)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\{1(12f+2f^2)\text{erf}(\varphi _0)f(1f)(\text{erf}(\varphi _1)`$ (17)
$`+\text{erf}(\varphi _2))\}.`$
where $`\text{erf}(y)=\frac{2}{\sqrt{\pi }}_0^y\mathrm{exp}(u^2)𝑑u`$, and $`\varphi _0=\frac{\theta }{\sqrt{2}\sigma (t1)}`$, $`\varphi _1=\frac{m^{t1}(t1)+\theta }{\sqrt{2}\sigma (t1)}`$, $`\varphi _2=\frac{m^{t1}(t1)+\theta }{\sqrt{2}\sigma (t1)}`$, $`{}_{b}{}^{}\text{C}_{a}^{}=\frac{b!}{a!(ba)!}`$ and $`a!`$ is the factorial with positive integer $`a`$. The detail derivation of the dynamical equations is shown in Appendix A.

## IV Results

Fig.2 shows the dependence of the overlap $`m^t(t)`$ on the loading rate $`\alpha `$ when the mean firing rate of the memory pattern is $`f=0.1`$, and the threshold is $`\theta =0.52`$, which is optimized to maximize the storage capacity. (a) is the case at $`\delta =0.0`$, (b) is at $`\delta =1.0`$, and (c) is at $`\delta =2.0`$.

The lines show analytical results obtained by the dynamical equations (1417). The upper line denotes the steady-state values of the overlap $`m^t(t)`$ in retrieval of the pattern sequence. $`m^t(t)`$ is obtained by setting the initial state of the network at the first memory pattern: $`𝒙(1)=𝝃^1`$. Setting the initial values at $`m^1(1)=1`$, $`\sigma ^2(1)=2\alpha f+\frac{\alpha \delta ^2f}{(1f)^2}`$, $`U(1)=0`$ and $`q(1)=f`$ and using the dynamical equations (1417), $`m^t(t)`$ is obtained. When the overlap at the steady state is smaller than $`0.5`$, the critical loading rate $`\alpha `$ is regarded as the storage capacity $`\alpha _C`$. The storage capacity $`\alpha _C`$ is $`0.27`$(a), $`0.178`$(b), and $`0.087`$(c). The lower line indicates the dependence of an initial critical overlap $`m_C`$ on $`\alpha `$. The stored pattern sequence is retrievable when the initial overlap $`m^1(1)`$ is greater than the critical value $`m_C`$. The region in which $`m^1(1)`$ is larger than $`m_C`$ represents the basin of attraction for the retrieval of the pattern sequence. $`m_C`$ is obtained by setting the initial state of the network at $`𝝃^1`$ with additional noise. We employ the following method to add noise. $`100s\%`$ of the minority components ($`x_i(1)=1`$) are flipped, while the same number of majority components ($`x_i(1)=0`$) are also flipped. The initial overlap $`m^1(1)`$ is given as $`1\frac{2s}{1f}`$. Then the mean firing rate of the network is kept equal to that of the memory pattern, $`f`$. The other initial values are equivalent to the upper line case. When the overlap at the steady state is smaller than $`0.5`$, the initial overlap $`m^1(1)`$ is regarded as the initial critical overlap $`m_C`$. The data points and error bars indicate the results of computer simulations of $`11`$ trials with $`5000`$ neurons: $`N=5000`$. The results are obtained from the equations (17). The data points indicate median values and both ends of the error bars does 1/4 and 3/4 deviations. A discrepancy between the values of $`m_C`$ obtained by the computer simulations and the analytical results is originated from the finite size effect of the computer simulations Matsumoto and Okada (2002).

Fig.3 shows the dependence of the storage capacity $`\alpha _C`$ on the standard derivation $`\delta `$ of the fluctuation of LTD at $`f=0.1`$ and $`\theta =0.52`$.

The solid line shows the analytical results obtained by the same procedure to obtain $`\alpha _C`$ in Fig.2. The data points and error bars show the results of computer simulation of $`10`$ trials at $`N=5000`$. The means and standard deviations of $`\alpha _C`$ of $`10`$ trials are plotted as the data points and the error bars, respectively. As the variance $`\delta ^2`$ increases, $`\alpha _C`$ decreases. In other words, the model is robust against the imbalance between LTP and LTD of STDP. Thus, the balance does not need to be maintained precisely, but must simply be maintained on average.

Fig.4 shows an asymptotic property of $`\alpha _C`$ in a large limit of $`\delta ^2`$.

The solid line shows the analytical results obtained by the same procedure to obtain $`\alpha _C`$ in Fig.2 at $`f=0.1`$ and $`\theta =0.52`$ while the dashed line shows $`\mathrm{log}_{10}\alpha _C=\mathrm{log}_{10}\frac{2}{\pi \delta ^2}1.13`$. This figure indicates that $`\alpha _C`$ converges to $`0`$ as order of $`\frac{1}{\delta ^2}`$ with respect to $`\delta `$, $`O(\frac{1}{\delta ^2})`$, in the large $`\delta ^2`$ limit Mimura et al. .

Next, we discuss the dependency of the basin of attraction on $`\delta `$. Each region between the upper line and the lower line in Fig. 2(a),(b),(c) shows the basin of attraction at $`\delta =0.0,1.0,2.0`$, respectively. The basin becomes smaller as a value of $`\delta `$ increases. However, the initial critical overlap $`m_C`$ is unchanged. To introduce a threshold control scheme is known to enlarge the basin of attraction Okada (1996); Kitano and Aoyagi (1998); Matsumoto and Okada (2002).

## V Discussion

In this paper, we investigated the impact of the imbalance between LTP and LTD of STDP on the retrieval properties of spatio-temporal patterns, employing an associative memory network. We analytically investigated the retrieval properties using the statistical neurodynamics. When the fluctuation of LTD is assumed to obey the Gaussian distribution with mean $`0`$ and variance $`\delta ^2`$, the storage capacity takes a finite value even at large $`\delta `$. This implies that the balance does not need to be maintained precisely, but must be maintained on average. This mechanism might work in the brain. Furthermore, the storage capacity converges to $`0`$ as order $`O(\frac{1}{\delta ^2})`$ in the large $`\delta ^2`$ limit. Finally, we found that the basin of attraction becomes smaller as the fluctuation of LTD increases while the initial critical overlap remains unchanged.

We found that the storage capacity takes a finite value even at large $`\delta `$. When $`\delta `$ is larger than $`1.0`$, LTD sometimes disappears in the learning process. The spatio-temporal patterns seem not to be retrievable. Surprisingly, even in this situation, the patterns are retrievable. This implies that the present model achieves strong robustness against the imbalance between LTP and LTD.

## Appendix A Derivation of dynamical equations by the statistical neurodynamics

The detail derivation of the dynamical equations (1417) is given in this appendix. At first, we give a sketch of the derivation. The main point in this derivation is to divide an internal potential $`u_i(t)`$ into two parts, a signal term for a retrieval of a target pattern and a cross-talk noise term that represents contributions from non-target patterns and prevents the target pattern from being retrieved. We evaluate the cross-talk noise term. Specifically, the internal potential $`u_i(t)`$ of $`i`$-th neuron at time $`t`$ is expressed as (see equation (12))

$`u_i(t)`$ $`=`$ $`{\displaystyle \underset{\mu =1}{\overset{p}{}}}(\overline{\xi }_i^{\mu +1}\overline{\xi }_i^{\mu 1})m^\mu (t)`$ (19)
$`{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\mu =1}{\overset{p}{}}}ϵ_{ij}^{\mu 1}\xi _i^{\mu 1}\xi _j^\mu x_j(t)`$
$`=`$ $`(\overline{\xi }_i^{t+1}\overline{\xi }_i^{t1})m^t(t)+{\displaystyle \underset{\mu t}{\overset{p}{}}}(\overline{\xi }_i^{\mu +1}\overline{\xi }_i^{\mu 1})m^\mu (t)`$
$`{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\mu =1}{\overset{p}{}}}ϵ_{ij}^{\mu 1}\xi _i^{\mu 1}\xi _j^\mu x_j(t)`$
$`=`$ $`(\overline{\xi }_i^{t+1}\overline{\xi }_i^{t1})m^t(t)+z_i(t),`$ (20)
where $`\overline{\xi }_i^\mu =\xi _i^\mu f`$, $`m^\mu (t)`$ is the overlap between $`𝝃^\mu (t)`$ and $`𝒙(t)`$ and is given by

$$m^\mu (t)=\frac{1}{Nf(1f)}\underset{i=1}{\overset{N}{}}\overline{\xi }_i^\mu x_i(t),$$
(21)
and

$`z_i(t)`$ $`=`$ $`{\displaystyle \underset{\mu t}{\overset{p}{}}}(\overline{\xi }_i^{\mu +1}\overline{\xi }_i^{\mu 1})m^\mu (t)`$ (22)
$`{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\mu =1}{\overset{p}{}}}ϵ_{ij}^{\mu 1}\xi _i^{\mu 1}\xi _j^\mu x_j(t).`$
The first term in equation (20) is the signal term and the second term is the cross-talk noise term. Since $`x_i(t)`$ in equation (21) depends on $`\xi _i^\mu `$, the distribution of the cross-talk noise term $`z_i(t)`$ is unknown. The dependence on $`\xi _i^\mu `$ is extracted from $`x_i(t)`$ using the Taylor expansion (see equation (27)). In the thermodynamic limit, $`N\mathrm{}`$, $`m^\mu (t)`$ tends to be deterministic. Therefore, $`x_i^{\{\mu \}}(t)`$, which denotes that it does not include $`\xi _i^\mu `$, is independent of $`\xi _i^\mu `$. This enables us to assume that the cross-talk noise term $`z_i(t)`$ obeys a Gaussian distribution with mean $`0`$ and variance $`\sigma ^2`$ Okada (1996); Matsumoto and Okada (2002). Since the distribution of the cross-talk noise term is known, the recursive equation of the overlap is obtained (see equation (52)).

To extract the dependence on $`\xi _i^\mu `$ from $`x_i(t)`$, the state of the $`i`$-th neuron at time $`t+1`$ is transformed:

$`x_i(t+1)`$ $`=`$ $`F(u_i(t)\theta )`$
$`=`$ $`F({\displaystyle \underset{\nu =1}{\overset{p}{}}}(\overline{\xi }_i^{\nu +1}\overline{\xi }_i^{\nu 1})m^\nu (t)`$
$`{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\nu =1}{\overset{p}{}}}ϵ_{ij}^{\nu 1}\xi _i^{\nu 1}\xi _j^\nu x_j(t)\theta ).`$
The first term in the function $`F()`$ of equation (LABEL:eq.xf) is transformed using the periodic boundary condition of $`\xi _i^{p+1}=\xi _i^1`$ and $`\xi _i^0=\xi _i^p`$:

$`{\displaystyle \underset{\nu =1}{\overset{p}{}}}(\overline{\xi }_i^{\nu +1}\overline{\xi }_i^{\nu 1})m^\nu (t)`$ (24)
$`=`$ $`{\displaystyle \underset{\nu =1}{\overset{p}{}}}\overline{\xi }_i^{\nu +1}m^\nu (t){\displaystyle \underset{\nu =1}{\overset{p}{}}}\overline{\xi }_i^{\nu 1}m^\nu (t)`$
$`=`$ $`{\displaystyle \underset{\nu ^{}=1}{\overset{p}{}}}\overline{\xi }_i^\nu ^{}m^{\nu ^{}1}(t){\displaystyle \underset{\nu ^{}=1}{\overset{p}{}}}\overline{\xi }_i^\nu ^{}m^{\nu +1}(t)`$
$`=`$ $`{\displaystyle \underset{\nu =1}{\overset{p}{}}}\overline{\xi }_i^\nu \{m^{\nu 1}(t)m^{\nu +1}(t)\}.`$
The second term in the function $`F()`$ of equation (LABEL:eq.xf) is transformed using the periodic boundary condition:

$`{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\nu =1}{\overset{p}{}}}ϵ_{ij}^{\nu 1}\xi _i^{\nu 1}\xi _j^\nu x_j(t)`$ (25)
$`=`$ $`{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\nu ^{}=1}{\overset{p}{}}}ϵ_{ij}^\nu ^{}\xi _i^\nu ^{}\xi _j^{\nu ^{}+1}x_j(t).`$
To extract the dependency on $`\xi _i^\mu `$ from $`x_i(t+1)`$, using the equation (24) and (25), $`x(t+1)`$ is divided into two parts, the terms which include $`\xi _i^\mu `$ and the terms which do not include $`\xi _i^\mu `$:

$`x_i(t+1)`$ $`=`$ $`F({\displaystyle \underset{\nu =1}{\overset{p}{}}}\overline{\xi }_i^\nu \{m^{\nu 1}(t)m^{\nu +1}(t)\}`$
$`{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\nu =1}{\overset{p}{}}}ϵ_{ij}^\nu \xi _i^\nu \xi _j^{\nu +1}x_j(t)\theta )`$
$`=`$ $`F(\overline{\xi }_i^\mu \{m^{\mu 1}(t)m^{\mu +1}(t)\}`$
$`{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{j=1}{\overset{N}{}}}ϵ_{ij}^\mu \xi _i^\mu \xi _j^{\mu +1}x_j(t)`$
$`+{\displaystyle \underset{\nu \mu }{\overset{p}{}}}\overline{\xi }_i^\nu \left\{m^{\nu 1}(t)m^{\nu +1}(t)\right\}`$
$`{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\nu \mu }{\overset{p}{}}}ϵ_{ij}^\nu \xi _i^\nu \xi _j^{\nu +1}x_j(t)\theta ).`$
At time $`t`$, the pattern $`𝝃^t`$ is designed to be retrieved. Therefore, we can assume that $`m^t(t)`$ is order of $`1`$ with respect to $`N`$, $`m^t(t)O(1)`$, and $`m^\mu (t)(\mu t)`$ is order of $`1/\sqrt{N}`$ with respect to $`N`$, $`m^\mu (t)O(1/\sqrt{N})`$. Since $`m^\mu (t)O(1/\sqrt{N})`$, $`m^{\mu 1}(t)`$ and $`m^{\mu +1}(t)`$ are order of $`1/\sqrt{N}`$ with respect to $`N`$. Since $`m^{\mu +1}(t)O(1/\sqrt{N})`$ and $`ϵ_{ij}^\mu O(1)`$, the second term in equation (LABEL:eq.taylor1) can be considered to be order of $`1/\sqrt{N}`$ with respect to $`N`$. In the thermodynamic limit, $`N\mathrm{}`$, the first and the second terms in equation (LABEL:eq.taylor1) are small. Using the Taylor expansion up to the first order of $`x_i(t+1)`$, $`x_i(t+1)`$ is transformed:

$$x_i(t+1)=x_i^{(\mu )}(t+1)+u_i^{\{\mu \}}(t)x_i^{(\mu )}(t+1),$$
(27)
where $`x_i^{(\mu )}(t+1)`$ is independent of $`\xi _i^\mu `$, $`x_i^{(\mu )}(t+1)`$ is differential of $`x_i^{(\mu )}(t+1)`$, and

$`u_i^{\{\mu \}}(t)`$ $`=`$ $`\overline{\xi }_i^\mu \left\{m^{\mu 1}(t)m^{\mu +1}(t)\right\}`$ (28)
$``$ $`{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{j=1}{\overset{N}{}}}ϵ_{ij}^\mu \xi _i^\mu \xi _j^{\mu +1}x_j(t),`$
$`x_i^{(\mu )}(t+1)`$ $`=`$ $`F(u_i(t)u_i^{\{\mu \}}(t)\theta ),`$ (29)
$`x_i^{(\mu )}(t+1)`$ $`=`$ $`F^{}(u_i(t)u_i^{\{\mu \}}(t)\theta ).`$ (30)
We assume that the function $`F()`$ is differentiable. This assumption is valid since the average of $`F()`$ over a Gaussian noise term will be taken in a later step (see equation (51)). For $`\mu t`$, the overlap $`m^\mu (t)`$ is expressed as

$`m^\mu (t)`$ $`=`$ $`{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{i=1}{\overset{N}{}}}(\xi _i^\mu f)x_i(t)={\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{i=1}{\overset{N}{}}}\overline{\xi }_i^\mu \left\{x_i^{(\mu )}(t)+u_i^{\{\mu \}}(t1)x_i^{(\mu )}(t)\right\}`$ (32)
$`=`$ $`{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{i=1}{\overset{N}{}}}\overline{\xi }_i^\mu x_i^{(\mu )}(t)+{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{i=1}{\overset{N}{}}}(\overline{\xi }_i^\mu )^2\left\{m^{\mu 1}(t1)m^{\mu +1}(t1)\right\}x_i^{(\mu )}(t)`$
$`\left({\displaystyle \frac{1}{Nf(1f)}}\right)^2{\displaystyle \underset{i=1}{\overset{N}{}}}\overline{\xi }_i^\mu {\displaystyle \underset{j=1}{\overset{N}{}}}ϵ_{ij}^\mu \xi _i^\mu \xi _j^{\mu +1}x_j(t1)x_i^{(\mu )}(t).`$
If the averages over $`\xi _j^\mu (\mu t)`$ and $`ϵ_{ji}^\mu (\mu t)`$ are taken in the right-hand side of equation (32), the third term vanishes since $`\text{E}[ϵ_{ji}^\mu ]=0`$. Since the third term including $`x_j(t1)`$ depends on both $`\xi _j^\mu `$ and $`ϵ_{ji}^\mu `$, $`\xi _j^\mu `$ and $`ϵ_{ji}^\mu `$ are extracted from $`x_j(t)`$ before the averages are taken:

$`x_j(t+1)`$ $`=`$ $`F(\overline{\xi }_j^\mu \{m^{\mu 1}(t)m^{\mu +1}(t)\}{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{ki}{\overset{N}{}}}ϵ_{jk}^\mu \xi _j^\mu \xi _k^{\mu +1}x_k(t){\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{\nu \mu }{\overset{p}{}}}ϵ_{ji}^\nu \xi _j^\nu \xi _i^{\nu +1}x_i(t)`$ (34)
$`{\displaystyle \frac{1}{Nf(1f)}}ϵ_{ji}^\mu \xi _j^\mu \xi _i^{\mu +1}x_i(t)+{\displaystyle \underset{\nu \mu }{\overset{p}{}}}\overline{\xi }_j^\nu \{m^{\nu 1}(t)m^{\nu +1}(t)\}+{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{ki}{\overset{N}{}}}{\displaystyle \underset{\nu \mu }{\overset{p}{}}}ϵ_{jk}^\nu \xi _j^\nu \xi _k^{\nu +1}x_k(t)\theta )`$
$`=`$ $`x_j^{(\mu )(ϵ_{ji}^\mu )}(t+1)+u_j^{\{\mu ,ϵ_{ji}^\mu \}}(t)x_j^{(\mu )(ϵ_{ji}^\mu )}(t+1),`$
where $`x_j^{(\mu )(ϵ_{ji}^\mu )}(t)`$, which is independent of both $`\xi _j^\mu (t)`$ and $`ϵ_{ji}^\mu `$, is the differential of $`x_j^{(\mu )(ϵ_{ji}^\mu )}(t)`$ and

$`u_j^{\{\mu ,ϵ_{ji}^\mu \}}(t)`$ $`=`$ $`\overline{\xi }_j^\mu \left\{m^{\mu 1}(t)m^{\mu +1}(t)\right\}{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{ki}{\overset{N}{}}}ϵ_{jk}^\mu \xi _j^\mu \xi _k^{\mu +1}x_k(t){\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{\nu \mu }{\overset{p}{}}}ϵ_{ji}^\nu \xi _j^\nu \xi _i^{\nu +1}x_i(t)`$ (35)
$`{\displaystyle \frac{1}{Nf(1f)}}ϵ_{ji}^\mu \xi _j^\mu \xi _i^{\mu +1}x_i(t),`$
$`x_j^{(\mu )(ϵ_{ji}^\mu )}(t+1)`$ $`=`$ $`F(u_j(t)u_j^{\{\mu ,ϵ_{ji}^\mu \}}(t)\theta ),`$ (36)
$`x_j^{(\mu )(ϵ_{ji}^\mu )}(t+1)`$ $`=`$ $`F^{}(u_j(t)u_j^{\{\mu ,ϵ_{ji}^\mu \}}(t)\theta ).`$ (37)
The first, second and third term in the right-hand side of equation (35) is order of $`1/\sqrt{N}`$ with respect to $`N`$ and the fourth term is order of $`1/N`$ with respect to $`N`$. In the thermodynamic limit, $`N\mathrm{}`$, $`u_j^{\{\mu ,ϵ_{ji}^\mu \}}(t)`$ is small and equation (LABEL:eq.xmuepsilon) equals equation (34). Substituting equation (34) into the right-hand side of equation (32) and averaging the resultant expressions over $`\xi _j^\mu `$ and $`ϵ_{ji}^\mu `$ and using $`\text{E}[ϵ_{ji}^\mu ]=0`$ yield to the following equation for $`m^\mu (t)(\mu t)`$:

$`m^\mu (t)`$ $`=`$ $`{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{i=1}{\overset{N}{}}}\overline{\xi }_i^\mu x_i^{(\mu )}(t)`$
$`+{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{i=1}{\overset{N}{}}}(\overline{\xi }_i^\mu )^2\{m^{\mu 1}(t1)`$
$`m^{\mu +1}(t1)\}x_i^{(\mu )}(t)`$
$`=`$ $`{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{i=1}{\overset{N}{}}}\overline{\xi }_i^\mu x_i^{(\mu )}(t)`$
$`+U(t)\left\{m^{\mu 1}(t1)m^{\mu +1}(t1)\right\},`$
where $`U(t)=\frac{1}{N}_{i=1}^Nx_i^{(\mu )}(t)`$. Substituting equation (LABEL:eq.mnott2) into equation (19) yields

$`u_i(t)`$ $`=`$ $`(\overline{\xi }_i^{t+1}\overline{\xi }_i^{t1})m^t(t)`$ (40)
$`+{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\mu t}{\overset{p}{}}}(\overline{\xi }_i^{\mu +1}\overline{\xi }_i^{\mu 1})\overline{\xi }_j^\mu x_j^{(\mu )}(t)`$
$`+{\displaystyle \underset{\mu t}{\overset{p}{}}}U(t)\{\overline{\xi }_i^{\mu +1}m^{\mu 1}(t1)`$
$`2\overline{\xi }_i^{\mu 1}m^{\mu 1}(t1)+\overline{\xi }_i^{\mu 1}m^{\mu +1}(t1)\}`$
$`{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\mu =1}{\overset{p}{}}}ϵ_{ij}^{\mu 1}\xi _i^{\mu 1}\xi _j^\mu x_j(t).`$
Substituting equation (34) into the last term of equation (40) yields the following expression for $`u_i(t)`$:

$`u_i(t)`$ $`=`$ $`(\overline{\xi }_i^{t+1}\overline{\xi }_i^{t1})m^t(t)+{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\mu t}{\overset{p}{}}}(\overline{\xi }_i^{\mu +1}\overline{\xi }_i^{\mu 1})\overline{\xi }_j^\mu x_j^{(\mu )}(t)`$ (42)
$`+{\displaystyle \underset{\mu t}{\overset{p}{}}}U(t)\left\{\overline{\xi }_i^{\mu +1}m^{\mu 1}(t1)2\overline{\xi }_i^{\mu 1}m^{\mu 1}(t1)+\overline{\xi }_i^{\mu 1}m^{\mu +1}(t1)\right\}`$
$`{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\mu =1}{\overset{p}{}}}ϵ_{ij}^{\mu 1}\xi _i^{\mu 1}\xi _j^\mu \left\{x_j^{(\mu )(ϵ_{ji}^\mu )}(t)+u_j^{\{\mu ,ϵ_{ji}^\mu \}}(t1)x_j^{(\mu )(ϵ_{ji}^\mu )}(t)\right\}`$
$`=`$ $`(\overline{\xi }_i^{t+1}\overline{\xi }_i^{t1})m^t(t)+{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\mu t}{\overset{p}{}}}(\overline{\xi }_i^{\mu +1}\overline{\xi }_i^{\mu 1})\overline{\xi }_j^\mu x_j^{(\mu )}(t)`$
$`+{\displaystyle \underset{\mu t}{\overset{p}{}}}U(t)\left\{\overline{\xi }_i^{\mu +1}m^{\mu 1}(t1)2\overline{\xi }_i^{\mu 1}m^{\mu 1}(t1)+\overline{\xi }_i^{\mu 1}m^{\mu +1}(t1)\right\}`$
$`{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\mu =1}{\overset{p}{}}}ϵ_{ij}^{\mu 1}\xi _i^{\mu 1}\xi _j^\mu x_j^{(\mu )(ϵ_{ji}^\mu )}(t)`$
$`{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\mu =1}{\overset{p}{}}}ϵ_{ij}^{\mu 1}\xi _i^{\mu 1}\xi _j^\mu \overline{\xi }_j^\mu \left\{m^{\mu 1}(t1)m^{\mu +1}(t1)\right\}x_j^{(\mu )(ϵ_{ji}^\mu )}(t)`$
$`+\left({\displaystyle \frac{1}{Nf(1f)}}\right)^2{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\mu =1}{\overset{p}{}}}ϵ_{ij}^{\mu 1}\xi _i^{\mu 1}\xi _j^\mu {\displaystyle \underset{ki}{\overset{N}{}}}ϵ_{jk}^\mu \xi _j^\mu \xi _k^{\mu +1}x_k(t1)x_j^{(\mu )(ϵ_{ji}^\mu )}(t)`$
$`+\left({\displaystyle \frac{1}{Nf(1f)}}\right)^2{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\mu =1}{\overset{p}{}}}ϵ_{ij}^{\mu 1}\xi _i^{\mu 1}\xi _j^\mu {\displaystyle \underset{\nu \mu }{\overset{p}{}}}ϵ_{ji}^\nu \xi _j^\nu \xi _i^{\nu +1}x_i(t1)x_j^{(\mu )(ϵ_{ji}^\mu )}(t)`$
$`+\left({\displaystyle \frac{1}{Nf(1f)}}\right)^2{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\mu =1}{\overset{p}{}}}ϵ_{ij}^{\mu 1}\xi _i^{\mu 1}\xi _j^\mu ϵ_{ji}^\mu \xi _j^\mu \xi _i^{\mu +1}x_i(t1)x_j^{(\mu )(ϵ_{ji}^\mu )}(t).`$
The fifth, sixth, seventh and eighth terms vanish since $`\text{E}[ϵ_{ij}^\mu ]=0`$, and this yields

$`u_i(t)`$ $`=`$ $`(\overline{\xi }_i^{t+1}\overline{\xi }_i^{t1})m^t(t)+{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\mu t}{\overset{p}{}}}(\overline{\xi }_i^{\mu +1}\overline{\xi }_i^{\mu 1})\overline{\xi }_j^\mu x_j^{(\mu )}(t)`$ (43)
$`+{\displaystyle \underset{\mu t}{\overset{p}{}}}U(t)\left\{\overline{\xi }_i^{\mu +1}m^{\mu 1}(t1)2\overline{\xi }_i^{\mu 1}m^{\mu 1}(t1)+\overline{\xi }_i^{\mu 1}m^{\mu +1}(t1)\right\}`$
$`{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\mu =1}{\overset{p}{}}}ϵ_{ij}^{\mu 1}\xi _i^{\mu 1}\xi _j^\mu x_j^{(\mu )(ϵ_{ji}^\mu )}(t)`$
$`=`$ $`(\overline{\xi }_i^{t+1}\overline{\xi }_i^{t1})m^t(t)+z_i(t),`$ (44)
where

$`z_i(t)`$ $`=`$ $`{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\mu t}{\overset{p}{}}}(\overline{\xi }_i^{\mu +1}\overline{\xi }_i^{\mu 1})\overline{\xi }_j^\mu x_j^{(\mu )}(t)+{\displaystyle \underset{\mu t}{\overset{p}{}}}U(t)\{\overline{\xi }_i^{\mu +1}m^{\mu 1}(t1)2\overline{\xi }_i^{\mu 1}m^{\mu 1}(t1)`$ (45)
$`+\overline{\xi }_i^{\mu 1}m^{\mu +1}(t1)\}{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle }_{j=1}^N{\displaystyle }_{\mu =1}^pϵ_{ij}^{\mu 1}\xi _i^{\mu 1}\xi _j^\mu x_j^{(\mu )(ϵ_{ji}^\mu )}(t).`$
$`z_i(t)`$ is the cross-talk noise term. We assume that the cross-talk noise obeys a Gaussian distribution with mean $`0`$ and time-dependent variance $`\sigma ^2(t)`$: $`\text{E}[z_i(t)]=0,\text{E}[(z_i(t))^2]=\sigma ^2(t)`$ Okada (1996); Matsumoto and Okada (2002). The first and second term of $`z_i(t)`$ are the same cross-talk noise term as that of our previous work Matsumoto and Okada (2002). The square of $`z_i(t)`$ is given by

$`\left\{z_i(t)\right\}^2`$ $`=`$ $`\left({\displaystyle \frac{1}{Nf(1f)}}\right)^2{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\mu t}{\overset{p}{}}}(\overline{\xi }_i^{\mu +1}\overline{\xi }_i^{\mu 1})^2(\overline{\xi }_j^\mu )^2\left\{x_j^{(\mu )}(t)\right\}^2`$ (47)
$`+{\displaystyle \underset{\nu t}{\overset{p}{}}}U(t)^2\left\{\overline{\xi }_i^{\nu +1}m^{\nu 1}(t1)2\overline{\xi }_i^{\nu 1}m^{\nu 1}(t1)+\overline{\xi }_i^{\nu 1}m^{\nu +1}(t1)\right\}^2`$
$`+\left({\displaystyle \frac{1}{Nf(1f)}}\right)^2{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\mu =1}{\overset{p}{}}}(ϵ_{ij}^{\mu 1})^2(\xi _i^{\mu 1})^2(\xi _j^\mu )^2\left\{x_j^{(\mu )(ϵ_{ji}^\mu )}(t)\right\}^2`$
$`+{\displaystyle \frac{2}{Nf(1f)}}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\mu t}{\overset{p}{}}}(\overline{\xi }_i^{\mu +1}\overline{\xi }_i^{\mu 1})\overline{\xi }_j^\mu x_j^{(\mu )}(t)`$
$`\times {\displaystyle \underset{\mu ^{}t}{\overset{p}{}}}U(t)\{\overline{\xi }_i^{\mu ^{}+1}m^{\mu ^{}1}(t1)2\overline{\xi }_i^{\mu ^{}1}m^{\mu ^{}1}(t1)+\overline{\xi }_i^{\mu ^{}1}m^{\mu ^{}+1}(t1)\}`$
$`+2{\displaystyle \underset{\mu t}{\overset{p}{}}}U(t)\left\{\overline{\xi }_i^{\mu +1}m^{\mu 1}(t1)2\overline{\xi }_i^{\mu 1}m^{\mu 1}(t1)+\overline{\xi }_i^{\mu 1}m^{\mu +1}(t1)\right\}`$
$`\times {\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\mu =1}{\overset{p}{}}}ϵ_{ij}^{\mu 1}\xi _i^{\mu 1}\xi _j^\mu x_j^{(\mu )(ϵ_{ji}^\mu )}(t)`$
$`+2\left({\displaystyle \frac{1}{Nf(1f)}}\right)^2{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\mu ^{}=1}{\overset{p}{}}}ϵ_{ij}^{\mu ^{}1}\xi _i^{\mu ^{}1}\xi _j^\mu ^{}x_j^{(\mu ^{})(ϵ_{ji}^\mu ^{})}(t){\displaystyle \underset{k=1}{\overset{N}{}}}{\displaystyle \underset{\mu ^{}t}{\overset{p}{}}}(\overline{\xi }_i^{\mu ^{}+1}\overline{\xi }_i^{\mu ^{}1})\overline{\xi }_k^\mu ^{}x_k^{(\mu ^{})}(t)`$
$`=`$ $`\left({\displaystyle \frac{1}{Nf(1f)}}\right)^2{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\mu t}{\overset{p}{}}}(\overline{\xi }_i^{\mu +1}\overline{\xi }_i^{\mu 1})^2(\overline{\xi }_j^\mu )^2\left\{x_j^{(\mu )}(t)\right\}^2`$
$`+{\displaystyle \underset{\nu t}{\overset{p}{}}}U(t)^2\left\{\overline{\xi }_i^{\nu +1}m^{\nu 1}(t1)2\overline{\xi }_i^{\nu 1}m^{\nu 1}(t1)+\overline{\xi }_i^{\nu 1}m^{\nu +1}(t1)\right\}^2`$
$`+\left({\displaystyle \frac{1}{Nf(1f)}}\right)^2{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{\mu =1}{\overset{p}{}}}(ϵ_{ij}^{\mu 1})^2(\xi _i^{\mu 1})^2(\xi _j^\mu )^2\left\{x_j^{(\mu )(ϵ_{ji}^\mu )}(t)\right\}^2`$
$`=`$ $`{\displaystyle \underset{a=0}{\overset{t}{}}}{}_{2(a+1)}{}^{}\text{C}_{(a+1)}^{}\alpha q(ta){\displaystyle \underset{b=1}{\overset{a}{}}}U^2(tb+1)+{\displaystyle \frac{\alpha \delta ^2}{(1f)^2}}q(t),`$ (48)
where $`p=\alpha N`$, $`q(t)=\frac{1}{N}_{i=1}^N\left\{x_i^{(\mu )}(t)\right\}^2`$, $`{}_{b}{}^{}\text{C}_{a}^{}=\frac{b!}{a!(ba)!}`$, $`a!`$ is the factorial with positive integer $`a`$ and $`\text{E}[(ϵ_{ij}^\mu )^2]=\delta ^2`$. Since $`\text{E}[ϵ_{ij}^\mu ]=0`$, the fourth, fifth and sixth terms in equation (47) vanish. We applied the relationship $`_{a=0}^b({}_{b}{}^{}\text{C}_{a}^{})^2={}_{2b}{}^{}\text{C}_{b}^{}`$ in this derivation. Since $`\text{E}[z_i(t)]=0`$, the variance $`\sigma ^2(t)`$ is equal to $`\text{E}[\{z_i(t)\}^2]`$. We then get the recursive equation for $`\sigma ^2(t)`$:

$`\sigma ^2(t)`$ $`=`$ $`{\displaystyle \underset{a=0}{\overset{t}{}}}{}_{2(a+1)}{}^{}\text{C}_{(a+1)}^{}\alpha q(ta){\displaystyle \underset{b=1}{\overset{a}{}}}U^2(tb+1)`$ (49)
$`+{\displaystyle \frac{\alpha \delta ^2}{(1f)^2}}q(t).`$
The overlap between the state $`𝒙(t)`$ and the retrieval pattern $`𝝃^t`$ is given by

$`m^t(t)`$ $`=`$ $`{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{i=1}{\overset{N}{}}}(\xi _i^tf)x_i(t)`$
$`=`$ $`{\displaystyle \frac{1}{Nf(1f)}}{\displaystyle \underset{i=1}{\overset{N}{}}}(\xi _i^tf)`$
$`\times F\left((\xi _i^t\xi _i^{t2})m^{t1}(t1)+z_i(t1)\theta \right).`$
Since $`u_i(t)`$ is independent and identical distribution (i.i.d.), by the law of large numbers, the average over $`i`$ can be replaced by an average over the memory patterns $`\xi ^\mu `$ and the Gaussian noise term $`z𝒩(0,\sigma ^2)`$. Then, the recursive equation for the overlap $`m^t(t)`$ is transformed:

$`m^t(t)`$ $`=`$ $`{\displaystyle \frac{1}{f(1f)}}{\displaystyle \frac{1}{\sqrt{2\pi }\sigma }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}dze^{\frac{z^2}{2\sigma ^2}}(\xi ^tf)`$
$`\times F((\xi ^t\xi ^{t2})m^{t1}(t1)+z\theta )`$
$`=`$ $`{\displaystyle \frac{1}{f(1f)}}{\displaystyle \frac{1}{\sqrt{2\pi }}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}d\stackrel{~}{z}e^{\frac{\stackrel{~}{z}^2}{2}}(\xi ^tf)`$
$`\times F((\xi ^t\xi ^{t2})m^{t1}(t1)+\sigma (t1)\stackrel{~}{z}\theta )`$
(51)
$`=`$ $`{\displaystyle \frac{12f}{2}}\text{erf}(\varphi _0){\displaystyle \frac{1f}{2}}\text{erf}(\varphi _1)+{\displaystyle \frac{f}{2}}\text{erf}(\varphi _2),`$ (52)
where $`\stackrel{~}{z}=z/\sigma `$, $`\text{erf}(y)=\frac{2}{\sqrt{\pi }}_0^y\mathrm{exp}(u^2)𝑑u`$, and $`\varphi _0=\frac{\theta }{\sqrt{2}\sigma (t1)}`$, $`\varphi _1=\frac{m^{t1}(t1)+\theta }{\sqrt{2}\sigma (t1)}`$, $`\varphi _2=\frac{m^{t1}(t1)+\theta }{\sqrt{2}\sigma (t1)}`$, and $``$ denotes an average over the memory pattern $`\xi ^\mu `$. Since $`x_i^{}(t)x_i^{(\mu )}(t)O(\frac{1}{\sqrt{N}})`$ and the thermodynamic limit, $`N\mathrm{}`$, is considered, $`x_i^{(\mu )}(t)=x_i^{}(t)`$. Using this relationship, we derive $`U(t)`$:

$`U(t)`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}x_i^{(\mu )}(t)={\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}x_i^{}(t)`$ (56)
$`=`$ $`{\displaystyle \frac{1}{f(1f)}}{\displaystyle \frac{1}{\sqrt{2\pi }}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑ze^{\frac{z^2}{2}}`$
$`\times F^{}((\xi ^t\xi ^{t2})m^{t1}(t1)`$
$`+\sigma (t1)z\theta )`$
$`=`$ $`{\displaystyle \frac{1}{f(1f)}}{\displaystyle \frac{1}{\sqrt{2\pi }}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑ze^{\frac{z^2}{2}}z`$
$`\times F((\xi ^t\xi ^{t2})m^{t1}(t1)`$
$`+\sigma (t1)z\theta )`$
$`=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi }\sigma (t1)}}\{(12f+2f^2)e^{\varphi _0^2}`$
$`+f(1f)(e^{\varphi _1^2}+e^{\varphi _2^2})\}.`$
Since $`x_i(t)x_i^{(\mu )}(t)O(\frac{1}{\sqrt{N}})`$ and $`N\mathrm{}`$, $`x_i^{(\mu )}(t)=x_i(t)`$. Using this relationship, we derive $`q(t)`$:

$`q(t)`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}\{x_i^{(\mu )}(t)\}^2={\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}\{x_i(t)\}^2`$ (59)
$`=`$ $`{\displaystyle \frac{1}{f(1f)}}{\displaystyle \frac{1}{\sqrt{2\pi }}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑ze^{\frac{z^2}{2}}`$
$`\times F^2((\xi ^t\xi ^{t2})m^{t1}(t1)`$
$`+\sigma (t1)z\theta )`$
$`=`$ $`{\displaystyle \frac{1}{2}}\{1(12f+2f^2)\text{erf}(\varphi _0)f(1f)(\text{erf}(\varphi _1)`$
$`+\text{erf}(\varphi _2))\}.`$