# Multiplexing rhythmic information by spike timing dependent plasticity

Nimrod Sherf<EMAIL_ADDRESS>Physics Department, Ben-Gurion
University of the Negev  Zlotowski Center for Neuroscience, Ben-Gurion
University of the Negev  Maoz Shamir<EMAIL_ADDRESS>Physics Department,
Ben-Gurion University of the Negev  Zlotowski Center for Neuroscience, Ben-
Gurion University of the Negev  Dept. of Physiology and Cell-biology, Faculty
of Health Sciences, Ben-Gurion University of the Negev

###### Abstract

Rhythmic activity has been associated with a wide range of cognitive processes
including the encoding of sensory information, navigation, the transfer of
emotional information and others. Previous studies have shown that spike-
timing-dependent plasticity (STDP) can facilitate the transfer of rhythmic
activity downstream the information processing pathway. However, STDP has also
been known to generate strong winner-take-all like competitions between
subgroups of correlated synaptic inputs. Consequently, one might expect that
STDP would induce strong competition between different rhythmicity channels
thus preventing the multiplexing of information across different frequency
channels. This study explored whether STDP facilitates the multiplexing of
information across multiple frequency channels, and if so, under what
conditions. We investigated the STDP dynamics in the framework of a model
consisting of two competing sub-populations of neurons that synapse in a
feedforward manner onto a single post-synaptic neuron. Each sub-population was
assumed to oscillate in an independent manner and in a different frequency
band. To investigate the STDP dynamics, a mean field Fokker-Planck theory was
developed in the limit of the slow learning rate. Surprisingly, our theory
predicted limited interactions between the different sub-groups. Our analysis
further revealed that the interaction between these channels was mainly
mediated by the shared component of the mean activity. Next, we generalized
these results beyond the simplistic model using numerical simulations. We
found that for a wide range of parameters, the system converged to a solution
in which the post-synaptic neuron responded to both rhythms. Nevertheless, all
the synaptic weights remained dynamic and did not converge to a fixed point.
These findings imply that STDP can support the multiplexing of rhythmic
information, and demonstrate how functionality (multiplexing of information)
can be retained in the face of continuous remodeling of all the synaptic
weights.

Key words

Hebbian learning; synaptic updating; WTA competition; multiplexing:

###### pacs:

Valid PACS appear here

##  I Introduction

Neuronal oscillations have been described and studied for more than a century
coenen2014adolf ; Haas2003 ; steriade1985abolition ; buzsaki1992high ;
gray1994synchronous ; bragin1999high ; burgess2002functional ;
buzsaki2006rhythms ; shamir2009representation ; buzsaki2015editorial ;
ray2015gamma ; bocchio2017synaptic ; Proskovec2018 ; taub2018oscillations .
Rhythmic activity in the central nervous system has been associated with:
attention, learning, encoding of external stimuli, consolidation of memory and
motor output engel1992temporal ; singer1995visual ; engel2001dynamic ;
fries2005mechanism ; buzsaki2006rhythms ; jensen2007human ;
knyazev2007motivation ; bocchio2017synaptic ; hobson2002cognitive ;
shamir2009representation ; engel2010beta ; buzsaki2015editorial ;
taub2018oscillations ; storchi2017modulation . For example, Bocchio et al.
bocchio2017synaptic  described how situations of perceived threat or perceived
safety are translated into rhythmic activity in different frequency bands in
different brain regions. The transfer of oscillatory signal between large
neuronal populations is not trivial and requires some mechanism to prevent
destructive interference and maintain the rhythmic component.

Recently, it was suggested that synaptic plasticity, and especially spike-
timing-dependent plasticity (STDP), can provide such a mechanism. STDP can be
thought of as a generalization of Hebb’s rule that neurons that “fire together
wire together” Hebb  to the temporal domain. In STDP, the amount of
potentiation (increase in the synaptic strength) and depression (decrease in
the synaptic strength) depends on the temporal relation between the pre- and
post-synaptic firing111 Neuronal interaction is not symmetric. A synaptic
terminal connects two neurons: a transmitting neuron on the pre-synaptic side
and a receiving neuron on the post-synaptic side. When a spike that propagates
along the axon of the pre-synaptic neuron reaches the synaptic terminal it
induces a change in the membrane potential of the post-synaptic neuron.. Luz
and Shamir luz2016oscillations  demonstrated that STDP can facilitate the
transfer of rhythmic activity downstream, and analyzed the features of the
STDP that govern this ability.

In their study, Luz and Shamir focused on the transfer of rhythmic activity
from a single population, thus ignoring the possibility of converging
information from different brain regions. However, it has been suggested that
rhythmic activity in different frequency bands is used for multiplexing
information CFCExample1 ; CFCExample2 ; CFCExample3 ; CFCExample4 ;
saleem2017subcortical . On the other hand, STDP has been shown to generate a
winner-take-all like competition between subgroups of correlated pools of
neurons gilson2009emergence ; gutig2003learning ; morrison2008phenomenological
; song2000competitive ; shamir2019theories . This raises the question whether
STDP can enable the transmission of rhythmic activity in different frequency
bands and facilitate multiplexing of information.

We address this question here in a modelling study. In sections II.1, II.2,
II.3 and II.4, we define the network model and the STDP learning rule. Next,
in section II.5, we derive a mean-field approximation for the STDP dynamics in
the limit of a slow learning rate. Analysis of the STDP dynamics revealed
constraints on the transmission of multi-frequency oscillations. In section
II.6, we extend our results using numerical simulations beyond the analytical
toy model. Finally, in section III, we summarize our results and discuss their
implications.

##  II Results

###  II.1 The pre-synaptic populations

We model a system of two excitatory populations of $N$ neurons each, synapsing
in a feed-forward manner onto the same downstream neuron. Both populations are
assumed to exhibit rhythmic activity in a different frequency, representing a
different feature of the external stimulus. Each neuron is further
characterized by a preferred phase, to which it is more likely to fire during
the cycle. The preferred phases are assumed to be evenly spaced on the ring.
Thus, it is convenient to think of the neurons in each population as organized
on a ring according to their preferred phases of firing, fig. 1.

The spiking activity of neuron $k\in\left\\{1,...N\right\\}$ in population
$\eta\in\left\\{1,2\right\\}$,
$\rho_{\eta,k}(t)=\sum_{i}\delta(t-t_{k,i}^{\eta})$ is a doubly stochastic
process, where $\\{t_{k,i}^{\eta}\\}_{i=1}^{\infty}$ are the spike times.
Given the ‘intensity’ $D_{\eta}$ of feature $\eta\in\left\\{1,2\right\\}$ of
the stimulus, $\rho_{\eta,k}(t)$, follows an independent inhomogeneous Poisson
process statistics with a mean rate that is given by:

$\langle\rho_{\eta,k}(t)\rangle=D_{\eta}(1+\gamma\cos[\nu_{\eta}t-\phi_{\eta,k}]),\
\phi_{\eta,k}=2\pi k/N.$ (1)

where $\nu_{\eta}$ is the angular frequency of oscillations for neurons in
population $\eta$, $\gamma$ is the modulation to the mean ratio of the firing
rate, and $\phi_{\eta,k}$ is the phase of the $k$th neuron from population
$\eta$. As the intensity parameters, $D_{\eta}$, represent features of the
external stimulus they also fluctuate on a timescale which is typically longer
than the characteristic timescale of the neural response. For simplicity we
assumed that $D_{1}$ and $D_{2}$ are independent random variables with an
identical distribution:

$\displaystyle\langle D_{\eta}\rangle=D,$ (2a) $\displaystyle\langle
D_{\eta}D_{\xi}\rangle=D^{2}(1+\sigma^{2}\delta_{\eta\xi}),$ (2b)

where $\langle...\rangle$ denotes averaging with respect to the neuronal noise
and stimulus statistics. The essence of multiplexing is to enable the
transmission of different information channels; hence, the assumption of
independence of $D_{1}$ and $D_{2}$ represents fluctuations of different
stimulus features. This assumption also drives the winner-take-all competition
between the two populations.

###  II.2 Post-synaptic neuron model

Spike time correlations are the driving force of STDP dynamics luz2014effect .
Correlated pairs are more likely to affect the firing of the post-synaptic
neuron, and as a result, to modify their synaptic connection gutig2003learning
; shamir2019theories . To analyze the STDP dynamics we need a simplified model
for the post-synaptic firing that will enable us to compute the pre-post
cross-correlations, and in particular, their dependence on the synaptic
weights.

Following the work of luz2016oscillations ; morrison2008 ; song2000competitive
; kempter2001intrinsic ; kempter1999hebbian ; luz2012balancing , the post-
synaptic neuron is modeled as a linear Poisson neuron with a characteristic
delay $d>0$. The mean firing rate of the post-synaptic neuron at time $t$,
$r_{post}(t)$,is given by

$r_{post}(t)=\frac{1}{N}\sum_{\eta=1}^{2}\sum_{k=1}^{N}w_{\eta,k}\rho_{\eta,k}(t-d),$
(3)

where $w_{\eta,k}$ is the synaptic weight of the $k$th neuron of population
$\eta$.

###  II.3 Temporal correlations & order parameters

The utility of the linear neuron model is that the pre-post correlations are
given as a linear combination of the correlations of the pre-synaptic
populations. The cross-correlation between pre-synaptic neurons is given by:

$\begin{split}\Gamma_{(\eta,j),(\xi,k)}(\Delta
t)&=\langle\rho_{\eta,j}(t)\rho_{\xi,k}(t+\Delta t)\rangle\\\
&=\delta_{\eta\xi}\bigg{(}D^{2}(1+\sigma^{2})\big{(}1+\frac{\gamma^{2}}{2}\cos[\nu_{\eta}\Delta
t\\\ &+\phi_{\eta,j}-\phi_{\xi,k}]\big{)}+\delta_{jk}D\delta(\Delta
t)\bigg{)}\\\ &+D^{2}(1-\delta_{\eta\xi}).\end{split}$ (4)

The correlation between the $j$th neuron from population 1 and the post-
synaptic neuron can therefore be written as

$\begin{split}\Gamma_{(1,j),\text{ post}}(\Delta t)=&\frac{D}{N}\delta(\Delta
t-d)w_{1,j}+\frac{D^{2}}{N}\sum_{k=1}^{N}w_{1,k}(1+\\\
&\sigma^{2})\bigg{(}1+\frac{\gamma^{2}}{2}\cos[\nu_{1}(\Delta
t-d)+\phi_{1,j}-\\\
&\phi_{1,k}]\bigg{)}+\frac{D^{2}}{N}\sum_{l=1}^{N}w_{2,l}.\end{split}$ (5)

Thus, the correlations are determined by the global order parameters $\bar{w}$
and $\tilde{w}e^{i\psi}$, where $\bar{w}$ is the mean synaptic weight and
$\tilde{w}e^{i\psi}$ is its first Fourier component. For $N\gg 1$ these
parameters are defined as follows

$\bar{w}_{\eta}(t)=\int_{0}^{2\pi}w_{\eta}(\phi,t)\frac{d\phi}{2\pi}$ (6)

and

$\tilde{w}_{\eta}(t)e^{i\psi_{\eta}}=\int_{0}^{2\pi}w_{\eta}(\phi,t)e^{i\phi}\frac{d\phi}{2\pi}.$
(7)

The phase $\psi_{\eta}$ is determined by the condition that $\tilde{w}_{\eta}$
is real non-negative. Note that the coupling between the two populations is
only expressed through the last term of the correlation function, eq. 5.

###  II.4 The STDP rule

Based on luz2014effect ; gutig2003learning ; luz2012balancing  we model the
synaptic modification $\Delta w$ following either a pre- or post-synaptic
spike as:

$\Delta{w}=\lambda[f_{+}(w)K_{+}(\Delta t)-f_{-}(w)K_{-}(\Delta t)],$ (8)

The STDP rule, eq. 8, is written as a sum of two processes: potentiation (+,
increase in the synaptic weight) and depression (-, decrease). We further
assume a separation of variables by writing each process as a product of a
weight dependent function, $f_{\pm}(w)$, and a temporal kernel,
$K_{\pm}(\Delta t)$. The term $\Delta t=t_{post}-t_{pre}$ is the time
difference between pre- and post-synaptic spiking. Here we assumed, for
simplicity, that all pairs of pre and post spike times contribute additively
to the learning process via eq. 8. Note, however, that the temporal kernels of
the STDP rule, $K_{\pm}(\Delta t)$ have a finite support. Here we normalized
the kernels, $\int K_{\pm}(\Delta t)d\Delta t=1$. The parameter $\lambda$ is
the learning rate. It is assumed that the learning process is slower than the
neuronal spiking activity and the timescale of changes in the external
stimulus. Here, we used the synaptic weight dependent functions of the form of
gutig2003learning :

$\displaystyle f_{+}(w)$ $\displaystyle=(1-w)^{\mu}$ (9) $\displaystyle
f_{-}(w)$ $\displaystyle=\alpha w^{\mu},$ (10)

where $\alpha>0$ is the relative strength of depression and $\mu\in[0,1]$
controls the non-linearity of the learning rule; thus, $f(w)_{\pm}$ ensures
that the synaptic weights are confined to the boundaries $w\in[0,1]$. Gütig
and colleagues gutig2003learning  showed that the relevant parameter regime
for the emergence of a non-trivial structure is $\alpha>1$ and small $\mu$.

Empirical studies reported a large repertoire of temporal kernels for STDP
rules markram1997regulation ; bi1998synaptic ; sjostrom2001rate ;
zhang1998critical ; Abbott2000 ; froemke2006contribution ;
Nishiyama2000CalciumSR ; shouval2002unified ; WOODIN2003807 . Here we focus on
two families of STDP rules:  1. A temporally asymmetric kernel bi1998synaptic
; zhang1998critical ; Abbott2000 ; froemke2006contribution .  2. A temporally
symmetric kernel Nishiyama2000CalciumSR ; Abbott2000 ; shouval2002unified ;
WOODIN2003807 .

For the temporally asymmetric kernel we use the exponential model:

$K_{\pm}(\Delta t)=\frac{e^{\mp\Delta
t/\tau_{\pm}}}{\tau_{\pm}}\Theta(\pm\Delta t),$ (11)

where $\Delta t=t_{post}-t_{pre}$, $\Theta(x)$ is the Heaviside function, and
$\tau_{\pm}$ is the characteristic timescale of the potentiation $(+)$ or
depression $(-)$. We assume that $\tau_{-}>\tau_{+}$ as typically reported.

For the temporally symmetric learning rule we use a difference of Gaussians
model:

$K_{\pm}(\Delta t)=\frac{1}{\tau_{\pm}\sqrt{2\pi}}e^{-\frac{1}{2}(\frac{\Delta
t}{\tau_{\pm}})^{2}},$ (12)

where $\tau_{\pm}$ is the temporal width. In this case, the order of firing is
not important; only the absolute time difference. We further assume, in both
models, that $\tau_{+}<\tau_{-}$, as is typically reported.

Figure 1: A schematic description of the network architecture showing two pre-
synaptic populations, each oscillating at a different frequency. The output of
these pre-synaptic neurons, serves as a feed-forward input to a single post-
synaptic neuron.

###  II.5 STDP dynamics in the limit of slow learning

Due to noisy neuronal activities, the learning dynamics is stochastic .
However, in the limit of a slow learning rate, $\lambda\rightarrow 0$, the
fluctuations become negligible and one can obtain deterministic dynamic
equations for the (mean) synaptic weights (see luz2014effect  for a detailed
derivation)

$\frac{\dot{w}_{\eta,j}(t)}{\lambda}=I^{+}_{(\eta,j)}(t)-I^{-}_{(\eta,j)}(t)$
(13)

where $\eta=1,2$ and

$I^{\pm}_{(\eta,j)}(t)=f_{\pm}(w_{\eta,j}(t))\int_{-\infty}^{\infty}\Gamma_{(\eta,j),\text{
post}}(\Delta)K_{\pm}(\Delta)d\Delta.$ (14)

Using the correlation structure, eq. 5, eq. 14 yields

$\begin{split}&I^{\pm}_{(\eta,j)}(t)\equiv
f_{\pm}(w_{\eta,j}(t))D^{2}\bar{K}_{\pm}\bigg{(}\bar{w}_{\eta}(t)(1+\sigma^{2})+\bar{w}_{\xi}+\frac{\gamma^{2}}{2}\\\
&(1+\sigma^{2})\frac{\tilde{K}_{\pm}}{\bar{K}_{\pm}}\tilde{w}_{\eta}(t)\cos[\phi_{\eta,j}-\Omega^{\eta}_{\pm}-\nu_{\eta}d-\psi_{\eta}]+\frac{1}{ND\bar{K}_{\pm}}\\\
&K_{\pm}(d)w_{\eta,j}\bigg{)},\end{split}$ (15)

where $\bar{K}_{\pm}$ and $\tilde{K}_{\pm}e^{i\Omega^{\eta}_{\pm}}$ are the
Fourier transforms of the STDP kernels

$\displaystyle\bar{K}_{\pm}=\int_{-\infty}^{\infty}K_{\pm}(\Delta)d\Delta,$
(16)
$\displaystyle\tilde{K}_{\pm}e^{i\Omega^{\eta}_{\pm}}=\int_{-\infty}^{\infty}K_{\pm}(\Delta)e^{-i\nu_{\eta}\Delta}d\Delta.$
(17)

Note that for our specific choice of kernels, $\bar{K}_{\pm}=1$, by
construction.

The dynamics of the synaptic weights can be written in terms of the order
parameters, $\bar{w}$ and $\tilde{w}$ (see eqs. 6 and 7). In the continuum
limit, eq. 13 becomes

$\begin{split}\frac{\dot{w}_{\eta}(\phi,t)}{\lambda}=&F_{\eta,d}(\phi,t)+\bar{w}_{\eta}(t)F_{\eta,0}(\phi,t)(1+\sigma^{2})+\\\
&\tilde{w}_{\eta}(t)F_{\eta,1}(\phi,t)+\bar{w}_{\xi}(t)F_{\eta,0}(\phi,t),\end{split}$
(18)

where

$\displaystyle\begin{split}F_{\eta,d}(\phi,t)=&w_{\eta}(\phi,t)\frac{D}{N}\bigg{(}f_{+}(w_{\eta}(\phi,t))K_{+}(d)-\\\
&f_{-}(w_{\eta}(\phi,t))K_{-}(d)\bigg{)},\end{split}$ (19a) $\displaystyle
F_{\eta,0}(\phi,t)=$ $\displaystyle
D^{2}\bigg{(}\bar{K}_{+}f_{+}(w_{\eta}(\phi,t))-\bar{K}_{-}f_{-}(w_{\eta}(\phi,t))\bigg{)},$
(19b)
$\displaystyle\begin{split}F_{\eta,1}(\phi,t)=&\frac{\gamma^{2}}{2}(1+\sigma^{2})\bigg{(}\tilde{K}_{+}f_{+}(w_{\eta}(\phi,t))\cos[\phi-\Omega^{\eta}_{+}-\\\
&\nu_{\eta}d-\psi_{\eta}]-\tilde{K}_{-}f_{-}(w_{\eta}(\phi,t))\cos[\phi-\Omega^{\eta}_{-}\\\
&-\nu_{\eta}d-\psi_{\eta}]\bigg{)}.\end{split}$ (19c)

Integrating eq. 18 over $\phi$ yields the dynamics of the order parameters

$\begin{split}\frac{\dot{\bar{w}}_{\eta}(t)}{\lambda}=&\bar{F}_{\eta,d}(t)+\bar{w}_{\eta}(t)\bar{F}_{\eta,0}(t)(1+\sigma^{2})+\\\
&\tilde{w}_{\eta}(t)\bar{F}_{\eta,1}(t)+\bar{w}_{\xi}(t)\bar{F}_{\eta,0}(t),\end{split}$
(20)
$\begin{split}&\frac{1}{\lambda}\frac{d}{dt}(\tilde{w}_{\eta}(t)e^{i\psi_{\eta}})=\frac{e^{i\psi_{\eta}}}{\lambda}\big{(}\dot{\tilde{w}}_{\eta}(t)+i\tilde{w}_{\eta}(t)\dot{\psi}_{\eta}(t)\big{)}=\\\
&\tilde{F}_{\eta,d}(t)e^{i\Phi_{\eta,d}}+\bar{w}_{\eta}(t)(1+\sigma^{2})\tilde{F}_{\eta,0}(t)e^{i\Phi_{\eta,0}}+\tilde{w}_{\eta}\tilde{F}_{\eta,1}(t)\\\
&e^{i\Phi_{\eta,1}}+\bar{w}_{\xi}(t)\tilde{F}_{\eta,0}(t)e^{i\Phi_{\eta,0}},\end{split}$
(21)

where

$\displaystyle\bar{F}_{\eta,x}(t)=\int_{-\pi}^{\pi}F_{\eta,x}(\phi,t)\frac{d\phi}{2\pi},$
(22a)
$\displaystyle\tilde{F}_{\eta,x}(t)e^{i\Phi_{\eta,x}}=\int_{-\pi}^{\pi}e^{i\phi}F_{\eta,x}(\phi,t)\frac{d\phi}{2\pi},\
x=d,0,1.$ (22b)

Note that in eq. 21, $\tilde{w}_{2}$ does not appear explicitly in the
dynamics of $\tilde{w}_{1}$. This results from the linearity of the post
synaptic neuron model we chose, eq. 3.

Equations 20 and 21 describe high dimensional coupled non-linear dynamics for
the synaptic weights. Studying the development of rhythmic activity is thus
not a trivial task. To this end, we took an indirect approach. Below we
analyze the homogeneous solution, $w_{\eta,i}=w_{\eta}$ for all $i$, in which
the rhythmic activity does not propagate downstream, and investigate its
stability. We are interested in investigating the conditions in which the
homogeneous solution is unstable, and the STDP dynamics can evolve to a
solution that has the capacity to transmit rhythmic information in both
channels downstream.

#### The homogeneous solution.

The symmetry of the STDP dynamics, eq. 18, with respect to rotation guarantees
the existence of a uniform solution where $w_{\eta,j}(t)=w_{\eta}^{*}$
$\forall j\in\left\\{1,...N\right\\}$ and $\tilde{w}_{\eta}(t)=0$ with
$\eta=1,2$.

Solving the fixed point equation for the homogeneous solution yields

$\frac{f_{-}(w^{*})}{f_{+}(w^{*})}=\frac{1+X_{+}}{1+X_{-}}\equiv\alpha_{c},$
(23)

where

$X_{\pm}\equiv\frac{1}{(2+\sigma^{2})ND}K_{\pm}(d)\geq 0.$ (24)

Due to the scaling of $X_{\pm}$ with $N$, $\alpha_{c}$ is not expected to be
far from $1$. From symmetry $w_{1}^{*}=w_{2}^{*}=w^{*}$:

$w^{*}=\Bigg{(}1+\bigg{(}\frac{\alpha}{\alpha_{c}}\bigg{)}^{1/\mu}\Bigg{)}^{-1}.$
(25)

Substituting the homogeneous solution into the post-synaptic firing rate
equation, eq. 3 yields

$\langle\rho_{post}\rangle=2Dw^{*}.$ (26)

Thus, in the homogeneous solution, the post-synaptic neuron will fire at a
constant rate in time and the rhythmic information will not be relayed
downstream.

#### Stability analysis of the homogeneous solution.

Performing standard stability analysis, we consider small fluctuations around
the homogeneous fixed point, $w_{\eta,j}=w^{*}+\delta{w_{\eta,j}}$, and expand
to first order in the fluctuations:

$\delta\bm{\dot{w}}=\lambda D^{2}\bm{M}\delta\bm{w},$ (27)

where $\bm{M}$ is the stability matrix.

Using eq. 15, the fluctuations can be written as

$\delta\dot{w}_{\eta,j}=\delta I_{(\eta,j)}^{+}-\delta I_{(\eta,j)}^{-}.$ (28)

with

$\begin{split}&\delta I_{(\eta,j)}^{\pm}=\frac{\partial
f_{\pm}(w^{*})}{\partial w_{\eta,j}}D^{2}(2+\sigma^{2})(1+X_{\pm})w^{*}\delta
w_{\eta,j}+f_{\pm}D^{2}\bigg{(}(\\\
&1+\sigma^{2})\delta\bar{w}_{\eta}+\delta\bar{w}_{\xi}+\frac{\gamma^{2}}{2}(1+\sigma^{2})\tilde{K}_{\pm}\cos[\phi_{\eta,j}-\Omega_{\pm}^{\eta}-\nu_{\eta}d\\\
&-\psi_{\eta}]\delta\tilde{w}_{\eta}\bigg{)}+\frac{f_{\pm}(w^{*})}{N}DK_{\pm}(d)\delta
w_{\eta,j}.\end{split}$ (29)

Without loss of generality, taking $\eta=1$, $\xi=2$ yields

$\begin{split}&\delta\dot{w}_{1,j}=-\hat{g}_{0}\delta w_{1,j}-\Delta
f(w^{*})(\delta\bar{w}_{1}+\delta\bar{w}_{2})-\sigma^{2}\Delta
f(w^{*})\delta\bar{w}_{1}\\\
&+\gamma^{2}(1+\sigma^{2})\big{(}f_{+}(w^{*})\tilde{K}_{+}(\nu_{1})\cos[\phi_{1,j}-\Omega_{+}^{1}-\nu_{1}d-\psi_{1}]\\\
&-f_{-}(w^{*})\tilde{K}_{-}(\nu_{1})\cos[\phi_{1,j}-\Omega_{-}^{1}-\nu_{1}d-\psi_{1}]\big{)}\delta\tilde{w}_{1}.\end{split}$
(30)

In the homogeneous fixed point, eq. 23:

$\begin{split}\hat{g_{0}}&=(2+\sigma^{2})\bigg{(}\alpha\mu(1+X_{-})\frac{w^{*\mu}}{1-w^{*}}+f_{+}(w^{*})-f_{-}(w^{*})\bigg{)}\\\
&=g_{0}-(2+\sigma^{2})\Delta f(w^{*}),\end{split}$ (31)

where

$\displaystyle
g_{0}\equiv\alpha\mu(2+\sigma^{2})(1+X_{-})\frac{w^{*\mu}}{1-w^{*}}$ (32a)
$\displaystyle\Delta f(w)\equiv f_{-}(w)-f_{+}(w).$ (32b)

Studying eq. 30, the stability analysis around the homogeneous fixed point
yields four prominent eigenvectors. Two are in the subspace of the uniform
directions of the two populations, and two are in directions of the first
Fourier modes of each population. As the uniform modes of fluctuations,
$\delta\bar{\bm{w}}^{\top}=(\delta\bar{w_{1}},\ \delta\bar{w_{2}})$, span an
invariant subspace of the stability matrix, $\bm{M}$, we can study the
restricted matrix, $\bar{\bm{M}}$, defined by:

$\delta\bm{\dot{\bar{w}}}=\lambda D^{2}\bar{\bm{M}}\delta\bm{\bar{w}}.$ (33)

The matrix $\bar{\bm{M}}$ has two eigenvectors,
$\bm{v}_{\text{u}}^{\top}=(1,1)$ and $\bm{v}_{\text{WTA}}^{\top}=(1,-1)$. The
corresponding eigenvalues are

$\displaystyle\bar{\lambda}_{\text{u}}=-g_{0}$ (34a)
$\displaystyle\begin{split}&\lambda_{\text{WTA}}=\bar{\lambda}_{\text{u}}+2\Delta
f(w^{*}).\end{split}$ (34b)

The first eigenvector represents the uniform mode of fluctuations and its
eigenvalue is always negative, fig. 2(a); hence, the homogeneous solution is
always stable with respect to uniform fluctuations. Furthermore,
$\bar{\lambda}_{\text{u}}$ serves as a stabilizing term in other modes of
fluctuations. We distinguish two regimes: $\alpha<\alpha_{c}$ and
$\alpha>\alpha_{c}$. For $\alpha<\alpha_{c}$, $\lim_{\mu\rightarrow
0+}\bar{\lambda}_{\text{u}}=-\infty$, and the uniform solution is expected to
remain stable. For $\alpha>\alpha_{c}$, $\lim_{\mu\rightarrow
0+}\bar{\lambda}_{\text{u}}=0$, and structure may emerge for sufficiently low
values of $\mu$; see legend color map in fig. 2(a).

The second eigenvalue represents a winner-take-all like mode of fluctuations,
in which synapses in one population suppress the other. Clearly, a winner-
take-all competition will prevent multiplexing. For $\alpha>\alpha_{c}$,
$\lim_{\mu\rightarrow 0+}\lambda_{\text{WTA}}=2(\alpha_{c}-1)$. For the
temporally asymmetric learning rule, eq. 11, $X_{-}=0$; consequently
$\alpha_{c}>1$. In the temporally symmetric difference of Gaussians STDP
model, eq. 12, $\alpha_{c}>1$ if and only if $\tau_{+}<\tau_{-}$, which is the
typical case. For $\alpha<\alpha_{c}$ in the limit of small $\mu$ the
divergence of $\bar{\lambda}_{\text{u}}$ stabilizes fluctuations in this mode.
In this case ($\alpha<\alpha_{c}$), $\lambda_{\text{WTA}}$ reaches its maximum
at an intermediate value of $\mu\in(0,1)$, Figure 2(b). For a small range of
$\alpha<\alpha_{c}$, $\alpha\approx\alpha_{c}$ this maximum can be positive,
see fig. 2(c). Figure 2(b) depicts $\lambda_{\text{WTA}}$ as a function of
$\mu$ for different values of $\alpha$ as shown by color. Note that
$\lambda_{\text{WTA}}$ depends on the temporal structure of the STDP rule
solely via the value of $\alpha_{c}$ and $X_{-}$; however, its sign is
independent of $X_{-}$. As can seen in the figure, for $\alpha>\alpha_{c}>1$,
$\lambda_{\text{WTA}}$ is a decreasing function of $\mu$ and the homogeneous
solution loses stability in the competitive winner-take-all direction in the
limit of small $\mu$. Note that these eigenvalues are not completely identical
in the asymmetric and symmetric cases due to the fact that
$X_{-}^{\text{Symmetric}}\neq X_{-}^{\text{Asymmetric}}=0$. However,
$X_{-}^{\text{Symmetric}}\sim 10^{-3}$; hence, both rules are qualitatively
similar.

(a)

(a)

(b)

(b)

(c)

(c)

Figure 2: The eigenvalues of the uniform and winner-take-all modes for the
asymmetric STDP rule, eq. 11. (a) The uniform eigenvalue
$\bar{\lambda}_{\text{u}}$ is shown as a function of $\mu$, for different
values of $\alpha/\alpha_{\text{c}}$ which are depicted by color. (b) The
eigenvalue of the competitive winner-take-all mode, $\lambda_{\text{WTA}}$, is
shown as a function of $\mu$, for different values of
$\alpha/\alpha_{\text{c}}$ as depicted by color. (Inset) Enlarged section of
the figure, showing the eigenvalue corresponding to the choice of parameters
in our simulations below, section II.6 -depicted by a blue star. (c) The
maximal value of $\lambda_{\text{WTA}}$ in the interval $\mu\in(0,1)$ is shown
as a function of $\alpha$. Note that for $\alpha\geq\alpha_{\text{c}}$,
$\lambda^{\text{max}}_{\text{WTA}}$ is obtained on the boundary $\mu=0$.
Unless stated otherwise in the legends, the parameters used in these figures
are $\sigma=0.6$, $D=10\text{spikes}/\text{sec}$, $N=120$,
$\tau_{-}=50\text{msec}$, $\tau_{+}=20\text{msec}$ and $d=10\text{msec}$.

The rhythmic modes are eigenvectors of the stability matrix $\bm{M}$ with
eigenvalues

$\displaystyle\tilde{\lambda}_{\nu_{\eta}}=$
$\displaystyle\bar{\lambda}_{\text{u}}+(2+\sigma^{2})\Delta f(w^{*})+$ (35)
$\displaystyle\gamma^{2}(1+\sigma^{2})f_{+}(w^{*})\tilde{Q}$
$\displaystyle\tilde{Q}=$
$\displaystyle\tilde{K}_{+}(\nu_{\eta})\cos[\Omega_{+}^{\eta}+\nu_{\eta}d]-$
(36)
$\displaystyle\alpha_{c}\tilde{K}_{-}(\nu_{\eta})\cos[\Omega_{-}^{\eta}+\nu_{\eta}d],$

where $\eta=1,2$. The first two terms which do not depend on the temporal
structure of learning rule or on the frequency, contain the stabilizing term
$\bar{\lambda}_{\text{u}}\leq 0$, and their dependence on $\alpha$ and $\mu$
is similar to that of $\lambda_{\text{WTA}}$; compare with LABEL:eq:LambdaWTA.
The last term depends on the real part of the Fourier transform of the delayed
STDP rule at the specific frequency of oscillations, $\nu_{\eta}$. This last
term can destabilize the system in a direction that will enable the
propagation of rhythmic activity downstream while keeping the competitive WTA
mode stable depending on the interplay between the rhythmic activity and the
temporal structure of the STDP rule. For $\tilde{Q}>0$,
$\tilde{\lambda}_{\nu_{\eta}}$ is an increasing function of the modulation to
the mean ratio $\gamma$. If in addition $\alpha>\alpha_{c}>1$ then
$\tilde{\lambda}_{\nu_{\eta}}$ is an increasing function of $\sigma$ as well.

In the low frequency limit, $\lim_{\nu\rightarrow 0}\tilde{Q}=1-\alpha_{c}$,
and depends on the characteristic timescales of $\tau_{+}$, $\tau_{-}$, and
$d$ only via $\alpha_{c}$. For large frequencies
$\lim_{\nu\rightarrow\infty}\tilde{Q}=0$. In this limit the STDP dynamics
loses its sensitivity to the rhythmic activity. Consequently, the resultant
modulation of the synaptic weights profile, $\tilde{w}$, will become
negligible; hence, effectively rhythmic information will not propagate
downstream even if the rhythmic eigenvalue is unstable luz2016oscillations .
Thus, the intermediate frequency region is the most relevant for multiplexing.

In the case of the temporally symmetric kernel, eq. 12, the value of
$\tilde{Q}$ is given by

$\displaystyle\tilde{Q}=\cos[\nu_{\eta}d]\Big{(}e^{-\frac{1}{2}(\nu_{\eta}\tau_{+})^{2}}-\alpha_{c}e^{-\frac{1}{2}(\nu_{\eta}\tau_{-})^{2}}\Big{)},$
(37)

where $\Omega_{\pm}^{\eta}=0$. Figure 3(a) shows the rhythmic eigenvalue,
$\tilde{\lambda}_{\nu}$, as a function of the oscillation frequency,
$\bar{\nu}\equiv\nu/(2\pi)$, for different values of the delay, $d$ as
depicted by color. Since typically, $\tau_{+}<\tau_{-}$, for finite $\nu$,
$\tilde{K}_{+}(\nu)>\tilde{K}_{-}(\nu)$. Consequently, $\tilde{Q}$ will be
dominated by the potentiation term,
$\cos[\nu_{\eta}d]e^{-\frac{1}{2}(\nu_{\eta}\tau_{+})^{2}}$, except for the
very low frequency range of $\nu\lesssim 1/\tau_{-}$. Typical values for the
delay, $d$, are 1-10$ms$, whereas typical values for the characteristic
timescales for the STDP, $\tau_{\pm}$, are tens of $ms$. As a result, the
specific value of the delay, $d$, does not affect the rhythmic eigenvalue
much, and the system becomes unstable in the rhythmic direction for
$1/\tau_{-}\lesssim\nu\lesssim 1/\tau_{+}$.

Increasing the relative strength of the depression, $\alpha$, strengthens the
stabilizing term $\bar{\lambda}_{\text{u}}$; however, $\Delta f(w^{*})$ scales
approximately linearly with $\alpha$ such that the rhythmic eigenvalue is
elevated, causing the frequency range in which $\tilde{\lambda}_{\nu}>0$ to
widen, fig. 3(b). Similarly, for $\alpha>\alpha_{c}>1$, increasing $\mu$
strengthens $\bar{\lambda}_{\text{u}}$ and reduces the frequency range in
which $\tilde{\lambda}_{\nu}>0$, fig. 3(c). Decreasing the characteristic
timescale of potentiation, $\tau_{+}$ increases the frequency region with an
unstable rhythmic eigenvalue; however, when $\tau_{+}$ becomes comparable to
the delay, $d$, the oscillatory nature of $\tilde{\lambda}_{\nu}$ in $\nu$ is
revealed, fig. 3(d).

In the case of the temporally asymmetric kernel, eq. 11, the value of
$\tilde{Q}$ is given by

$\tilde{Q}^{=}\frac{\cos[\Omega_{+}^{\eta}+\nu_{\eta}d]}{\sqrt{1+(\nu_{\eta}\tau_{+})^{2}}}-\alpha_{c}\frac{\cos[\Omega_{-}^{\eta}+\nu_{\eta}d]}{\sqrt{1+(\nu_{\eta}\tau_{-})^{2}}},$
(38)

with $\Omega_{\pm}^{\eta}=\mp\arctan(\nu_{\eta}\tau_{\pm})$.

The main difference between the temporally symmetric and the asymmetric rules
is that due to the discontinuity of the asymmetric STDP kernel,
$\tilde{K}_{\pm}$ decay algebraically rather than exponentially fast with
$\nu$. As a result, the phase $\cos[\Omega_{\pm}^{\eta}+\nu_{\eta}d]$, plays a
more central role in controlling the stability of the rhythmic eigenvalue. As
above, since typically, $\tau_{+}<\tau_{-}$, then
$\tilde{K}_{+}(\nu)>\tilde{K}_{-}(\nu)$. Figure 4(a) shows the rhythmic
eigenvalue, $\tilde{\lambda}_{\nu}$, for different values of $d$. The dashed
vertical lines depict the frequencies at which the potentiation term changes
its sign, $\Omega_{+}^{\eta}+\nu^{*}d=\pi/2$. As can be seen from the figure,
for this choice of parameters the upper cutoff of the central frequency range
in which the rhythmic eigenvalue is unstable is dominated by $\nu^{*}$, which
is governed by the delay.

The effects of parameters $\alpha$ and $\mu$ show similar trends as for the
symmetric STDP rule. Specifically, increasing $\mu$ or decreasing $\alpha$, in
general, shrinks the region in which fluctuations in the rhythmic direction
are unstable, figs. 4(b) and 4(c). Increasing the characteristic timescale of
potentiation, $\tau_{+}$ beyond that of the depression, makes the depression
term, $\tilde{K}_{-}(\nu_{\eta})\cos[\Omega_{-}^{\eta}+\nu_{\eta}d]$, more
dominant, fig. 4(d). In this case the lower frequency cutoff will be dominated
by the change of sign in the depression term; i.e., by the angular frequency
$\nu^{*}$, such that $\Omega_{-}^{\eta}+\nu^{*}d=\pi/2$, shown by the light
brown curve ($\tau_{+}=50\text{msec}$) in fig. 4(d).

From the above analysis, one would expect that multiplexing will emerge via
STDP if both rhythmic eigenvalues are unstable and the winner-take-all
eigenvalue is stable. However, this intuition relies on two major factors. The
first is that the, stability analysis is local and can only describe the
dynamics near the homogeneous fixed point. The second is the choice of a
linear neuron model, eq. 3, which resulted in the lack of explicit interaction
term between the two rhythms $\tilde{w}_{1}$ and $\tilde{w}_{2}$ in eq. 21.
Both of these issues are addressed below by studying the STDP dynamics using a
conductance based Hodgkin-Huxley type numerical model of a spiking post-
synaptic neuron.

(a)

(a)

(b)

(b)

(c)

(c)

(d)

(d)

Figure 3: The rhythmic eigenvalue, $\tilde{\lambda}_{\nu}$, in the case of the
symmetric learning rule, eq. 12. (a) $\tilde{\lambda}_{\nu}$ is shown as a
function of frequency, $\bar{\nu}=\nu/(2\pi)$, for different values of $d$ as
depicted by color. (b) $\tilde{\lambda}_{\nu}$ is shown as a function of
frequency for different values of $\alpha/\alpha_{\text{c}}$ as depicted by
color. (c) $\tilde{\lambda}_{\nu}$ as a function of frequency for different
values of $\mu$ \- by color. (d) $\tilde{\lambda}_{\nu}$ as a function of
frequency for different values of $\tau_{+}$ \- by color, with
$\tau_{-}=50\text{msec}$. The blue (11Hz) and red (14Hz) stars show the
eigenvalue corresponding to the parameters of the simulations, in section
II.6. Unless stated otherwise in the legends, the parameters used to produce
these figures are $\gamma=1$, $\sigma=0.6$, $D=10\text{spikes}/\text{sec}$,
$N=120$, $\tau_{-}=50\text{msec}$, $\tau_{+}=20\text{msec}$, $\mu=0.01$,
$\alpha=1.1$ and $d=10\text{msec}$.

(a)

(a)

(b)

(b)

(c)

(c)

(d)

(d)

Figure 4: The rhythmic eigenvalue, $\tilde{\lambda}_{\nu}$, in the case of the
asymmetric learning rule, eq. 11. (a) $\tilde{\lambda}_{\nu}$ is shown as a
function of the frequency, $\bar{\nu}=\nu/(2\pi)$, for different values of $d$
as depicted by color. (b) $\tilde{\lambda}_{\nu}$ is shown as a function of
the frequency, for different values of $\alpha/\alpha_{\text{c}}$ as depicted
by color. (c) $\tilde{\lambda}_{\nu}$ as a function of the frequency, for
different values of $\mu$ \- by color. (d) $\tilde{\lambda}_{\nu}$ as a
function of frequency, for different values of $\tau_{+}$ \- by color, with
$\tau_{-}=50\text{msec}$. The black (11Hz) and purple (14Hz) stars show the
eigenvalue with the corresponding parameters to the simulations in section
II.6. Unless stated otherwise in the legends, the parameters used to produce
these figures are $\gamma=1$, $\sigma=0.6$, $D=10\text{spikes}/\text{sec}$,
$N=120$, $\tau_{-}=50\text{msec}$, $\tau_{+}=20\text{msec}$, $\mu=0.01$,
$\alpha=1.1$ and $d=10\text{msec}$.

###  II.6 Conductance based neuron model

In our numerical analysis we used the conductance based model described in
Shriki et al. shriki2003rate  for the post-synaptic neuron. This choice was
motivated by the ability to control the linearity of the neuron’s response to
its inputs. Often the response of a neuron is quantified using an f-I curve,
which maps the frequency (f) of the neuronal spiking response to a certain
level of injected current (I). In the Shriki model shriki2003rate , a strong
transient potassium A-current yields a threshold linear f-I curve to a good
approximation.

The synaptic weights in all simulations were updated according to the STDP
rule presented in eq. 8 with all pre-post spike time pairs contribute
additively to the synaptic plasticity. Further details on the numerical
simulations can be found in section IV.

Asymmetric learning rule with a strong A-current. Figure 5 presents the
results of simulating the STDP dynamics with a conductance based post-synaptic
neuron with a strong A-current. In this regime, the f-I curve of the post-
synaptic neuron is well approximated by a threshold-linear function (see fig.
1 in shriki2003rate ). Consequently, it is reasonable to expect that our
analytical results will hold, in the limit of a slow learning rate. Indeed,
even though the initial conditions of all the synapses are uniform, the
uniform solution loses its stability, and a structure that shows phase
preference begins to emerge, figs. 5a and 5b. After about 1000 sec of
simulation time, the STDP dynamics of each sub-population converges to an
approximately periodic solution. The order parameters $\bar{w}$ and
$\tilde{w}$ appear to converge to a fixed point, fig. 5c, while the phases
$\psi_{1}$ and $\psi_{2}$ continue to drift with a relatively fixed velocity,
fig. 5d. For the specific choice of parameters used in this simulation, the
competitive winner-take-all eigenvalue is stable, see inset fig. 2(b), whereas
the rhythmic eigenvalue is unstable, see colored stars in fig. 4(d).

STDP dynamics, asymmetric learning rule

Figure 5: Numerical simulation of STDP dynamics of two populations with 120
excitatory synapses each with the temporally asymmteric STDP rule, eq. 11. (a)
and (b) The synaptic weights are shown as a function of time, for population 1
($11\text{Hz}$) and population 2 ($14\text{Hz}$) in a and b, respectively. The
colors represent the preferred phases of the pre-synaptic neurons on the ring.
(c) The dynamics of the order parameters, the mean, $\bar{w}$, (red and
magenta) and first Fourier component, $\tilde{w}$, (blue and green) are shown
as a function of time for both populations. (d) The dynamics of the phases,
the phase, $\psi_{i}$, of each population $i$ is shown as a function of time.
Note that the drift velocities are approximately constant in time.

Symmetric learning rule with a strong A-current.

A typical result of simulating the STDP dynamics for the symmetric learning
rule is shown in fig. 6. As above, the post-synaptic neuron was characterized
by a strong A-current, yielding a (threshold) linear f-I curve. In this
example, the uniform initial conditions of the synaptic weights evolve to an
approximately limit cycle solution for each population after $\sim 1000$ sec.
As in the asymmetric case, both populations coexist; hence, STDP facilitates
rhythmic information for both learning rules. As before, for the specific
choice of parameters used in this simulation the competitive winner-take-all
eigenvalue is stable and the rhythmic eigenvalue is unstable, see colored
stars in fig. 3(d).

STDP dynamics, symmetric learning rule

Figure 6: Numerical simulation of STDP dynamics of two populations with 120
excitatory synapses, each with the symmetric learning rule, eq. 12. (a) and
(b) The synaptic weights are shown as a function of time, for population 1
($11\text{Hz}$) and population 2 ($14\text{Hz}$) in a and b, respectively. (c)
The dynamics of the order parameters, the mean, $\bar{w}$, (red and magenta)
and first Fourier component, $\tilde{w}$, (blue and green) of each population
are shown as a function of time. (d) The dynamics of the phases, the phase,
$\psi_{i}$, of each population $i$ is shown as a function of time.

Asymmetric learning rule with a non linear neuron. Above we studied the STDP
dynamics with a linear post-synaptic neuron numerically. This choice limits
the interaction between the two populations. However, it is not uncommon for
neurons to have a non-linear f-I curve. In order to determine whether STDP can
facilitate multiplexing despite the expected interaction due to the non-linear
f-I curve, we used the Shriki model shriki2003rate  with no A-current. A
typical simulation result is shown in fig. 7.As can be seen from the figure,
the system converges to a dynamical solution that is qualitatively similar to
the case of the linear neuron model (fig. 5). Specifically, the order
parameters of each synaptic population converge to a solution that enables the
transmission of rhythmic activity downstream, whereas the synaptic weights
themselves in each population converge to a dynamic solution that appears to
be approximately a limit cycle.

Synaptic plasticity Vs time, asymmetric learning rule, non-linear neuron

Figure 7: Numerical simulation of the STDP dynamics of two populations with
120 excitatory synapses each with the asymmetric learning rule, eq. 11 and a
non-linear post-synaptic neuron. (a) and (b). The synaptic weights are shown
as a function of time, for population 1 ($11\text{Hz}$) and population 2
($14\text{Hz}$) in a and b, respectively. (c) The dynamics of the order
parameters- the mean, $\bar{w}$, (red and magenta) and first Fourier
component, $\tilde{w}$, (blue and green) of each population are shown as a
function of time. (d) The dynamics of the phases: the phase, $\psi_{i}$, of
each population $i$ is shown as a function of time.

##  III  Summary & Discussion

Rhythmic activity in the brain has fascinated and puzzled neuroscientists for
more than a century. Nevertheless, the utility of rhythmic activity remains
enigmatic. One explanation frequently put forward is that of the multiplexing
of information. Our work provides some measure of support for this hypothesis
from the theory of unsupervised learning.

We studied the computational implications of a microscopic learning rule,
namely STDP, in the absence of a reward or teacher signal. Previous work shows
that STDP generates a strong winner-take-all like competition between
subgroups of correlated neurons, thus effectively eliminating the possibility
of multiplexing gutig2003learning . Our work demonstrates that rhythmic
activity does not necessarily generate competition between different rhythmic
channels. Moreover, we found that under a wide range of parameters STDP
dynamics will develop spontaneously the capacity for multiplexing.

Not every learning rule; i.e., a choice of parameters that describes the STDP
update rule, will support multiplexing. This observation provides a natural
test for our theory. Clearly, if multiplexing has evolved via a process of
STDP, then the STDP rule must exhibit instability with respect to the rhythmic
modes and stability against fluctuations in the winner-take-all direction.
These constraints serve as the basic predictions for our theory.

In Luz & Shamir luz2016oscillations , due to the underlying U(1) symmetry, the
system converged to a limit cycle solution. For a single population, the order
parameter, $\tilde{w}e^{i\psi}$, will drift on the ring $|\mathbf{w}|=const$
with a constant velocity. Here, in the linear neuron model, in the limit of a
slow learning rate one expects the system to converge to the product space of
two limit cycles. As there is no reason to expect that the ratio of drift
velocities of the two populations will be rational, the order parameters
$\tilde{w}_{1}e^{i\psi_{1}}$, $\tilde{w}_{2}e^{i\psi_{2}}$ will, most likely,
cover the torus uniformly. Nevertheless, the reduced dynamics of each
population will exhibit a limit cycle. This intuition relies on the lack of
interaction between $\tilde{w}^{1}$ and $\tilde{w}^{2}$ in the dynamics of the
order parameters, eq. 21. Essentially, the interaction between the two
populations is mediated solely via their mean component, $\bar{w}$. However,
introducing non-linearity to the response properties of the post-synaptic
neuron will induce an interaction between the modes. Similarly, for any finite
learning rate, $\lambda\neq 0$, the rhythmic modes will not be orthogonal and
consequently will be correlated.

Introducing a post-synaptic neuron with a non-linear f-I curve and a finite
learning rate induces an interaction between the two populations.
Consequently, for a finite learning rate the order parameters
$\tilde{w}_{1}e^{i\psi_{1}}$, $\tilde{w}_{2}e^{i\psi_{2}}$ will not be
confined to a torus and $\tilde{w}_{\eta}e^{i\psi_{\eta}}$ will fluctuate
around (in contrast with on) the ring. Traces for this behaviour can be seen
by the fact that drift velocity in the numerical simulations is not constant
and the global order parameters exhibit small fluctuations around, see 7c $\&$
7d for example. Nevertheless, the post-synaptic neuron will respond to both
rhythms; hence, multiplexing will be retained.

Synaptic weights in the central nervous system are highly volatile and
demonstrate high turnover rates as well as considerable size changes that
correlate with the synaptic weight mongillo2017intrinsic ;
attardo2015impermanence ; grutzendler2002long ; zuo2005development ;
holtmaat2005transient ; holtmaat2006experience ; loewenstein2011multiplicative
; loewenstein2015predicting . How can the brain retain functionality in the
face of these considerable changes in synaptic connectivity? Our work
demonstrates how functionality, in terms of retaining the ability to transmit
downstream rhythmic information in several channels, can be retained even when
the entire synaptic population is modified throughout its entire dynamic
range. Here, robustness of function is ensured by stable dynamics for the
global order parameters.

##  IV Appendix

#### The conductance based model for the Post-synaptic neuron

We used the conductance based model of Shriki et al. shriki2003rate . The
model is fully defined and its response to synchronous inputs is discussed in
shriki2003rate .

Here, we introduce the main features and equations utilized in this model. The
voltage equation is

$C_{m}\dot{V}=-I_{L}-I_{\text{Na}}-I_{k}-I_{A}+\sum_{i}g_{i}(E_{i}-V),$ (39)

where $I_{L}$, the leak current is given by $I_{L}=g_{L}(V-E_{L})$. $I_{Na}$
and $I_{k}$ are the sodium and potassium currents respectively and are given
by $I_{Na}=\bar{g}_{Na}m_{\infty}^{3}h(V-E_{Na})$ and
$I_{k}=\bar{g}_{k}n^{4}h(V-E_{k})$. The relaxation equations of the the gating
variables $x=h,n$ are $dx/dt=(x_{\infty}-x)/\tau_{x}$. The time independent
functions $x_{\infty}=h_{\infty},n_{\infty},m_{\infty}$ and $\tau_{x}$ are:
$x_{\infty}=\alpha_{x}/(\alpha_{x}+\beta_{x})$ and
$\tau_{x}=0.1/(\alpha_{x}+\beta_{x})$, with
$\alpha_{m}=-0.1(V+30)/(exp(-0.1(V+30))-1)$, $\beta_{m}=4exp(-(V+55)/18)$,
$\alpha_{h}=0.07exp(-(V+44)/20)$ , $\beta_{h}=1/(exp(-0.1(V+14))+1)$,
$\alpha_{n}=-0.01(V+34)/(exp(-0.1(V+34))-1)$ and
$\beta_{n}=0.125exp(-(V+44)/80)$.

The A-current $I_{A}$ that linearizes the f-I relationship is
$I_{A}=\bar{g}_{A}a_{\infty}^{3}b(V-E_{k})$, where
$a_{\infty}=1/(exp(-(V+50)/20)+1)$ and $db/dt=(b_{\infty}-b)/\tau_{A}$. The
time independent function of $b_{\infty}$ is $b_{\infty}=1/(exp((V+80)/6)+1)$
with the voltage independent time constant $\tau_{A}$. In the case of the non-
linear f-I curve, fig. 7, we took this value to be 0.

The term $g_{i}$ is the total conductance of the $i$th pre-synaptic population
and can be written as follows

$g_{x}(t)=g_{x}^{0}\sum_{j=1}^{N_{x}}(w_{j}^{x}(t)\sum_{s}\frac{[t-t_{j}^{s}]_{+}}{\tau_{x}}e^{-(t-t_{j}^{s})/\tau_{x}}).$
(40)

Here, $N_{x}$ is the number of neurons in the $x$th population, $w_{j}^{x}(t)$
is the synaptic weight of the $j$th neuron from the $x$th population and
$[y]_{+}=max(0,y)$. The $s$ spike of the $j$th neuron is denoted by
$t_{j}^{s}$. We used $g_{x}^{0}=g_{x}^{R}S_{x}$ with $S_{x}=1000/N_{x}$ and
$g_{x}^{R}=90\text{mS/cm}^{2}$ (see gutig2003learning ; luz2016oscillations
for details).

The membrane capacity is $c_{m}=0.1\mu\text{F/cm}^{2}$. The sodium, potassium
and leak conductances are $\bar{g}_{Na}=100\text{mS/cm}^{2}$,
$\bar{g}_{k}=40\text{mS/cm}^{2}$ and $g_{L}=0.05\text{mS/cm}^{2}$
respectively. The conductance and characteristic decay times of the A-current
are $g_{A}=20\text{mS/cm}^{2}$ and $\tau_{A}=20\text{msec}$ respectively. The
reversal potentials of the ionic and synaptic currents are
$E_{Na}=55\text{mV}$, $E_{K}=−80\text{mV}$, $E_{L}=−65\text{mV}$,
$E_{e}=0\text{mV}$, and $E_{in}=−80\text{mV}$.

#### Modeling pre-synaptic activity

During the spiking neuron simulations, pre-synaptic activities were modeled by
independent inhomogeneous Poisson processes, with a modulated mean firing rate
given by eq. 1, where $\gamma=1$, and each second $D$ was independently
sampled from a uniform distribution with a minimum of
$7\,\text{spikes}\,\times\text{sec}^{-1}$ and a maximum of
$13\,\text{spikes}\,\times\text{sec}^{-1}$,
($D=7+\textbf{U}(0,6)\,\text{spikes}\,\times\text{sec}^{-1}$). Each pre-
synaptic neuron, spiked according to an approximated Bernoulli process, with a
probability of $p\approx r\Delta t$, where $r$ is the mean firing rate (eq. 1)
and $\Delta t=1\text{msec}$. The number of pre-synaptic neurons in each
population was $N=120$, all excitatory.

#### STDP

The learning rate of the simulations presented in this paper is
$\lambda=0.01$. The power of the synaptic weights in eq. 9 is $\mu=0.01$ and
the ratio of depression to potentiation is $\alpha=1.1$. In order to update
the synaptic weights we relied on the separation of time scales between the
synaptic dynamics of eq. 39; hence, the synaptic weights were updated every
$1\text{sec}$ of simulation.

  * •

Asymmetric learning rule The characteristic decay times were chosen to be
$\tau_{+}=20\text{ms}$ and $\tau_{-}=50\text{ms}$.

  * •

Symmetric learning rule Here, based on our analysis, we chose the ratio of
decay times to be $\sim 10$, $\tau_{+}=5\text{ms}$, $\tau_{-}=50\text{ms}$.

In addition, the initial conditions for all neurons were uniform; i.e.,
$w_{\eta}(\phi,t=0)=0.5$, $\eta=1,2$.

## ACKNOWLEDGMENTS

This research was supported by the Israel Science Foundation (ISF) grant
number 300/16, and in part by the National Science Foundation under Grant No.
NSF PHY-1748958 and the United States-Israel Binational Science Foundation
grant 2013204.

## References

  * [1] Anton Coenen, Edward Fine, and Oksana Zayachkivska.  Adolf beck: A forgotten pioneer in electroencephalography.  Journal of the History of the Neurosciences, 23(3):276–286, 2014\.
  * [2] L. Haas.  Hans berger (1873-1941), richard caton (1842-1926), and electroencephalography.  J Neurol Neurosurg Psychiatry, 74(1):9–9, Jan 2003.  12486257[pmid].
  * [3] M Steriade, M Deschenes, L Domich, and C Mulle.  Abolition of spindle oscillations in thalamic neurons disconnected from nucleus reticularis thalami.  Journal of neurophysiology, 54(6):1473–1497, 1985.
  * [4] Gyorgy Buzsaki, Zsolt Horvath, Ronald Urioste, Jamille Hetke, and Kensall Wise.  High-frequency network oscillation in the hippocampus.  Science, 256(5059):1025–1027, 1992.
  * [5] Charles M Gray.  Synchronous oscillations in neuronal systems: mechanisms and functions.  Journal of computational neuroscience, 1(1-2):11–38, 1994.
  * [6] Anatol Bragin, Jerome Engel Jr, Charles L Wilson, Itzhak Fried, and Gyorgy Buzsáki.  High-frequency oscillations in human brain.  Hippocampus, 9(2):137–142, 1999.
  * [7] Adrian P Burgess and Lia Ali.  Functional connectivity of gamma eeg activity is modulated at low frequency during conscious recollection.  International Journal of Psychophysiology, 46(2):91–100, 2002.
  * [8] Gyorgy Buzsaki.  Rhythms of the Brain.  Oxford University Press, 2006.
  * [9] Maoz Shamir, Oded Ghitza, Steven Epstein, and Nancy Kopell.  Representation of time-varying stimuli by a network exhibiting oscillations on a faster time scale.  PLoS computational biology, 5(5):e1000370, 2009.
  * [10] György Buzsáki and Walter Freeman.  Editorial overview: brain rhythms and dynamic coordination.  Current opinion in neurobiology, 31:v–ix, 2015.
  * [11] Supratim Ray and John HR Maunsell.  Do gamma oscillations play a role in cerebral cortex?  Trends in cognitive sciences, 19(2):78–85, 2015.
  * [12] Marco Bocchio, Sadegh Nabavi, and Marco Capogna.  Synaptic plasticity, engrams, and network oscillations in amygdala circuits for storage and retrieval of emotional memories.  Neuron, 94(4):731–743, 2017.
  * [13] Amy L. Proskovec, Alex I. Wiesman, Elizabeth Heinrichs-Graham, and Tony W. Wilson.  Beta oscillatory dynamics in the prefrontal and superior temporal cortices predict spatial working memory performance.  Scientific Reports, 8(1):8488, 2018.
  * [14] Aryeh Hai Taub, Rita Perets, Eilat Kahana, and Rony Paz.  Oscillations synchronize amygdala-to-prefrontal primate circuits during aversive learning.  Neuron, 97(2):291–298, 2018.
  * [15] Andreas K Engel, Peter König, Andreas K Kreiter, Thomas B Schillen, and Wolf Singer.  Temporal coding in the visual cortex: new vistas on integration in the nervous system.  Trends in neurosciences, 15(6):218–226, 1992.
  * [16] Wolf Singer and Charles M Gray.  Visual feature integration and the temporal correlation hypothesis.  Annual review of neuroscience, 18(1):555–586, 1995.
  * [17] Andreas K Engel, Pascal Fries, and Wolf Singer.  Dynamic predictions: oscillations and synchrony in top–down processing.  Nature Reviews Neuroscience, 2(10):704, 2001.
  * [18] Pascal Fries.  A mechanism for cognitive dynamics: neuronal communication through neuronal coherence.  Trends in cognitive sciences, 9(10):474–480, 2005.
  * [19] Ole Jensen, Jochen Kaiser, and Jean-Philippe Lachaux.  Human gamma-frequency oscillations associated with attention and memory.  Trends in neurosciences, 30(7):317–324, 2007.
  * [20] Gennady G Knyazev.  Motivation, emotion, and their inhibitory control mirrored in brain oscillations.  Neuroscience & Biobehavioral Reviews, 31(3):377–395, 2007.
  * [21] J Allan Hobson and Edward F Pace-Schott.  The cognitive neuroscience of sleep: neuronal systems, consciousness and learning.  Nature Reviews Neuroscience, 3(9):679, 2002.
  * [22] Andreas K Engel and Pascal Fries.  Beta-band oscillations—signalling the status quo?  Current opinion in neurobiology, 20(2):156–165, 2010.
  * [23] Riccardo Storchi, Robert A Bedford, Franck P Martial, Annette E Allen, Jonathan Wynne, Marcelo A Montemurro, Rasmus S Petersen, and Robert J Lucas.  Modulation of fast narrowband oscillations in the mouse retina and dlgn according to background light intensity.  Neuron, 93(2):299–307, 2017.
  * [24] Donald Olding Hebb.  The organization of behavior: A neuropsychological theory.  Psychology Press, 2005.
  * [25] Yotam Luz and Maoz Shamir.  Oscillations via spike-timing dependent plasticity in a feed-forward model.  PLoS computational biology, 12(4):e1004878, 2016.
  * [26] Takumi Sase, Yuichi Katori, Motomasa Komuro, and Kazuyuki Aihara.  Bifurcation analysis on phase-amplitude cross-frequency coupling in neural networks with dynamic synapses.  Frontiers in Computational Neuroscience, 11:18, 2017.
  * [27] Alexandre Hyafil, Anne-Lise Giraud, Lorenzo Fontolan, and Boris Gutkin.  Neural cross-frequency coupling: connecting architectures, mechanisms, and functions.  Trends in neurosciences, 38(11):725–740, 2015.
  * [28] Anita K Roopun, Mark A Kramer, Lucy M Carracedo, Marcus Kaiser, Ceri H Davies, Roger D Traub, Nancy J Kopell, and Miles A Whittington.  Temporal interactions between cortical rhythms.  Frontiers in neuroscience, 2:34, 2008.
  * [29] Felix Siebenhühner, Sheng H Wang, J Matias Palva, and Satu Palva.  Cross-frequency synchronization connects networks of fast and slow oscillations during visual working memory maintenance.  Elife, 5:e13451, 2016.
  * [30] Aman B Saleem, Anthony D Lien, Michael Krumin, Bilal Haider, Miroslav Roman Roson, Asli Ayaz, Kimberly Reinhold, Laura Busse, Matteo Carandini, and Kenneth D Harris.  Subcortical source and modulation of the narrowband gamma oscillation in mouse visual cortex.  Neuron, 93(2):315–322, 2017.
  * [31] Matthieu Gilson, Anthony N Burkitt, David B Grayden, Doreen A Thomas, and J Leo van Hemmen.  Emergence of network structure due to spike-timing-dependent plasticity in recurrent neuronal networks iii: Partially connected neurons driven by spontaneous activity.  Biological Cybernetics, 101(5-6):411, 2009.
  * [32] Robert Gütig, Ranit Aharonov, Stefan Rotter, and Haim Sompolinsky.  Learning input correlations through nonlinear temporally asymmetric hebbian plasticity.  Journal of Neuroscience, 23(9):3697–3714, 2003.
  * [33] Abigail Morrison, Markus Diesmann, and Wulfram Gerstner.  Phenomenological models of synaptic plasticity based on spike timing.  Biological cybernetics, 98(6):459–478, 2008.
  * [34] Sen Song, Kenneth D Miller, and Larry F Abbott.  Competitive hebbian learning through spike-timing-dependent synaptic plasticity.  Nature neuroscience, 3(9):919, 2000.
  * [35] Maoz Shamir.  Theories of rhythmogenesis.  Current opinion in neurobiology, 58:70–77, 2019.
  * [36] Yotam Luz and Maoz Shamir.  The effect of stdp temporal kernel structure on the learning dynamics of single excitatory and inhibitory synapses.  PloS one, 9(7):e101109, 2014.
  * [37] Abigail Morrison, Markus Diesmann, and Wulfram Gerstner.  Phenomenological models of synaptic plasticity based on spike timing.  Biological cybernetics, 98(6):459–478, 2008.
  * [38] Richard Kempter, Wulfram Gerstner, and J Leo van Hemmen.  Intrinsic stabilization of output rates by spike-based hebbian learning.  Neural computation, 13(12):2709–2741, 2001.
  * [39] Richard Kempter, Wulfram Gerstner, and J Leo Van Hemmen.  Hebbian learning and spiking neurons.  Physical Review E, 59(4):4498, 1999.
  * [40] Yotam Luz and Maoz Shamir.  Balancing feed-forward excitation and inhibition via hebbian inhibitory synaptic plasticity.  PLoS computational biology, 8(1):e1002334, 2012.
  * [41] Henry Markram, Joachim Lübke, Michael Frotscher, and Bert Sakmann.  Regulation of synaptic efficacy by coincidence of postsynaptic aps and epsps.  Science, 275(5297):213–215, 1997.
  * [42] Guo-qiang Bi and Mu-ming Poo.  Synaptic modifications in cultured hippocampal neurons: dependence on spike timing, synaptic strength, and postsynaptic cell type.  Journal of neuroscience, 18(24):10464–10472, 1998.
  * [43] Per Jesper Sjöström, Gina G Turrigiano, and Sacha B Nelson.  Rate, timing, and cooperativity jointly determine cortical synaptic plasticity.  Neuron, 32(6):1149–1164, 2001.
  * [44] Li I Zhang, Huizhong W Tao, Christine E Holt, William A Harris, and Mu-ming Poo.  A critical window for cooperation and competition among developing retinotectal synapses.  Nature, 395(6697):37, 1998.
  * [45] L. F. Abbott and Sacha B. Nelson.  Synaptic plasticity: taming the beast.  Nature Neuroscience, 3(11):1178–1183, 2000.
  * [46] Robert C Froemke, Ishan A Tsay, Mohamad Raad, John D Long, and Yang Dan.  Contribution of individual spikes in burst-induced long-term synaptic modification.  Journal of neurophysiology, 95(3):1620–1629, 2006.
  * [47] Makoto Nishiyama, Kyonsoo Hong, Katsuhiko Mikoshiba, Mu-Ming Poo, and Kunio Kato.  Calcium stores regulate the polarity and input specificity of synaptic modification.  Nature, 408:584–588, 2000.
  * [48] Harel Z Shouval, Mark F Bear, and Leon N Cooper.  A unified model of nmda receptor-dependent bidirectional synaptic plasticity.  Proceedings of the National Academy of Sciences, 99(16):10831–10836, 2002.
  * [49] Melanie A Woodin, Karunesh Ganguly, and Mu ming Poo.  Coincident pre- and postsynaptic activity modifies gabaergic synapses by postsynaptic changes in cl− transporter activity.  Neuron, 39(5):807 – 820, 2003.
  * [50] Oren Shriki, David Hansel, and Haim Sompolinsky.  Rate models for conductance-based cortical neuronal networks.  Neural computation, 15(8):1809–1841, 2003.
  * [51] Gianluigi Mongillo, Simon Rumpel, and Yonatan Loewenstein.  Intrinsic volatility of synaptic connections—a challenge to the synaptic trace theory of memory.  Current opinion in neurobiology, 46:7–13, 2017.
  * [52] Alessio Attardo, James E Fitzgerald, and Mark J Schnitzer.  Impermanence of dendritic spines in live adult ca1 hippocampus.  Nature, 523(7562):592, 2015.
  * [53] Jaime Grutzendler, Narayanan Kasthuri, and Wen-Biao Gan.  Long-term dendritic spine stability in the adult cortex.  Nature, 420(6917):812, 2002.
  * [54] Yi Zuo, Aerie Lin, Paul Chang, and Wen-Biao Gan.  Development of long-term dendritic spine stability in diverse regions of cerebral cortex.  Neuron, 46(2):181–189, 2005.
  * [55] Anthony JGD Holtmaat, Joshua T Trachtenberg, Linda Wilbrecht, Gordon M Shepherd, Xiaoqun Zhang, Graham W Knott, and Karel Svoboda.  Transient and persistent dendritic spines in the neocortex in vivo.  Neuron, 45(2):279–291, 2005.
  * [56] Anthony Holtmaat, Linda Wilbrecht, Graham W Knott, Egbert Welker, and Karel Svoboda.  Experience-dependent and cell-type-specific spine growth in the neocortex.  Nature, 441(7096):979, 2006.
  * [57] Yonatan Loewenstein, Annerose Kuras, and Simon Rumpel.  Multiplicative dynamics underlie the emergence of the log-normal distribution of spine sizes in the neocortex in vivo.  Journal of Neuroscience, 31(26):9481–9488, 2011.
  * [58] Yonatan Loewenstein, Uri Yanover, and Simon Rumpel.  Predicting the dynamics of network connectivity in the neocortex.  Journal of Neuroscience, 35(36):12535–12544, 2015.


# Enumeration of Random Walk Positions in $L_{1}$-norm ball in
$\mathbb{Z}^{d}$

Luchen Shi, Will McCance, Hongjie Zeng

(September 2022)

###### Abstract

In this paper, we mainly concerned about deriving the general formula to count
the possible positions of $n$ step random walk in $\mathbb{Z}^{d}$ with unit
length in each step, which we denoted as $|P_{n}^{d}|$. For our results, we
firstly propose a recurrence relation of the counting formula:
$|P_{n}^{d+1}|=|P_{n}^{d}|+2\sum_{k=0}^{n-1}|P_{k}^{d}|$. Next, we propose two
methods in deriving the explicit formula of $|P_{n}^{d}|$ using generating
functions and Faulhaber’s formula. Finally, we reached our main theorem in the
matrix representation of our formula.

##  1 Preliminaries and Notations

In the Euclidean space $\mathbb{R}^{d}$, we define the $L_{p}$-norm for point
$x=(x_{1},x_{2}\cdots,x_{d})$ to be

$||x||_{p}=(|x_{1}|^{p}+|x_{2}|^{p}+\cdots+|x_{d}|^{p})^{\frac{1}{p}}$

and an open ball around the origin with radius $r$ is given by the set

$B^{p}(r)=\\{x\in\mathbb{R}^{n}:||x||_{p}=(|x_{1}|^{p}+|x_{2}|^{p}+\cdots+|x_{d}|^{p})^{\frac{1}{p}}\leq
r\\}$

When $d=2$, ”balls” corresponding to $L_{1}$-norm (often called the taxicab or
Manhattan metric) are bounded by squares with their diagonals parallel to the
coordinate axes; When $d=3$, ”balls” are octahedra with axes-aligned body
diagonals corresponding to $L_{1}$-norm. We will use $\mathbb{Z}^{d}$ to
denote the $d$-dimensional integer grid, i.e., the set of $d$-tuples
$x=(x^{1},\cdot\cdot\cdot,x^{d})$ where $x^{i}$ are all integers. Each $x$ in
$\mathbb{Z}^{d}$ has $2d$ ”nearest neighbors”, points that are at a distance
one away from it. To do simple random walk in $\mathbb{Z}^{d}$, the walker
starts at the origin and at each integer time $n$ moves to one of the nearest
neighbors, each with probability $\frac{1}{2d}$. Let
$S_{n}=(S^{1}_{n},\cdot\cdot\cdot,S^{d}_{n})$ be the position of the walker
after $n$ steps, which is inside the $L_{1}$-norm ball $B^{1}_{n}$ with radius
$n$. Then

$S_{n}=X_{1}+\cdots+X_{n}$

where $X_{i}=(X^{1}_{i},\cdots,X^{d}_{i})$ and the $X_{i}$ are independent
random vectors with $\mathbb{P}\\{X_{i}=y\\}=(2d)^{-1}$ for each
$y\in\mathbb{Z}^{d}$ that is distance one from the origin. Define $W_{n}^{d}$
to be the set containing the sequences of steps in all possible random walks
with $n$ steps:

$W_{n}^{d}=\\{X_{i}=(X^{1}_{i},\cdots,X^{d}_{i})|\,X_{i}\in\\{+1,-1\\}^{d}\,\,\,(i=1,\ldots,n)\\}$

We have $|W_{n}^{d}|=(2d)^{n}$, since each of the $n$ steps has $2d$
possibilities and are independent. Let us call a random walk of $n$ steps an
“$n$-walk.” From $W_{n}^{d}$, we can construct the set $P_{n}^{d}$ containing
all possible positions in $\mathbb{Z}^{d}$ after an $n$-walk by summing each
unit-length step:

$P_{n}^{d}={{\\{}}\sum_{i=1}^{n}X_{i}|X_{i}=(X^{1}_{i},\cdots,X^{d}_{i}),X^{k}_{i}\in\\{+1,-1\\}{{\\}}}$

After these preparations, we are ready to derive the recursive formula for
$|P_{n}^{d}|$.

##  2 Recursive Formula for $|P_{n}^{d}|$

We first investigate the problem of finding a formula for $|P_{n}^{d}|$. Let
us consider the one-dimensional case first. In this case, random variable $X$
takes value $y$ in the set $Y=\\{+1,-1\\}$. For the illustrative purpose, we
may represent the $\mathbb{Z}$ random walk in $\mathbb{Z}^{2}$ as shown in
Figure 1, to show the possible positions after certain numbers of steps:

0-1+1-11-1+1+1-1-202-1+1+1-1-1+1-3-113-1+1+1-1-1+1+1-1-4-2024-1+1+1-1-1+1+1-1-1+1-5-3-1135-1+1+1-1-1+1+1-1-1+1+1-1-6-4-20246

… … …

Figure 1: Representation of $\mathbb{Z}$ Random Walk in $\mathbb{Z}^{2}$

Clearly, $|P_{n}^{1}|=n+1$. Now, suppose we want to see how many different
positions can be reached in $\mathbb{Z}^{2}$ in $n$ steps. After $1$ step,
there are $4$ possible positions. After $2$ steps, there are $9$ possible
positions. We can see this by partitioning $\mathbb{Z}^{2}$ into ”slices” of
$\mathbb{Z}$: on the $x$-axis there are $|P_{2}^{1}|=3$ possible positions; on
either $y=1$ or $y=-1$ there are $|P_{1}^{1}|=2$ possible positions; on either
$y=2$ or $y=-2$ there is $|P_{0}^{1}|=1$ possible position. Therefore

$|P_{2}^{2}|=|P_{2}^{1}|+2|P_{1}^{1}|+2|P_{0}^{1}|=3+2\cdot 2+2\cdot 1=9$
$y=1$$y=2$$y=-1$$y=-2$$x(y=0)$$y$$(-2,0)$$(0,0)$$(2,0)$$(-1,1)$$(1,1)$$(-1,-1)$$(1,-1)$$(0,2)$$(0,-2)$$|P_{2}^{1}|$$|P_{1}^{1}|$$|P_{0}^{1}|$$|P_{1}^{1}|$$|P_{0}^{1}|$

Figure 2: Slicing Partitions of $P^{2}_{2}$

This idea motivates us to generalize. Every position in $\mathbb{Z}^{d+1}$ can
be written as a $(d+1)\times 1$ matrix where each entry is in $\mathbb{Z}$. We
consider the possible positions after $n$ steps. Since each step is a unit
vector, the maximum absolute value of the last entry is $n$. Therefore, the
last entry takes value in the set

$\\{-n,-(n-1),\ldots,-2,-1,0,1,2,\ldots,(n-1),n\\}$

First, consider the positions in which the last entry is $n$ or $-n$. In each
case the number of possible positions is $|P_{0}^{d}|$ (since there can be no
other steps in any other dimension). Next, consider the positions in which the
last entry is $(n-1)$ or $-(n-1)$. In each case the number of possible
positions is $|P_{1}^{d}|$ (there is one step in one of the remaining $d$
dimensions). We do this until we get to the positions in which the last entry
is $1$ or $-1$. In each case the number of possible positions is
$|P_{n-1}^{d}|$, since the other $(n-1)$ steps are in the remaining $d$
dimensions. Finally, if the last entry is $0$, the number of possible
positions is $|P_{n}^{d}|$. Therefore we have the following relation:

$|P_{n}^{d+1}|=2|P_{0}^{d}|+2|P_{1}^{d}|+\ldots+2|P_{n-1}^{d}|+|P_{n}^{d}|$

which, if written more simply, becomes

$|P_{n}^{d+1}|=|P_{n}^{d}|+2\bigg{(}\sum_{k=0}^{n-1}|P_{k}^{d}|\bigg{)}=2\bigg{(}\sum_{k=0}^{n}|P_{k}^{d}|\bigg{)}-|P_{n}^{d}|$
(1)

Using this recurrence relation, we can calculate the following

$\displaystyle|P_{n}^{2}|$
$\displaystyle=2\bigg{(}\sum_{k=0}^{n}|P_{k}^{1}|\bigg{)}-|P_{n}^{1}|$
$\displaystyle=2\sum_{k=0}^{n}(k+1)-(n+1)$
$\displaystyle=2\sum_{k=0}^{n}k+2\sum_{k=0}^{n}1-(n+1)$
$\displaystyle=2\frac{n(n+1)}{2}+2(n+1)-(n+1)$ $\displaystyle=n(n+1)+(n+1)$
$\displaystyle=(n+1)^{2}$

We may guess that the pattern $(n+1)^{d}$ continues, but unfortunately this is
not the case:

$\displaystyle|P_{n}^{3}|$
$\displaystyle=2\bigg{(}\sum_{k=0}^{n}|P_{k}^{2}|\bigg{)}-|P_{n}^{2}|$
$\displaystyle=2\sum_{k=0}^{n}(k+1)^{2}-(n+1)^{2}$
$\displaystyle=2\sum_{k=1}^{n+1}k^{2}-(n+1)^{2}$
$\displaystyle=2\frac{(n+1)(n+1+1)(2(n+1)+1)}{6}-(n+1)^{2}$
$\displaystyle=\frac{(n+1)(n+2)(2n+3)-3(n+1)^{2}}{3}$
$\displaystyle=\frac{(n+1)(2n^{2}+7n+6-3(n+1))}{3}$
$\displaystyle=\frac{(n+1)(2n^{2}+4n+3)}{3}$
$\displaystyle=\frac{2n^{3}+6n^{2}+7n+3}{3}$
$\displaystyle=\frac{2}{3}n^{3}+2n^{2}+\frac{7}{3}n+1$

This suggests that $|P_{n}^{d}|$ has a more complicated form. However, one of
the things we noticed is that $|P_{n}^{d}|$ seems to always be a polynomial in
$n$ with degree $d$. Indeed, this can be verified by induction.

##  3 $|P_{n}^{d}|$ from generating functions

Our first approach to compute $|P^{d}_{n}|$ involves generating functions and
our recursive formula $(1)$. Our goal was to find a generating function
$g^{d}(x)=|P^{d}_{0}|x^{0}+|P^{d}_{1}|x^{1}+|P^{d}_{2}|x^{2}\cdots=\sum_{i=0}^{\infty}|P^{d}_{i}|x^{i}$
and to determine a formula to find $|P^{d}_{n}|$ by expanding $g^{d}(x)$.
First, we desire to know the formula for
$g^{1}(x)=|P^{1}_{0}|x^{0}+|P^{1}_{1}|x^{1}+\cdots=1+2x+3x^{2}+\cdots$. We
know that this is the derivative of $\frac{1}{1-x}=1+x+x^{2}+x^{3}\cdots$.
Therefore, we have

$g^{1}(x)=\frac{1}{(1-x)^{2}}$

Furthermore, we have a recurrence relation that gives terms of $|P^{d}_{n}|$
in terms of $|P^{d-1}_{i}|$ where $i$ is a natural number ranging from 0 to
$n$. Thus, it might be natural to try to find out what $g^{d+1}(x)$ is in
terms of $g^{d}(x)$. In fact, we may do just that.

$\displaystyle g^{d}(x)=$
$\displaystyle|P^{d}_{0}|x^{0}+\,\,|P^{d}_{1}|x^{1}+\,\,|P^{d}_{2}|x^{2}+\,\,|P^{d}_{3}|x^{3}\cdots$
$\displaystyle+2xg^{d}(x)=$
$\displaystyle\qquad\enspace+2|P^{d}_{0}|x^{1}+2|P^{d}_{1}|x^{2}+2|P^{d}_{2}|x^{3}\cdots$
$\displaystyle+2x^{2}g^{d}(x)=$
$\displaystyle\qquad\quad\enspace\qquad\quad\,\,\,+2|P^{d}_{0}|x^{2}+2|P^{d}_{1}|x^{3}\cdots$
$\displaystyle\vdots$

We can see that the coefficients of each expanded term up to $2x^{n}g^{d}(x)$
sums to $|P^{d+1}_{n}|$. Therefore,

$g^{d+1}(x)=g^{d}(x)+\sum^{\infty}_{i=1}2x^{i}g^{d}(x)$
$g^{d+1}(x)=g^{d}(x)\big{(}1+2\sum^{\infty}_{i=1}x^{i}\big{)}$
$g^{d+1}(x)=g^{d}(x)\bigg{(}1+2\frac{x}{1-x}\bigg{)}$
$g^{d+1}(x)=g^{d}(x)\bigg{(}\frac{1+x}{1-x}\bigg{)}$

Given our initial condition, we find:

$g^{d}(x)=\frac{(1+x)^{d-1}}{(1-x)^{d+1}}$ (2)

Next, we must expand this function back into the form of an infinite series
and compute the coefficients of each term of the function. We will examine
this function in two parts: the binomial polynomial in the numerator, and the
infinite series generated by the denominator. Using the binomial expansion
theorem:

$(1+x)^{d-1}=\displaystyle\sum^{d-1}_{l=0}\binom{d-1}{l}x^{l}$

We now turn to the denominator:

$\frac{1}{(1-x)^{d+1}}$

We will manipulate this into another form using our previous background in
generating functions.

$\frac{1}{(1-x)^{d+1}}=\left(\frac{1}{1-x}\right)^{d+1}=(1+x+x^{2}+\cdots)^{d+1}$
$=\underbrace{(1+x+x^{2}+\cdots)(1+x+x^{2}+\cdots)\dots(1+x+x^{2}+\cdots)}_{\text{$d+1$
times}}$

To determine the expansion of the denominator of the function, we must order
and label the $d+1$ polynomials of the form $(1+x+x^{2}+\cdots)$. We are able
to determine the coefficients of this infinite series by formalizing the
process we use to expand it. We choose one term of the form $x^{e_{i}}$ from
each of the polynomials of the form $(1+x+x^{2}+\cdots)$. Clearly, $e_{i}$ is
always a nonnegative integer. We express this polynomial as the product
$\prod^{d+1}_{i=1}x^{e_{i}}$ for all nonnegative integers $e_{i}\geq 0$. The
coefficient of the $x^{m}$ term is the number of terms of the form
$x^{e_{1}}x^{e_{2}}\cdots x^{e_{d+1}}=x^{m}$. By exponent laws, we get
$x^{e_{1}+e_{2}+\cdots+e_{d+1}}=x^{m}$. This is equivalent to asking how many
$e_{1}+e_{2}+\cdots+e_{d+1}=m$ when $\forall i,\;e_{i}\geq 0$. It is a well
known combinatorial result that this is $\displaystyle\binom{d+m}{d}$.
Therefore,

$\displaystyle\frac{1}{(1-x)^{d+1}}=\displaystyle\sum^{\infty}_{m=0}\binom{d+m}{d}x^{m}$

To derive our final expansion, we multiply

$(1+x)^{d-1}\frac{1}{(1-x)^{d+1}}=\displaystyle\sum^{d-1}_{l=0}\binom{d-1}{l}x^{l}\displaystyle\sum^{\infty}_{m=0}\binom{d+m}{d}x^{m}$

To determine the coefficient $|P^{d}_{n}|$ for any $n$, we must sum over all
products of terms $x^{l}x^{m}=x^{l+m}=x^{n}$. After calculating the expansion
of each coefficient, we obtain the final expansion of $|P^{d}_{n}|$ using
generating functions:

$|P^{d}_{n}|={{\sum}}^{d-1}_{k=0}\binom{d-1}{k}\binom{d+n-k}{d}$ (3)

##  4 $|P_{n}^{d}|$ from Faulhaber’s formula

Next, we present another method to find $|P_{n}^{d}|$. Since we know
$|P_{n}^{d}|$ is a polynomial in $n$ with degree $d$, we only need the
coefficients to fully determine $|P_{n}^{d}|$. We may assume

$|P_{n}^{d}|=c(d,d)n^{d}+c(d,d-1)n^{d-1}+\ldots+c(d,1)n+c(d,0)=\mathbf{c}_{d}\cdot\mathbf{n}_{d}$
(4)

where

$\mathbf{c}_{d}=\begin{pmatrix}c(d,d)\\\ c(d,d-1)\\\ \vdots\\\ c(d,1)\\\
c(d,0)\end{pmatrix},\mathbf{n}_{d}=\begin{pmatrix}n^{d}\\\ n^{d-1}\\\
\vdots\\\ n\\\ 1\end{pmatrix}$ (5)

Using this notation, we may rewrite recurrence relation $(6)$ as:

$c(d+1,d+1)n^{d+1}+c(d+1,d)n^{d}+\ldots+c(d+1,1)n+c(d+1,0)$
$\begin{split}&=2\sum_{k=0}^{n}\big{(}c(d,d)k^{d}+c(d,d-1)k^{d-1}+\ldots+c(d,1)k+c(d,0)\big{)}\\\
&-\big{(}c(d,d)n^{d}+c(d,d-1)n^{d-1}+\ldots+c(d,1)n+c(d,0)\big{)}\\\
&=2c(d,d)\sum_{k=0}^{n}k^{d}+2c(d,d-1)\sum_{k=0}^{n}k^{d-1}+\ldots+2c(d,1)\sum_{k=0}^{n}k+2c(d,0)\sum_{k=0}^{n}1\\\
&-\big{(}c(d,d)n^{d}+c(d,d-1)n^{d-1}+\ldots+c(d,1)n+c(d,0)\big{)}\end{split}$

We notice that sums of the form $\displaystyle\sum_{k=0}^{n}k^{d}$ appear. The
explicit expansion of this sum is given by Faulhaber’s formula:

$\sum_{k=0}^{n}k^{d}=1^{d}+2^{d}+\ldots+n^{d}=\frac{n^{d+1}}{d+1}+\frac{1}{2}n^{d}+\frac{1}{d+1}\sum_{k=2}^{d}B_{k}\binom{d+1}{k}n^{d-k+1}$
(6)

where $B_{k}$ represents the $k$th Bernoulli number. We use $(6)$ to find the
relationship between $\mathbf{c}_{d+1}$ and $\mathbf{c}_{d}$:

$\begin{split}&c(d+1,d+1)n^{d+1}+c(d+1,d)n^{d}+\ldots+c(d+1,1)n+c(d+1,0)\\\
&=2c(d,d)\sum_{k=0}^{n}k^{d}+2c(d,d-1)\sum_{k=0}^{n}k^{d-1}+\ldots+2c(d,1)\sum_{k=0}^{n}k+2c(d,0)\sum_{k=0}^{n}1\\\
&-\big{(}c(d,d)n^{d}+c(d,d-1)n^{d-1}+\ldots+c(d,1)n+c(d,0)\big{)}\\\
&=2c(d,d)\bigg{(}\frac{n^{d+1}}{d+1}+\frac{1}{2}n^{d}+\frac{1}{d+1}\sum_{k=2}^{d}B_{k}\binom{d+1}{k}n^{d-k+1}\bigg{)}\\\
&+2c(d,d-1)\bigg{(}\frac{n^{d}}{d}+\frac{1}{2}n^{d-1}+\frac{1}{d}\sum_{k=2}^{d-1}B_{k}\binom{d}{k}n^{d-k+1}\bigg{)}\\\
&+\ldots\\\
&+2c(d,1)\bigg{(}\frac{n^{2}}{2}+\frac{1}{2}n\bigg{)}+2c(d,0)(n+1)\\\
&-c(d,d)n^{d}-c(d,d-1)n^{d-1}-\ldots-c(d,1)n-c(d,0)\\\
&=2c(d,d)\bigg{(}\frac{n^{d+1}}{d+1}+\frac{1}{d+1}\sum_{k=2}^{d}B_{k}\binom{d+1}{k}n^{d-k+1}\bigg{)}\\\
&+2c(d,d-1)\bigg{(}\frac{n^{d}}{d}+\frac{1}{d}\sum_{k=2}^{d-1}B_{k}\binom{d}{k}n^{d-k}\bigg{)}\\\
&+2c(d,d-2)\bigg{(}\frac{n^{d-1}}{d-1}+\frac{1}{d-1}\sum_{k=2}^{d-2}B_{k}\binom{d-1}{k}n^{d-k-1}\bigg{)}\\\
&+2c(d,d-3)\bigg{(}\frac{n^{d-2}}{d-2}+\frac{1}{d-2}\sum_{k=2}^{d-3}B_{k}\binom{d-2}{k}n^{d-k-2}\bigg{)}\\\
&+\ldots\\\
&+2c(d,4)\bigg{(}\frac{n^{5}}{5}+\frac{1}{5}\sum_{k=2}^{4}B_{k}\binom{5}{k}n^{5-k}\bigg{)}\\\
&+2c(d,3)\bigg{(}\frac{n^{4}}{4}+\frac{1}{4}\sum_{k=2}^{3}B_{k}\binom{4}{k}n^{4-k}\bigg{)}\\\
&+2c(d,2)\bigg{(}\frac{n^{3}}{3}+\frac{1}{3}\sum_{k=2}^{2}B_{k}\binom{3}{k}n^{3-k}\bigg{)}\\\
&+2c(d,1)\bigg{(}\frac{n^{2}}{2}\bigg{)}+2c(d,0)(n+1)-c(d,0)\end{split}$

Comparing coefficients of $n^{d+1},n^{d},\ldots,n,1$, we have the following:

$\displaystyle c(d+1,d+1)$ $\displaystyle=\frac{2c(d,d)}{d+1}$ $\displaystyle
c(d+1,d)$ $\displaystyle=\frac{2c(d,d-1)}{d}$ $\displaystyle c(d+1,d-1)$
$\displaystyle=\frac{2B_{2}}{d+1}\binom{d+1}{2}c(d,d)+\frac{2c(d,d-2)}{d-1}$
$\displaystyle c(d+1,d-2)$
$\displaystyle=\frac{2B_{3}}{d+1}\binom{d+1}{3}c(d,d)+\frac{2B_{2}}{d}\binom{d}{2}c(d,d-1)+\frac{2c(d,d-3)}{d-2}$
$\displaystyle\vdots$ $\displaystyle c(d+1,d-j)$
$\displaystyle=\frac{2B_{j+1}}{d+1}\binom{d+1}{j+1}c(d,d)+\frac{2B_{j}}{d}\binom{d}{j}c(d,d-1)+\frac{2B_{j-1}}{d-1}\binom{d-1}{j-1}c(d,d-2)+\ldots$
$\displaystyle+\frac{2B_{2}}{d-j+2}\binom{d-j+2}{2}c(d,d-j+1)+\frac{2c(d,d-j-1)}{d-j}$
$\displaystyle\vdots$
$\displaystyle+\frac{2B_{2}}{5}\binom{5}{2}c(d,4)+\frac{2c(d,2)}{3}$
$\displaystyle c(d+1,2)$
$\displaystyle=\frac{2B_{d-1}}{d+1}\binom{d+1}{d-1}c(d,d)+\frac{2B_{d-2}}{d}\binom{d}{d-2}c(d,d-1)+\frac{2B_{d-3}}{d-1}\binom{d-1}{d-3}c(d,d-2)+\ldots$
$\displaystyle+\frac{2B_{3}}{5}\binom{5}{3}c(d,4)+\frac{2B_{2}}{4}\binom{4}{2}c(d,3)+\frac{2c(d,1)}{2}$
$\displaystyle c(d+1,1)$
$\displaystyle=\frac{2B_{d}}{d+1}\binom{d+1}{d}c(d,d)+\frac{2B_{d-1}}{d}\binom{d}{d-1}c(d,d-1)+\frac{2B_{d-2}}{d-1}\binom{d-1}{d-2}c(d,d-2)+\ldots$
$\displaystyle+\frac{2B_{4}}{5}\binom{5}{4}c(d,4)+\frac{2B_{3}}{4}\binom{4}{3}c(d,3)+\frac{2B_{2}}{3}\binom{3}{2}c(d,2)+2c(d,0)$
$\displaystyle c(d+1,0)$ $\displaystyle=c(d,0)$

Using matrix form we may write this more compactly:

$\mathbf{c}_{d+1}=M_{d}^{d+1}\mathbf{c}_{d}$ (7)

where

$M_{d}^{d+1}=\setcounter{MaxMatrixCols}{11}\begin{pmatrix}\frac{2}{d+1}&0&0&0&0&\cdots&0&\cdots&0&0&0\\\
0&\frac{2}{d}&0&0&0&\cdots&0&\cdots&0&0&0\\\
\frac{2B_{2}}{d+1}\binom{d+1}{2}&0&\frac{2}{d-1}&0&0&\cdots&0&\cdots&0&0&0\\\
\frac{2B_{3}}{d+1}\binom{d+1}{3}&\frac{2B_{2}}{d}\binom{d}{2}&0&\frac{2}{d-2}&0&\cdots&0&\cdots&0&0&0\\\
\frac{2B_{4}}{d+1}\binom{d+1}{4}&\frac{2B_{3}}{d}\binom{d}{3}&\frac{2B_{2}}{d-1}\binom{d-1}{2}&0&\frac{2}{d-3}&\cdots&0&\cdots&0&0&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&&\vdots&\vdots&\vdots\\\
\frac{2B_{j+1}}{d+1}\binom{d+1}{j+1}&\frac{2B_{j}}{d}\binom{d}{j}&\frac{2B_{j-1}}{d-1}\binom{d-1}{j-1}&\cdots&\frac{2B_{2}}{d-j+2}\binom{d-j+2}{2}&0&\frac{2}{d-j}&\cdots&0&0&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots\\\
\frac{2B_{d-2}}{d+1}\binom{d+1}{d-2}&\frac{2B_{d-3}}{d}\binom{d}{d-3}&\frac{2B_{d-4}}{d-1}\binom{d-1}{d-4}&\cdots&\cdots&\cdots&\frac{2B_{2}}{5}\binom{5}{2}&0&\frac{2}{3}&0&0\\\
\frac{2B_{d-1}}{d+1}\binom{d+1}{d-1}&\frac{2B_{d-2}}{d}\binom{d}{d-2}&\frac{2B_{d-3}}{d-1}\binom{d-1}{d-3}&\cdots&\cdots&\cdots&\cdots&\frac{2B_{2}}{4}\binom{4}{2}&0&\frac{2}{2}&0\\\
\frac{2B_{d}}{d+1}\binom{d+1}{d}&\frac{2B_{d-1}}{d}\binom{d}{d-1}&\frac{2B_{d-2}}{d-1}\binom{d-1}{d-2}&\cdots&\cdots&\cdots&\frac{2B_{4}}{5}\binom{5}{4}&\frac{2B_{3}}{4}\binom{4}{3}&\frac{2B_{2}}{3}\binom{3}{2}&0&2\\\
0&0&0&0&0&\cdots&0&\cdots&0&0&1\end{pmatrix}$

is a $(d+2)\times(d+1)$ matrix.
For example,

$\displaystyle\mathbf{c}_{1}$ $\displaystyle=\begin{pmatrix}c(1,1)\\\
c(1,0)\end{pmatrix}=\begin{pmatrix}1\\\ 1\end{pmatrix}$
$\displaystyle\mathbf{c}_{2}$
$\displaystyle=M_{1}^{2}\mathbf{c}_{1}=\begin{pmatrix}1&0\\\ 0&2\\\ 0&1\\\
\end{pmatrix}\begin{pmatrix}1\\\ 1\end{pmatrix}=\begin{pmatrix}1\\\ 2\\\
1\end{pmatrix}$ $\displaystyle\mathbf{c}_{3}$
$\displaystyle=M_{2}^{3}\mathbf{c}_{2}=\begin{pmatrix}\frac{2}{3}&0&0\\\
0&1&0\\\ \frac{1}{3}&0&2\\\ 0&0&1\\\ \end{pmatrix}\begin{pmatrix}1\\\ 2\\\
1\end{pmatrix}=\begin{pmatrix}\frac{2}{3}\\\ 2\\\ \frac{7}{3}\\\
1\end{pmatrix}$ $\displaystyle\mathbf{c}_{4}$
$\displaystyle=M_{3}^{4}\mathbf{c}_{3}=\begin{pmatrix}\frac{1}{2}&0&0&0\\\
0&\frac{2}{3}&0&0\\\ \frac{1}{2}&0&1&0\\\ 0&\frac{1}{3}&0&2\\\ 0&0&0&1\\\
\end{pmatrix}\begin{pmatrix}\frac{2}{3}\\\ 2\\\ \frac{7}{3}\\\
1\end{pmatrix}=\begin{pmatrix}1/3\\\ 4/3\\\ 8/3\\\ 8/3\\\ 1\end{pmatrix}$

which gives, correspondingly,

$\displaystyle|P_{n}^{1}|$ $\displaystyle=n+1$ $\displaystyle|P_{n}^{2}|$
$\displaystyle=n^{2}+2n+1$ $\displaystyle|P_{n}^{3}|$
$\displaystyle=\frac{2}{3}n^{3}+2n^{2}+\frac{7}{3}n+1$
$\displaystyle|P_{n}^{4}|$
$\displaystyle=\frac{1}{3}n^{4}+\frac{4}{3}n^{3}+\frac{8}{3}n^{2}+\frac{8}{3}n+1$

This is summarized in the following theorem:

######  Theorem 1.

The cardinality of the set of possible positions in $\mathbb{Z}^{d}$ after a
$n$-walk, $|P_{n}^{d}|$, is given by

$|P_{n}^{d}|=\mathbf{c}_{d}\cdot\mathbf{n}_{d}$ (8)

where $\mathbf{c}_{d}$ and $\mathbf{n}_{d}$ are as given in $(11)$, and
$\mathbf{c}_{d}$ is calculated recursively by

$\mathbf{c}_{d}=M_{d-1}^{d}M_{d-2}^{d-1}\cdots
M_{2}^{3}M_{1}^{2}\mathbf{c}_{1}$ (9)

where

$\mathbf{c}_{1}=\begin{pmatrix}1\\\ 1\\\ \end{pmatrix}$

and $M_{d}^{d+1}$ is given above.

##  5 Comparing $c(d,d),c(d,d-1),c(d,d-2)$ from two approaches

Since our two methods of finding $|P_{n}^{d}|$ both originate from the
recurrence relation (1), they must yield the same result for any nonnegrative
integer $d$. In fact, this was proved to be true from $d=1$ to $5$. During the
process of proving this, we found that although our second method involving
matrices seems to be complicated and does not provide a formula for
$|P_{n}^{d}|$, when actually carrying out calculation, it is much simpler than
expanding the sum obtained by our first method, as only basic arithmetic of
rational numbers is involved. We try to derive explicitly from both methods
the coefficient of each power of $n$ in $|P_{n}^{d}|$ and show that they are
equal. Unfortunately, this turns out to be difficult for coefficient of the
most general form $c(d,d-j)$ in $|P_{n}^{d}|$ and we pose this as an open
problem. However, for $c(d,d),c(d,d-1)$ and $c(d,d-2)$ we are able to prove
that both methods give the same result,as we now show.

###  5.1 $c(d,d)$

The first row of $M_{d}^{d+1}$ tells us

$c(d+1,d+1)=\frac{2}{d+1}c(d,d)$

which, if we develop further, yields

$\displaystyle c(d+1,d+1)$ $\displaystyle=\frac{2}{d+1}c(d,d)$
$\displaystyle=\frac{2}{d+1}\frac{2}{d}c(d-1,d-1)$ $\displaystyle=\cdots$
$\displaystyle=\frac{2}{d+1}\frac{2}{d}\cdots\frac{2}{2}c(1,1)$
$\displaystyle=\frac{2^{d}}{(d+1)!}$

Therefore,

$c(d,d)=\frac{2^{d-1}}{d\,!}$ (10)

Expanding $(1)$ gives

$|P_{n}^{d}|={\sum}_{k=0}^{d-1}\binom{d-1}{k}\frac{(n-k+d)\,!}{d\,!\,(n-k)\,!}=\frac{1}{d\,!}{\sum}_{k=0}^{d-1}\binom{d-1}{k}(n-(k-1))(n-(k-2))\cdots(n-(k-d))$
(11)

From this, and by the binomial theorem, the coefficient of $n^{d}$ is

$c(d,d)=\frac{1}{d\,!}{\sum}_{k=0}^{d-1}\binom{d-1}{k}=\frac{2^{d-1}}{d\,!}$

which agrees with $(10)$.

###  5.2 $c(d,d-1)$

The second row of $M_{d}^{d+1}$ tells us

$c(d+1,d)=\frac{2}{d}c(d,d-1)$

which we may develop as before:

$\displaystyle c(d+1,d)$ $\displaystyle=\frac{2}{d}c(d,d-1)$
$\displaystyle=\frac{2}{d}\frac{2}{d-1}c(d-1,d-2)$ $\displaystyle=\cdots$
$\displaystyle=\frac{2}{d}\frac{2}{d-1}\cdots\frac{2}{2}\frac{2}{1}c(1,0)$
$\displaystyle=\frac{2^{d}}{d\,!}$

and thus

$c(d,d-1)=\frac{2^{d-1}}{(d-1)\,!}$ (12)

Now we look at $(19)$ again and try to determine the coefficient of $n^{d-1}$.
To contribute to $c(d,d-1)n^{d-1}$, we may take $-(k-i)$ out of the $i$th
factor and $n$ out of the other $(d-1)$ factors; since $i$ ranges from $1$ to
$d$, there are $d$ ways to do this, and we can write

$\displaystyle c(d,d-1)$
$\displaystyle=\frac{1}{d\,!}{\sum}_{k=0}^{d-1}\binom{d-1}{k}{\sum}_{i=1}^{d}(-(k-i))$
$\displaystyle=\frac{1}{d\,!}{\sum}_{k=0}^{d-1}\binom{d-1}{k}\bigg{(}\frac{d(d+1)}{2}-kd\bigg{)}$
$\displaystyle=\frac{d(d+1)2^{d-1}}{2d\,!}-\frac{1}{(d-1)\,!}{\sum}_{k=0}^{d-1}k\binom{d-1}{k}$
$\displaystyle=\frac{(d+1)2^{d-2}}{(d-1)\,!}-\frac{2^{d-2}}{(d-2)\,!}$
$\displaystyle=\frac{2^{d-2}}{(d-2)\,!}\bigg{(}\frac{d+1}{d-1}-1\bigg{)}$
$\displaystyle=\frac{2^{d-1}}{(d-1)\,!}$

(In the process we used the identity
$k\displaystyle\binom{n}{k}=n\displaystyle\binom{n-1}{k-1}$)
which agrees with $(12)$.

###  5.3 $c(d,d-2)$

The third row of $M_{d}^{d+1}$ tells us

$c(d+1,d-1)=\frac{2B_{2}}{d+1}\binom{d+1}{2}c(d,d)+\frac{2}{d-1}c(d,d-2)$

which we may develop further to obtain

$\displaystyle c(d+1,d-1)$
$\displaystyle=\frac{d}{6}c(d,d)+\frac{2}{d-1}c(d,d-2)$
$\displaystyle=\frac{2^{d-1}}{6(d-1)\,!}+\frac{2}{d-1}c(d,d-2)$
$\displaystyle=\cdots$
$\displaystyle=(d-1)\frac{2^{d-1}}{6(d-1)\,!}+\frac{2^{d-1}}{(d-1)(d-2)\cdots(d-(d-1))}c(2,0)$
$\displaystyle=\frac{2^{d-1}}{6(d-2)\,!}+\frac{2^{d-1}}{(d-1)\,!}$
$\displaystyle=\frac{2^{d-1}}{6(d-2)\,!}\bigg{(}1+\frac{6}{d-1}\bigg{)}$
$\displaystyle=\frac{2^{d-1}(d+5)}{6(d-1)\,!}$

and thus

$c(d,d-2)=\frac{2^{d-2}(d+4)}{6(d-2)\,!}$ (13)

Before we turn to finding $c(d,d-2)$ from $(11)$, we first prove three lemmas
which are going to be used.

######  Lemma 1.

$\displaystyle\sum_{k=1}^{n}k\binom{n}{k}=n2^{n-1}$ and
$\displaystyle\sum_{k=1}^{n}k^{2}\binom{n}{k}=n(n+1)2^{n-2}$

###### Proof.

We begin with a special case of the Binomial Theorem:

$(1+x)^{n}=\sum_{k=0}^{n}\binom{n}{k}x^{k}$ (14)

Differentiating both sides of $(14)$ with respect to $x$, we get

$n(1+x)^{n-1}=\sum_{k=1}^{n}k\binom{n}{k}x^{k-1}$ (15)

Substituting $x=1$ in $(15)$, we get

$\sum_{k=1}^{n}k\binom{n}{k}=n2^{n-1}$ (16)

which is the first identity. Multiplying both sides of $(15)$ by $x$, we get

$nx(1+x)^{n-1}=\sum_{k=1}^{n}k\binom{n}{k}x^{k}$ (17)

Differentiating both sides of $(17)$ with respect to $x$, we get

$n\big{(}(1+x)^{n-1}+(n-1)x(1+x)^{n-2}\big{)}=\sum_{k=1}^{n}k^{2}\binom{n}{k}x^{k-1}$
(18)

Substituting $x=1$ in $(18)$ finally gives

$\sum_{k=1}^{n}k^{2}\binom{n}{k}=n\big{(}2^{n-1}+(n-1)2^{n-2}\big{)}=n(n+1)2^{n-2}$
(19)

∎

######  Lemma 2.

$\displaystyle\sum_{1\leq i<j\leq n}(i+j)=\frac{n(n+1)(n-1)}{2}$

######  Lemma 3.

$\displaystyle\sum_{1\leq i<j\leq n}ij=\frac{n(n+1)(3n+2)(n-1)}{24}$

We are now ready to proceed to finding $c(d,d-2)$ from $(11)$.

$\displaystyle c(d,d-2)$
$\displaystyle=\frac{1}{d\,!}{\sum}_{k=0}^{d-1}\binom{d-1}{k}(-1)^{2}{\sum}_{1\leq
i<j\leq d}(k-i)(k-j)$
$\displaystyle=\frac{1}{d\,!}{\sum}_{k=0}^{d-1}\binom{d-1}{k}{\sum}_{1\leq
i<j\leq d}{(}k^{2}-(i+j)k+ij{)}$
$\displaystyle=\frac{1}{d\,!}{\sum}_{k=0}^{d-1}\binom{d-1}{k}\bigg{(}\binom{d}{2}k^{2}-k\cdot\frac{d(d+1)(d-1)}{2}+\frac{d(d+1)(3d+2)(d-1)}{24}\bigg{)}$
$\displaystyle=\frac{1}{d\,!}\bigg{(}\binom{d}{2}{\sum}_{k=0}^{d-1}k^{2}\binom{d-1}{k}-\frac{d(d+1)(d-1)}{2}{\sum}_{k=0}^{d-1}k\binom{d-1}{k}+\frac{d(d+1)(3d+2)(d-1)}{24}{\sum}_{k=0}^{d-1}\binom{d-1}{k}\bigg{)}$
$\displaystyle=\frac{2^{d-3}}{6(d-2)\,!}(2d+8)=\frac{2^{d-2}(d+4)}{6(d-2)\,!}$

which agrees with $(13)$.

###  5.4 What’s next?

Using similar procedure, from the matrix $M_{d}^{d+1}$ we can show that

$\displaystyle c(d,d-3)$ $\displaystyle=\frac{2^{d-2}d}{(d-3)\,!}$ (20)
$\displaystyle c(d,d-4)$
$\displaystyle=\frac{2^{d-6}(5d^{2}+33d-32)}{45(d-4)\,!}$ (21)

However, if we expand $(19)$, we get the following:

$\displaystyle c(d,d-3)$
$\displaystyle=\frac{1}{d\,!}{\sum}_{k=0}^{d-1}\binom{d-1}{k}(-1)^{3}{\sum}_{1\leq
i_{1}<i_{2}<i_{3}\leq d}(k-i_{1})(k-i_{2})(k-i_{3})$ (22) $\displaystyle
c(d,d-4)$
$\displaystyle=\frac{1}{d\,!}{\sum}_{k=0}^{d-1}\binom{d-1}{k}(-1)^{4}{\sum}_{1\leq
i_{1}<i_{2}<i_{3}<i_{4}\leq d}(k-i_{1})(k-i_{2})(k-i_{3})(k-i_{4})$ (23)

Of course, we have reason to believe that $(22)$ ends up being $(20)$, and
$(23)$ ends up being $(21)$. However, a major difficulty appears, which is to
evaluate the sums of the form

${\sum}_{1\leq i_{1}<i_{2}<\cdots<i_{j}\leq
d}(k-i_{1})(k-i_{2})\cdots(k-i_{j})$ (24)

which appears in the expression for $c(d,d-j)$ obtained from $(11)$:

$c(d,d-j)=\frac{1}{d\,!}{\sum}_{k=0}^{d-1}\binom{d-1}{k}(-1)^{j}{\sum}_{1\leq
i_{1}<i_{2}<\cdots<i_{j}\leq d}(k-i_{1})(k-i_{2})\cdots(k-i_{j})$ (25)

One possible way to simplify the sum $(24)$ is by the substitution

$i_{1}^{\prime}=k-i_{j},i_{2}^{\prime}=k-i_{j-1},\cdots,i_{j}^{\prime}=k-i_{1}$

So that

${\sum}_{1\leq i_{1}<i_{2}<\cdots<i_{j}\leq
d}(k-i_{1})(k-i_{2})\cdots(k-i_{j})={\sum}_{k-d\leq
i_{1}^{\prime}<i_{2}^{\prime}<\cdots<i_{j}^{\prime}\leq
k-1}i_{1}^{\prime}i_{2}^{\prime}\cdots i_{j}^{\prime}$

This requires us to do sums of the form

${\sum}_{a\leq i_{1}<i_{2}<\cdots<i_{j}\leq b}i_{1}i_{2}\cdots
i_{j}\,\,\,\,\,(a,b\in\mathbb{Z},a<b)$

Of course, we have

$\displaystyle{\sum}_{a\leq i_{1}<i_{2}<\cdots<i_{j}\leq b}i_{1}i_{2}\cdots
i_{j}$ $\displaystyle=a{\sum}_{a+1\leq i_{2}<\cdots<i_{j}\leq b}i_{2}\cdots
i_{j}$ $\displaystyle+(a+1){\sum}_{a+2\leq i_{2}<\cdots<i_{j}\leq
b}i_{2}\cdots i_{j}$ $\displaystyle+\cdots+(b-j){\sum}_{b-j+1\leq
i_{2}<\cdots<i_{j}\leq b}i_{2}\cdots i_{j}$
$\displaystyle+(b-j+1){\sum}_{b-j+2\leq i_{2}<\cdots<i_{j}\leq b}i_{2}\cdots
i_{j}$ $\displaystyle={\sum}_{s=a}^{b-j+1}s{\sum}_{s+1\leq
i_{2}<\cdots<i_{j}\leq b}i_{2}\cdots i_{j}$

which gives a way to recursively evaluate the sum from the previous sum of the
same form. However, when carried out in practice, this turns out to be
increasingly complicated as $j$ gets larger. An interesting phenomenon that
appears is that Faulhaber’s formula is also involved in such calculations,
which might be related to the fact that the formula is also an indispensable
part in obtaining the matrix $M_{d}^{d+1}$.
On the other hand, when we try to find more coefficients from the matrix
$M_{d}^{d+1}$, the main difficulty comes from two aspects: Firstly, we need
information on coefficients $c(d,d),\cdots c(d,d-j+2)$ in order to find
$c(d,d-j)$; Secondly, since there is no closed-form formula for the Bernoulli
numbers, as they are involved in calculations, it is difficult to spot the
pattern in which the coefficients evolve. We currently do not have a good way
to overcome these difficulties and consequently propose the following open
problem.

######  Open Problem 1.

Derive expressions for $c(d,d-j)$ from both the matrix $M_{d}^{d+1}$ and
$(11)$ and prove that they are indeed equal.

## References

  * [1] Donald E. Knuth (1993). ”Johann Faulhaber and sums of powers”. Mathematics of Computation. 61 (203): 277–294. arXiv:math.CA/9207222
  * [2] Introductory Combinatorics (5th Edition), Richard A.Brualdi, Pearson
  * [3] Lectures on Contemporary Probability, Gregory F.Lawler, Lester N.Coyle, American Mathematical Society, Institute for Advanced Study
  * [4] Lectures on Contemporary Probability, Gregory F.Lawler, Lester N.Coyle, American Mathematical Society, Institute for Advanced Study
  * [5] Bourbaki, Nicolas (1987) [1981]. Topological Vector Spaces: Chapters 1–5. Éléments de mathématique. Translated by Eggleston, H.G.; Madan, S. Berlin New York: Springer-Verlag. ISBN 3-540-13627-4. OCLC 17499190.
  * [6] Khaleelulla, S. M. (1982). Counterexamples in Topological Vector Spaces. Lecture Notes in Mathematics. Vol. 936. Berlin, Heidelberg, New York: Springer-Verlag. ISBN 978-3-540-11565-6. OCLC 8588370.
  * [7] Narici, Lawrence; Beckenstein, Edward (2011). Topological Vector Spaces. Pure and applied mathematics (Second ed.). Boca Raton, FL: CRC Press. ISBN 978-1584888666. OCLC 144216834.
  * [8] Schaefer, Helmut H.; Wolff, Manfred P. (1999). Topological Vector Spaces. GTM. Vol. 8 (Second ed.). New York, NY: Springer New York Imprint Springer. ISBN 978-1-4612-7155-0. OCLC 840278135.
  * [9] Trèves, François (2006) [1967]. Topological Vector Spaces, Distributions and Kernels. Mineola, N.Y.: Dover Publications. ISBN 978-0-486-45352-1. OCLC 853623322.
  * [10] Wilansky, Albert (2013). Modern Methods in Topological Vector Spaces. Mineola, New York: Dover Publications, Inc. ISBN 978-0-486-49353-4. OCLC 849801114.