# Mathematical derivation of wave propagation properties
in hierarchical neural networks
with predictive coding feedback dynamics

Grégory Faye111<EMAIL_ADDRESS>Institut de Mathématiques
de Toulouse ; UMR5219, Université de Toulouse ; UPS IMT, F-31062 Toulouse
Cedex 9 France  ANITI, Université de Toulouse, 31062, Toulouse, France
Guilhem Fouilhé  Institut de Mathématiques de Toulouse ; UMR5219, Université
de Toulouse ; UPS IMT, F-31062 Toulouse Cedex 9 France  CerCo, CNRS, 31052
Toulouse, France  ANITI, Université de Toulouse, 31062, Toulouse, France
Rufin VanRullen  CerCo, CNRS, 31052 Toulouse, France  ANITI, Université de
Toulouse, 31062, Toulouse, France

###### Abstract

Sensory perception (e.g. vision) relies on a hierarchy of cortical areas, in
which neural activity propagates in both directions, to convey information not
only about sensory inputs but also about cognitive states, expectations and
predictions. At the macroscopic scale, neurophysiological experiments have
described the corresponding neural signals as both forward and backward-
travelling waves, sometimes with characteristic oscillatory signatures. It
remains unclear, however, how such activity patterns relate to specific
functional properties of the perceptual apparatus. Here, we present a
mathematical framework, inspired by neural network models of predictive
coding, to systematically investigate neural dynamics in a hierarchical
perceptual system. We show that stability of the system can be systematically
derived from the values of hyper-parameters controlling the different signals
(related to bottom-up inputs, top-down prediction and error correction).
Similarly, it is possible to determine in which direction, and at what speed
neural activity propagates in the system. Different neural assemblies
(reflecting distinct eigenvectors of the connectivity matrices) can
simultaneously and independently display different properties in terms of
stability, propagation speed or direction. We also derive continuous-limit
versions of the system, both in time and in neural space. Finally, we analyze
the possible influence of transmission delays between layers, and reveal the
emergence of oscillations at biologically plausible frequencies.

##  1 Introduction

The brain’s anatomy is characterized by a strongly hierarchical architecture,
with a succession of brain regions that process increasingly complex
information. This functional strategy is mirrored by the succession of
processing layers found in modern deep neural networks (and for this reason,
we use the term “layer” in this work to denote one particular brain region in
this hierarchy, rather than the laminar organization of cortex that is well-
known to neuroscientists). The hierarchical structure is especially obvious in
the organization of the visual system [15], starting from the retina through
primary visual cortex (V1) and various extra-striate regions, and culminating
in temporal lobe regions for object recognition and in parietal regions for
motion and location processing.
In this hierarchy of brain regions, the flow of information is clearly
bidirectional: there are comparable number of fibers sending neural signals
down (from higher to lower levels of the hierachy) as there are going up [8].
While the bottom-up or “feed-forward” propagation of information is easily
understood as integration of sensory input (and matches the functional
structure found in artificial deep learning networks), the opposite feedback
direction of propagation is more mysterious, and its functional role remains
unknown.
Predictive coding is one dominant theory to explain the function of cortical
feedback [27]. Briefly, the theory states that each layer in the cortical
hierarchy generates predictions about what caused their own activity; these
predictions are sent to the immediately preceding layer, where a prediction
error can be computed, and carried forward to the original layer, which can
then iteratively update its prediction. Over time (and as long as the sensory
input does not change), the system settles into a state where top-down
predictions agree with bottom-up inputs, and no prediction error is
transmitted. Like any large-scale theory of brain function, the predictive
coding theory is heavily debated [22]. But macroscopic (EEG) experiments have
revealed characteristic propagation signatures that could be hallmarks of
predictive coding. For instance, Alamia and VanRullen [1] showed evidence for
alpha-band (7-15Hz) oscillatory travelling waves propagating in both
directions (feed-forward and feedback); the oscillation frequency and dynamics
were compatible with a simplistic hierarchical model that included a
biologically plausible time delay for transmitting signals between layers, and
was also confirmed by a rudimentary mathematical model. In another study,
Bastos et al [4, 5] found that beta (15-30Hz) and gamma-frequency (30-100Hz)
oscillations could reflect, respectively, the predictions and prediction
errors signals carried by backward and forward connections.
More recently, predictive coding has been explored in the context of deep
neural networks [32, 9, 25]. For instance, Choksi et al [9] augmented existing
deep convolutional networks with feedback connections and a mechanism for
computing and minimizing prediction errors, and found that the augmented
system displayed more robust perception, better aligned with human abilities.
In another study, Pang et al [25] used a similar system and reported the
emergence of illusory contour perception comparable to what humans (but not
standard deep neural networks) would typically perceive.
While the concept of predictive coding is potentially fundamental for
understanding brain function, and its large-scale implementation in deep
artificial neural networks provides empirical support for its potential
functional relevance, there is a gap of theoretical knowledge about the type
of brain activity that predictive coding could engender, and the potential
conditions for its stability. Here, we propose a mathematical framework where
a potentially infinite number of neuronal layers exchange signals in both
directions according to predictive coding principles. The stable propagation
of information in such a system can be explored analytically as a function of
its initial state, its internal parameters (controlling the strength of
inputs, predictions, and error signals) and its connectivity (e.g. convolution
kernels). Our approach considers both a discrete approximation of the system,
as well as continuous abstractions. We demonstrate the practical relevance of
our findings by applying them to a ring model of orientation processing.
Finally, we extend our analytical framework to a more biologically plausible
situation with communication delays between successive layers. This gives rise
to oscillatory signals resembling those observed in the brain.

##  2 Model description

Our initial model is inspired by the generic formulation of predictive coding
proposed in the context of deep learning models by Choksi et al. [9]. This
formulation considers different update terms at each time step: feed-forward
inputs, memory term, feedback- and feed-forward prediction error corrections.
By modulating the hyper-parameters controlling each of these terms, the model
can be reconciled with different formulations of predictive coding (for
instance, the Rao and Ballard model [27] by setting the feed-forward input
term to zero) or other models of hierarchical brain function (e.g. similar to
Heeger’s model [19] by setting the feed-forward error correction to zero).
Indeed, our objective is precisely to characterize the propagation dynamics
inside the network as a function of the relative value of these hyper-
parameters, which in turn alters the model’s functionality.

We consider the following recurrence equation where
$\mathcal{E}_{j}^{n}\in\mathbb{R}^{d}$ represents an encoder at step $n$ and
layer $j$

$\mathcal{E}_{j}^{n+1}=\beta\mathcal{W}^{f}\mathcal{E}_{j-1}^{n+1}+(1-\beta)\mathcal{E}_{j}^{n}-\alpha\mathcal{F}_{j-1}^{n}-\lambda\mathcal{B}_{j}^{n},\quad
j=1,\cdots,J-1\,,$

where $\mathcal{W}^{f}\in\mathscr{M}_{d}(\mathbb{R})$ is a $d\times d$ square
matrix representing the weights of feedforward connections which we assume to
be the same for each layer such that $\mathcal{W}^{f}\mathcal{E}_{j-1}^{n+1}$
models an instantaneous feedforward drive from layer $j-1$ to layer $j$,
controlled by hyper-parameter $\beta$. The term $\mathcal{F}_{j-1}^{n}$
encodes a feedforward error correction process, controlled by hyper-parameter
$\alpha$, where the reconstruction error $\mathcal{R}_{j-1}^{n}$ at layer
$j-1$, defined as the square error between the representation
$\mathcal{E}_{j-1}^{n}$ and the predicted reconstruction
$\mathcal{W}^{b}\mathcal{E}_{j}^{n}$, that is

$\mathcal{R}_{j-1}^{n}:=\frac{1}{2}\|\mathcal{E}_{j-1}^{n}-\mathcal{W}^{b}\mathcal{E}_{j}^{n}\|^{2},$

propagates to the layer $j$ to update its representation. Here,
$\mathcal{W}^{b}\in\mathscr{M}_{d}(\mathbb{R})$ is a $d\times d$ square matrix
representing the weights of feedback connections which we assume to be the
same for each layer. Following [27, 9, 32, 1], the contribution
$\mathcal{F}_{j-1}^{n}$ is then taken to be the gradient of
$\mathcal{R}_{j-1}^{n}$ with respect to $\mathcal{E}_{j}^{n}$, that is

$\mathcal{F}_{j-1}^{n}=\nabla\mathcal{R}_{j-1}^{n}=-(\mathcal{W}^{b})^{\mathbf{t}}\mathcal{E}_{j-1}^{n}+(\mathcal{W}^{b})^{\mathbf{t}}\mathcal{W}^{b}\mathcal{E}_{j}^{n}.$

On the other hand, $\mathcal{B}_{j}^{n}$ incorporates a top-down prediction to
update the representation at layer $j$. This term thus reflects a feedback
error correction process, controlled by hyper-parameter $\lambda$. Similar to
the feedforward process, $\mathcal{B}_{j}^{n}$ is defined as the the gradient
of $\mathcal{R}_{j}^{n}$ with respect to $\mathcal{E}_{j}^{n}$, that is

$\mathcal{B}_{j}^{n}=\nabla\mathcal{R}_{j}^{n}=-\mathcal{W}^{b}\mathcal{E}_{j+1}^{n}+\mathcal{E}_{j}^{n}.$

As a consequence, our model reads

$\mathcal{E}_{j}^{n+1}=\beta\mathcal{W}^{f}\mathcal{E}_{j-1}^{n+1}+\alpha(\mathcal{W}^{b})^{\mathbf{t}}\mathcal{E}_{j-1}^{n}+\left[(1-\beta-\lambda)\mathbf{I}_{d}-\alpha(\mathcal{W}^{b})^{\mathbf{t}}\mathcal{W}^{b}\right]\mathcal{E}_{j}^{n}+\lambda\mathcal{W}^{b}\mathcal{E}_{j+1}^{n},$
(2.1)

for each $j=1,\cdots,J-1$ and $n\geq 0$ where we denoted $\mathbf{I}_{d}$ the
identity matrix of $\mathscr{M}_{d}(\mathbb{R})$. We supplement the recurrence
equation (2.1) with the following boundary conditions at layer $j=0$ and layer
$j=J$. First, at layer $j=0$, we impose

$\mathcal{E}_{0}^{n}=\mathcal{S}_{0}^{n},\quad n\geq 0,$ (2.2)

where $\mathcal{S}_{0}^{n}\in\mathbb{R}^{d}$ is a given source term, which can
be understood as the network’s constant visual input. At the final layer
$j=J$, there is no possibility of incoming top-down signal, and thus one gets

$\mathcal{E}_{J}^{n+1}=\beta\mathcal{W}^{f}\mathcal{E}_{J-1}^{n+1}+\alpha(\mathcal{W}^{b})^{\mathbf{t}}\mathcal{E}_{J-1}^{n}+\left[(1-\beta)\mathbf{I}_{d}-\alpha(\mathcal{W}^{b})^{\mathbf{t}}\mathcal{W}^{b}\right]\mathcal{E}_{J}^{n},\quad
n\geq 0.$ (2.3)

Finally, at the initial step $n=0$, we set

$\mathcal{E}^{0}_{j}=\mathcal{H}_{j},\quad j=0,\cdots,J,$ (2.4)

for some given initial sequence $(\mathcal{H}_{j})_{0,\cdots,J}$. For
instance, in Choksi et al [9], $\mathcal{H}_{j}$ was initialized by a first
feedforward pass through the system, i.e. $\beta>0$ and $\alpha=\lambda=0$.
Throughout we assume the natural following compatibility condition between the
source terms and the initial condition, namely

$\mathcal{S}_{0}^{0}=\mathcal{H}_{0}.$ (2.5)

Regarding the hyper-parameters of the problem we assume that

$0\leq\beta<1,\quad\text{ with }\quad 0\leq\alpha+\lambda\leq 1.$ (2.6)

Our key objective is to characterize the behavior of the solutions of the
above recurrence equation (2.1) as a function of the hyper-parameters and the
feedforward and feedback connections matrices $\mathcal{W}^{f}$ and
$\mathcal{W}^{b}$. We would like to stay as general as possible to encompass
as many situations as possible, keeping in mind that we already made strong
assumptions by imposing that the weight matrices of feedforward and feedback
connections are identical from one layer to another. Motivated by concrete
applications, we will mainly consider matrices $\mathcal{W}^{f}$ and
$\mathcal{W}^{b}$ which act as convolutions on $\mathbb{R}^{d}$.

$n+1\quad$$n\quad$$j-1$$j$$j+1$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$
Figure 1: Schematic illustration of the network structure of model (2.1) where
each point represents a given neuronal layer index $j$ (x-axis) at a
particular time step $n$ (y-axis), and the red arrows indicate the
contributions leading to the update of $\mathcal{E}_{j}^{n+1}$.

##  3 The identity case

It turns out that we will gain much information by first treating the
simplified case where $\mathcal{W}^{f}$ and $\mathcal{W}^{b}$ are both
identity. That is, from now on, and throughout this section we assume that

$\mathcal{W}^{f}=\mathcal{W}^{b}=\mathbf{I}_{d}.$

That is, each neuron in a layer is only connected to the corresponding neuron
in the immediately preceding and following layer, with unit weight in each
direction. Under such a setting, the recurrence equation (2.1) reduces to a
scalar equation, that is

$e_{j}^{n+1}=\beta e_{j-1}^{n+1}+\alpha
e_{j-1}^{n}+(1-\beta-\lambda-\alpha)e_{j}^{n}+\lambda e_{j+1}^{n},\quad
j=1,\cdots,J-1,$ (3.1)

with this time the unknown $e_{j}^{n}\in\mathbb{R}$, together with

$e_{0}^{n}=s_{0}^{n},\quad n\geq 0,$ (3.2)

and

$e_{J}^{n+1}=\beta e_{J-1}^{n+1}+\alpha
e_{J-1}^{n}+(1-\beta-\alpha)e_{J}^{n},\quad n\geq 0.$ (3.3)

with

$e^{0}_{j}=h_{j},\quad j=0,\cdots,J.$ (3.4)

###  3.1 Wave propagation on an infinite depth network

It will be first useful to consider the above problem set on an infinite
domain and look at

$e_{j}^{n+1}=\beta e_{j-1}^{n+1}+\alpha
e_{j-1}^{n}+(1-\beta-\lambda-\alpha)e_{j}^{n}+\lambda e_{j+1}^{n},\quad
j\in\mathbb{Z},$ (3.5)

given some initial sequence

$e^{0}_{j}=h_{j},\quad j\in\mathbb{Z}.$

This situation has no direct equivalent in the brain, where the number of
hierarchically connected layers is necessarily finite; but it is a useful
mathematical construct. Indeed, such recurrence equations set on the integers
$\mathbb{Z}$ are relatively well understood from the mathematical numerical
analysis community. The behavior of the solution sequence
$(e_{j}^{n})_{j\in\mathbb{Z}}$ can be read out from the so-called
amplification factor function defined as

$\rho(\theta):=\frac{\alpha\left(e^{-\mathbf{i}\theta}-1\right)+1-\beta+\lambda\left(e^{\mathbf{i}\theta}-1\right)}{1-\beta
e^{-\mathbf{i}\theta}},\quad\theta\in[-\pi,\pi],$ (3.6)

and which relates spatial and temporal modes. Indeed, formally, the sequence
$(\rho(\theta)^{n}e^{\mathbf{i}j\theta})_{j\in\mathbb{Z}}$ is an explicit
solution to (3.5) for each $\theta\in[-\pi,\pi]$. Actually one can be much
precise and almost explicit in the sense that one can relate the expression of
the solutions to (3.5) starting from some initial sequence
$(h_{j})_{j\in\mathbb{Z}}$ to the properties of $\rho$ in a systematic way
that we now briefly explain.

Let us first denote by
$\mathcal{G}^{n}=(\mathcal{G}^{n}_{j})_{j\in\mathbb{Z}}$ the sequence which is
the fundamental solution of (3.5) in the special case where
$(\mathcal{H}_{j})_{j\in\mathbb{Z}}$ is the Dirac delta sequence
$\boldsymbol{\delta}$. The Dirac delta sequence $\boldsymbol{\delta}$ is
defined as $\boldsymbol{\delta}_{0}=1$ and $\boldsymbol{\delta}_{j}=0$ for all
$j\in\mathbb{Z}\backslash\\{0\\}$. As a consequence, we have
$\mathcal{G}^{0}=\boldsymbol{\delta}$ and for each $n\geq 0$

$\mathcal{G}_{j}^{n+1}-\beta\mathcal{G}_{j-1}^{n+1}=\alpha\mathcal{G}_{j-1}^{n}+(1-\beta-\lambda-\alpha)\mathcal{G}_{j}^{n}+\lambda\mathcal{G}_{j+1}^{n},\quad
j\in\mathbb{Z}.$

The starting point of the analysis is the following representation formula,
obtained via inverse Fourier transform, which reads

$\mathcal{G}^{n}_{j}=\frac{1}{2\pi}\int_{-\pi}^{\pi}e^{\mathbf{i}j\theta}\rho(\theta)^{n}\mathrm{d}\theta,\quad
n\geq 1,\quad j\in\mathbb{Z}.$ (3.7)

Then, given any initial sequence $(h_{j})_{j\in\mathbb{Z}}$, the solution
$(e_{j}^{n})_{j\in\mathbb{Z}}$ to (3.5) can be represented as the convolution
product between the initial sequence and the fundamental solution, namely

$e_{j}^{n}=\sum_{\ell\in\mathbb{Z}}\mathcal{G}_{j-\ell}^{n}h_{\ell},\quad
j\in\mathbb{Z},\quad n\geq 1.$ (3.8)

That is, having characterized the fundamental solution for a simple input
pattern ($\boldsymbol{\delta}$), with a unitary impulse provided to a single
layer, we can now easily generalize to any arbitrary input pattern, by
applying the (translated) fundamental solution to each layer.

Our aim is to understand under which conditions on the hyper-parameters we can
ensure that the solutions of (3.5) given through (3.8) remain bounded for all
$n\geq 1$ independently of the choice of the initial sequence
$(h_{j})_{j\in\mathbb{Z}}$. More precisely, we introduce the following
terminology. We say that the recurrence equation is stable if for each bounded
initial sequence $(h_{j})_{j\in\mathbb{Z}}\in\ell^{\infty}(\mathbb{Z})$, the
corresponding solution $(e_{j}^{n})_{j\in\mathbb{Z}}$ given by (3.8) satisfies

$\underset{j\in\mathbb{Z}}{\sup}|e_{j}^{n}|\underset{n\rightarrow\infty}{\longrightarrow}0.$

On the other hand, we say that the recurrence equation is unstable if one can
find a bounded initial sequence
$(h_{j})_{j\in\mathbb{Z}}\in\ell^{\infty}(\mathbb{Z})$ such that the
corresponding solution $(e_{j}^{n})_{j\in\mathbb{Z}}$ given by (3.8) satisfies

$\underset{j\in\mathbb{Z}}{\sup}|e_{j}^{n}|\underset{n\rightarrow\infty}{\longrightarrow}+\infty.$

Finally, we say that the recurrence equation is marginally stable if there
exists a universal constant $C>0$ such that for each bounded initial sequence
$(h_{j})_{j\in\mathbb{Z}}\in\ell^{\infty}(\mathbb{Z})$, the corresponding
solution $(e_{j}^{n})_{j\in\mathbb{Z}}$ given by (3.8) satisfies

$\underset{j\in\mathbb{Z}}{\sup}|e_{j}^{n}|\leq
C~{}\underset{j\in\mathbb{Z}}{\sup}|h_{j}|,\quad n\geq 1.$

It turns out that one can determine the stability properties of the recurrence
equation by solely looking at the amplification factor function. Indeed, from
[29], we know that

$\underset{n\rightarrow\infty}{\lim}\|\mathcal{G}^{n}\|_{\ell^{1}(\mathbb{Z})}^{1/n}=\underset{\theta\in[-\pi,\pi]}{\max}|\rho(\theta)|,$

where we have set

$\|\mathcal{G}^{n}\|_{\ell^{1}(\mathbb{Z})}:=\sum_{j\in\mathbb{Z}}|\mathcal{G}_{j}^{n}|.$

As a consequence, we directly deduce that the recurrence equation is stable
when $|\rho(\theta)|<1$ for all $\theta\in[-\pi,\pi]$, whereas it is unstable
if there exists $\theta_{0}\in[-\pi,\pi]$ such that $|\rho(\theta_{0})|>1$.
The limiting case occurs precisely when
$\underset{\theta\in[-\pi,\pi]}{\max}|\rho(\theta)|=1$ and there is actually a
long history of works [30, 14, 26, 12, 10] that have studied the marginal
stability of the recurrence equation in that case. All such results rely on a
very precise understanding of the amplification factor function and lead to
the following statement.

######  Theorem 1 ([30, 14, 26, 12, 10]).

Suppose that there exist finitely many
$\theta_{1},\cdots,\theta_{K}\in[-\pi,\pi]$ such that for all
$\theta\in[-\pi,\pi]\backslash\left\\{\theta_{1},\cdots,\theta_{k}\right\\}$
one has $|\rho(\theta)|<1$ and $|\rho(\theta_{k})|=1$ for each $k=1,\cdots,K$.
Furthermore, assume that there exist $c_{k}\in\mathbb{R}$,
$\sigma_{k}\in\mathbb{C}$ with $\mathrm{Re}(\sigma_{k})>0$ and an integer
$\mu_{k}\geq 1$ such that

$\frac{\rho(\theta_{k}+\theta)}{\rho(\theta_{k})}=\exp\left(-\mathbf{i}c_{k}\theta-\sigma_{k}\theta^{2\mu_{k}}+\mathcal{O}(|\theta|^{2\mu_{k}+1})\right),\text{
as }\theta\rightarrow 0.$

Then the recurrence equation is marginally stable.

Based on the above notions of stability/instability, we see that the only
interesting situation is when the recurrence equation is marginally stable,
and thus when the amplification function is contained in the unit disk with
finitely many tangent points to the unit circle with prescribed asymptotic
expansions. This is also the only interesting situation from a biological
standpoint, as it ensures that the network remains active, yet without runaway
activations.

####  3.1.1 Study of the amplification factor function

Since we assumed that $0\leq\beta<1$, the denominator in (3.6) never vanishes
and is well-defined. Next, we crucially remark that we always have

$\rho(0)=1.$

We will now check under which conditions $|\rho(\theta)|\leq 1$ for all
$\theta\in[-\pi,\pi]$ to guarantee marginal stability of the recurrence
equation.

To assess stability, we compute

$\displaystyle\left|\rho(\theta)\right|^{2}$
$\displaystyle=\frac{\left((\lambda+\alpha)(\cos(\theta)-1)+1-\beta\right)^{2}+(\lambda-\alpha)^{2}\sin(\theta)^{2}}{1-2\beta\cos(\theta)+\beta^{2}}$
$\displaystyle=\frac{(\lambda+\alpha)^{2}(\cos(\theta)-1)^{2}+2(1-\beta)(\lambda+\alpha)(\cos(\theta)-1)+\left(1-\beta\right)^{2}+(\lambda-\alpha)^{2}(1-\cos(\theta)^{2})}{(1-\beta)^{2}+2\beta(1-\cos(\theta))}$
$\displaystyle=\frac{(1-\cos(\theta))\left((\lambda+\alpha)^{2}(1-\cos(\theta))-2(1-\beta)(\lambda+\alpha)+(\lambda-\alpha)^{2}(1+\cos(\theta))\right)+\left(1-\beta\right)^{2}}{(1-\beta)^{2}+2\beta(1-\cos(\theta))}$
$\displaystyle=\frac{(1-\cos(\theta))\left(-4\alpha\lambda\cos(\theta)-2(1-\beta)(\lambda+\alpha)+2(\lambda^{2}+\alpha^{2})\right)+\left(1-\beta\right)^{2}}{(1-\beta)^{2}+2\beta(1-\cos(\theta))}$

such that $\left|\rho(\theta)\right|^{2}\leq 1$ is equivalent to

$(1-\cos(\theta))\left(2\beta+4\alpha\lambda\cos(\theta)+2(1-\beta)(\lambda+\alpha)-2(\lambda^{2}+\alpha^{2})\right)\geq
0,\quad\theta\in[-\pi,\pi],$

and since $1-\cos(\theta)\geq 0$ we need to ensure

$\beta+2\alpha\lambda\cos(\theta)+(1-\beta)(\lambda+\alpha)-\lambda^{2}-\alpha^{2}\geq
0,\quad\theta\in[-\pi,\pi],$

and evaluating at $\pm\pi$ the above inequality we get

$\beta+(1-\beta)(\lambda+\alpha)-(\lambda+\alpha)^{2}\geq 0.$

But we remark that the above expression can be factored as

$\left(\beta+\lambda+\alpha\right)\left(1-\lambda-\alpha\right)\geq 0.$

As a consequence, $\left|\rho(\theta)\right|^{2}\leq 1$ if and only if
$\lambda+\alpha\leq 1$. This is precisely the condition that we made in (2.6).
We can actually track cases of equality which are those values of
$\theta\in[-\pi,\pi]$ for which we have

$(1-\cos(\theta))\left(2\beta+4\alpha\lambda\cos(\theta)+2(1-\beta)(\lambda+\alpha)-2(\lambda^{2}+\alpha^{2})\right)=0.$

We readily recover that at $\theta=0$ we have $|\rho(0)|=1$. So, now assuming
that $\theta\neq 0$, we need to solve

$\beta+2\alpha\lambda\cos(\theta)+(1-\beta)(\lambda+\alpha)-\lambda^{2}-\alpha^{2}=0,$

which we write as

$\beta-2\alpha\lambda+(1-\beta)(\lambda+\alpha)-\lambda^{2}-\alpha^{2}+2\alpha\lambda\left(\cos(\theta)+1)\right)=0,$

and using the previous factorization we get

$\left(\beta+\lambda+\alpha\right)\left(1-\lambda-\alpha\right)+2\alpha\lambda\left(\cos(\theta)+1)\right)=0,$

and we necessarily get that both $1+\cos(\theta)=0$ and $1-\lambda-\alpha=0$
must be satisfied. As consequence, $|\rho(\pm\pi)|=1$ if and only if
$1=\lambda+\alpha$.

As a summary we have obtained that:

  * •

if $0\leq\lambda+\alpha<1$ and $0\leq\beta<1$, then $|\rho(\theta)|<1$ for all
$\theta\in[-\pi,\pi]\backslash\\{0\\}$ with $\rho(0)=1$;

  * •

if $\lambda+\alpha=1$ and $0\leq\beta<1$, then $|\rho(\theta)|<1$ for all
$\theta\in(-\pi,\pi)\backslash\\{0\\}$ with $\rho(0)=1$ and $\rho(\pm\pi)=-1$.

We present in Figure 2 several representative illustrations of the spectral
curves $\rho(\theta)$ for various values of the hyper-parameters recovering
the results explained above.

(a) $\beta=0$ and $\alpha+\lambda<1$. (b) $\beta>0$ and $\alpha+\lambda<1$.
(c) $\beta=0$ and $\alpha+\lambda=1$. (d) $\beta>0$ and $\alpha+\lambda=1$.

Figure 2: Several representative illustration of the curve
$\theta\mapsto\rho(\theta)$ for $\theta\in[-\pi,\pi]$ in the case where
$\beta=0$ and $\beta\neq 0$. (a)-(b) Amplification factor function
$\rho(\theta)$ (blue curve) with a unique tangency point on the unit circle at
$z=1$ corresponding $\theta=0$. (c)-(f) When $\alpha+\lambda=1$ the function
$\rho(\theta)$ (blue curve) has two tangency points on the unit circle at
$z=1$ corresponding $\theta=0$ and $z=-1$ corresponding to $\theta=\pm\pi$.

Furthermore, near $\theta\sim 0$, we get that $\rho$ admits the following
asymptotic expansion

$\rho(\theta)=\exp\left(-\mathbf{i}\frac{\beta+\alpha-\lambda}{1-\beta}\theta-\frac{\beta(1-\alpha-\lambda)+\alpha+\lambda-(\lambda-\alpha)^{2}}{2(1-\beta)^{2}}\theta^{2}+\mathcal{O}(|\theta|^{3})\right),\text{
as }\theta\rightarrow 0,$

provided that

$\beta(1-\alpha-\lambda)+\alpha+\lambda-(\lambda-\alpha)^{2}\neq 0.$

In fact, since $-(\lambda-\alpha)^{2}\geq-(\lambda+\alpha)^{2}$ as both
$\alpha$ and $\lambda$ are positive, we remark that

$\beta(1-\alpha-\lambda)+\alpha+\lambda-(\lambda-\alpha)^{2}\geq\beta(1-\alpha-\lambda)+\alpha+\lambda-(\lambda+\alpha)^{2}=\left(\beta+\lambda+\alpha\right)\left(1-\lambda-\alpha\right)\geq
0.$

Finally, we remark that when $\alpha+\lambda=1$ we have

$\rho(\theta+\pi)=-\exp\left(-\mathbf{i}\frac{-\beta+\alpha-\lambda}{1+\beta}\theta-\frac{1-(\alpha-\lambda)^{2}}{2(1+\beta)^{2}}\theta^{2}+\mathcal{O}(|\theta|^{3})\right),\text{
as }\theta\rightarrow 0.$

From now on, we denote

$\displaystyle c_{0}:=\frac{\beta+\alpha-\lambda}{1-\beta},$
$\displaystyle\quad\sigma_{0}:=\frac{\beta(1-\alpha-\lambda)+\alpha+\lambda-(\lambda-\alpha)^{2}}{2(1-\beta)^{2}},$
$\displaystyle c_{\pi}:=\frac{-\beta+\alpha-\lambda}{1+\beta},$
$\displaystyle\quad\sigma_{\pi}:=\frac{1-(\alpha-\lambda)^{2}}{2(1+\beta)^{2}},$

and we always assume that

$\sigma_{0}>0,\text{ and }\sigma_{\pi}>0,$

which is equivalent to assume that $0<\alpha<1$ and $0<\lambda<1$.

Here, $(c_{0},\sigma_{0})$ and $(c_{\pi},\sigma_{\pi})$ are derived,
respectively, from the asymptotic expansions of the amplification factor
function $\rho(\theta)$ near $\theta=0$ and $\theta=\pi$, as defined above. On
the one hand $c_{0}$ reflects the propagation speed of the solution associated
with $\rho(0)$, while $\sigma_{0}$ can be understood as its spatio-temporal
spread (and similarly for the solution potentially associated with
$\rho(\pi)$). In the following, we explore the fundamental solutions of this
system for various values of its hyper-parameters.

####  3.1.2 Turning off the instantaneous feedforward connections: case
$\beta=0$

We first investigate the case where there is no instantaneous feedforward
connections in the network, that is we set $\beta=0$. This case, although less
generic, is compatible with the prominent Rao-Ballard formulation of
predictive coding [27], in which feedforward connections—after contributing to
setting the initial network activity—only convey prediction errors, as
captured by the hyper-parameter $\alpha$. In that case, the model is fully
explicit: the update at time step $n+1$ only depends on the internal states at
the previous step $n$ since we simply have

$e_{j}^{n+1}=\alpha e_{j-1}^{n}+(1-\lambda-\alpha)e_{j}^{n}+\lambda
e_{j+1}^{n},\quad j\in\mathbb{Z}.$

As we assumed that $\alpha+\lambda\leq 1$, the right-hand side of the
recurrence equation is a positive linear combination of elements of the
sequence $(e_{j}^{n})$ such that we have positivity principle of the solution,
namely

$\forall j\in\mathbb{Z},\quad 0\leq h_{j}\quad\Longrightarrow\quad\forall
j\in\mathbb{Z},\quad n\geq 1,\quad 0\leq e_{j}^{n}.$

Furthermore, since the recurrence equation is explicit, we have finite speed
propagation, in the following sense. Recall that when $\beta=0$, the
fundamental solution $\mathcal{G}^{n}$ is solution to

$\mathcal{G}_{j}^{n+1}=\alpha\mathcal{G}_{j-1}^{n}+(1-\lambda-\alpha)\mathcal{G}_{j}^{n}+\lambda\mathcal{G}_{j+1}^{n},\quad
n\geq 1,\quad\quad j\in\mathbb{Z},$

starting from $\mathcal{G}^{0}=\boldsymbol{\delta}$. Finite speed of
propagation then refers to the property that

$\mathcal{G}_{j}^{n}=0,\quad|j|>n.$

This in turn implies that necessarily $c_{0}\in(-1,1)$ which is readily seen
from the explicit formula $c_{0}=\alpha-\lambda$ in that case. Actually, it is
possible to be more precise and to give a general expression for the
fundamental solution. Roughly speaking, each $\mathcal{G}_{j}^{n}$ ressembles
a discrete Gaussian distribution centered at $j=c_{0}n$ and we refer to the
recent theoretical results of [14, 26, 10, 12] for a rigorous justification.

Essentially, the results can be divided into two cases depending on whether or
not $\alpha+\lambda=1$. As can be seen above, the special case
$\alpha+\lambda=1$ results in a cancellation of the “memory” term, such that a
neuronal layer $j$’s activity does not depend on its own activity at the
previous time step, but only on the activity of its immediate neighbors $j-1$
and $j+1$. More precisely, we have the following:

(a) $\beta=0$, $\lambda<\alpha$ and $\alpha+\lambda<1$. (b) $\beta=0$,
$\alpha<\lambda$ and $\alpha+\lambda=1$.

Figure 3: Illustration of the evolution of the fundamental solution
$\mathcal{G}^{n}$ starting from the Dirac delta sequence at $j=0$ in the case
$\beta=0$ at several time iterations. The circles represent the numerically
computed solution while the plain lines represent the Gaussian approximation.
(a) When $\lambda<\alpha$ there is a rightward propagation along a Gaussian
profile whose leading profile is given by
$\frac{1}{\sqrt{4\pi\sigma_{0}n}}\exp\left(-\frac{|j-c_{0}n|^{2}}{4\sigma_{0}n}\right)$.
(b) When $\alpha<\lambda$ there is a leftward propagation along a Gaussian
profile whose leading profile is given by
$\frac{1+(-1)^{n+j}}{\sqrt{4\pi\sigma_{0}n}}\exp\left(-\frac{|j-c_{0}n|^{2}}{4\sigma_{0}n}\right)$
which vanishes whenever $n+j$ is odd.

  * •

Case: $0\leq\lambda+\alpha<1$. The fundamental solution can be decomposed as

$\mathcal{G}_{j}^{n}=\frac{1}{\sqrt{4\pi\sigma_{0}n}}\exp\left(-\frac{|j-c_{0}n|^{2}}{4\sigma_{0}n}\right)+\mathcal{N}_{j}^{n},\quad|j|\leq
n,$

where the remainder term satisfies an estimate

$\left|\mathcal{N}_{j}^{n}\right|\leq\frac{C}{n}\exp\left(-\kappa\frac{|j-c_{0}n|^{2}}{n}\right),\quad|j|\leq
n,$

for some universal constants $C,\kappa>0$ which only depend on the hyper-
parameters and not $n$ and $j$. In Figure 3(a), we represented the fundamental
solution $\mathcal{G}_{j}^{n}$ at different time iterations (circles) in the
case $\lambda<\alpha$ where there is rightward propagation with $c_{0}>0$ and
compared it with the leading order fixed Gaussian profile centered at
$j=c_{0}n$ (plain line). On the other hand, in Figure 4, panels (a)-(b)-(c),
we illustrate the above results by presenting a space-time color plot of the
fundamental solution rescaled by a factor $\sqrt{n}$. We observe rightward
(respectively leftward) propagation with $c_{0}>0$ (respectively $c_{0}<0$)
when $\lambda<\alpha$ (respectively $\alpha<\lambda$), while when
$\alpha=\lambda$ we have $c_{0}=0$ and no propagation occurs.

  * •

Case: $\lambda+\alpha=1$. In this case, we first note that we have
$c_{0}=c_{\pi}$ together with $\sigma_{0}=\sigma_{\pi}$ and

$\mathcal{G}_{j}^{n}=\frac{1+(-1)^{n+j}}{\sqrt{4\pi\sigma_{0}n}}\exp\left(-\frac{|j-c_{0}n|^{2}}{4\sigma_{0}n}\right)+\mathcal{N}_{j}^{n},\quad|j|\leq
n,$

where the remainder term satisfies an estimate

$\left|\mathcal{N}_{j}^{n}\right|\leq\frac{C}{n}\exp\left(-\kappa\frac{|j-c_{0}n|^{2}}{n}\right),\quad|j|\leq
n,$

for some universal constants $C,\kappa>0$. In Figure 3(b), we represented the
fundamental solution $\mathcal{G}_{j}^{n}$ at different time iterations
(circles) in the case $\alpha<\lambda$ where there is leftward propagation
with $c_{0}<0$ and compared it with the leading order fixed Gaussian profile
centered at $j=c_{0}n$ (plain line). Similarly as in the previous case, in
Figure 4, panels (d)-(e)-(f), we illustrate the above results by presenting a
space-time color plot of the fundamental solution rescaled by a factor
$\sqrt{n}$. The direction of propagation still depends on the sign of $c_{0}$
and whether or not $\lambda\lessgtr\alpha$. Unlike the case
$\alpha+\lambda<1$, we observe a tiled pattern where $\mathcal{G}_{j}^{n}=0$
for even or odd integers alternatively for each time step.

(a) $\beta=0$, $\lambda<\alpha$ and $\alpha+\lambda<1$. (b) $\beta=0$,
$\lambda=\alpha$ and $\alpha+\lambda<1$. (c) $\beta=0$, $\alpha<\lambda$ and
$\alpha+\lambda<1$. (d) $\beta=0$, $\lambda<\alpha$ and $\alpha+\lambda=1$.
(e) $\beta=0$, $\lambda=\alpha$ and $\alpha+\lambda=1$. (f) $\beta=0$,
$\alpha<\lambda$ and $\alpha+\lambda=1$.

Figure 4: Illustration of the evolution of the rescaled solution sequence
$(\sqrt{n}\,\mathcal{G}_{j}^{n})_{j\in\mathbb{Z}}$ starting from the Dirac
delta sequence at $j=0$ in the case $\beta=0$. First row: $\alpha+\lambda<1$
and second row: $\alpha+\lambda=1$. When $\lambda\lessgtr\alpha$, we observe a
rightward/leftward propagation while when $\alpha=\lambda$ no propagation
occurs. In all panels, the pink curve is given by $j=nc_{0}$, clearly
illustrating the fact that $c_{0}$ measures the propagation speed of the
solution. Note that in the case $\beta=0$ and $\alpha+\lambda=1$, we have
$c_{0}=c_{\pi}$ which results in the tiled patterns observed in panels
(d)-(e)-(f).

As a partial intermediate summary, we note that the sign of $c_{0}$ (directly
related to the sign of $\alpha-\lambda$) always indicates in which direction
the associated Gaussian profile propagates. Namely if $\alpha>\lambda$ and
$c_{0}>0$ (resp. $\alpha<\lambda$ and $c_{0}<0$) there is rightward (resp.
leftward) propagation. Intuitively, this behavior reflects the functional role
of each hyper-parameter, with $\alpha$ and $\lambda$ controlling feed-forward
and feed-back prediction error correction, respectively. When
$\alpha=\lambda$, the two terms are equally strong, and there is no dominant
direction of propagation. In addition, when $\lambda+\alpha=1$, the Gaussian
profile is _oscillating_ because of the presence of $(-1)^{n+j}$. As will be
seen later when considering continuous versions of our model, this oscillatory
pattern might not be truly related to neural oscillations observed in the
brain, but could instead arise here as a consequence of discrete updating.

Finally, we note that the fundamental solution sequence
$(\mathcal{G}^{n}_{j})_{j\in\mathbb{Z}}$ is uniformly integrable for all
values of the parameters, that is there exists some universal constant $C>0$,
depending only on the hyper-parameters such that

$\|\mathcal{G}^{n}\|_{\ell^{1}(\mathbb{Z})}:=\sum_{j\in\mathbb{Z}}|\mathcal{G}_{j}^{n}|\leq
C,\quad n\geq 1.$

As a consequence, since given any bounded initial sequence
$(h_{j})_{j\in\mathbb{Z}}\in\ell^{\infty}(\mathbb{Z})$, the solution
$(e_{j}^{n})_{j\in\mathbb{Z}}$ to (3.5) can be represented as the convolution
product between the initial sequence and the fundamental solution, namely

$e_{j}^{n}=\sum_{\ell\in\mathbb{Z}}\mathcal{G}_{j-\ell}^{n}h_{\ell},\quad
j\in\mathbb{Z},\quad n\geq 1,$

we readily deduce that the solution $(e_{j}^{n})_{j\in\mathbb{Z}}$ is
uniformly bounded with respect to $n$, that is there exists some universal
constant denoted $C>0$, such that

$\underset{j\in\mathbb{Z}}{\sup}\left|e_{j}^{n}\right|\leq
C~{}\underset{j\in\mathbb{Z}}{\sup}\left|h_{j}\right|,\quad n\geq 1.$

This is exactly our definition of marginal stability.

####  3.1.3 Turning on the instantaneous feedforward connections: case
$\beta>0$

We now turn to the general case where $\beta>0$. That is, the feed-forward
connections continue to convey sensory inputs at each time step following the
network initializing, and $\beta$ controls the strength of these signals. In
that case, the recurrence equation is no longer explicit but implicit and the
positivity property together with the finite speed propagation no longer hold
true in general. Indeed, upon introducing the shift operator

$\mathbf{S}:(u_{j})_{j\in\mathbb{Z}}\mapsto(u_{j+1})_{j\in\mathbb{Z}},$

we remark that equation (3.5) can be written

$\left(\mathrm{Id}-\beta\mathbf{S}^{-1}\right)e^{n+1}=\alpha\mathbf{S}^{-1}e^{n}+(1-\beta-\lambda-\alpha)e^{n}+\lambda\mathbf{S}e^{n},\quad
n\geq 0,$

with $e^{n}=(e_{j}^{n})_{j\in\mathbb{Z}}$. Since $0<\beta<1$ and
$\||\mathbf{S}^{-1}\||_{\ell^{q}(\mathbb{Z})\to\ell^{q}(\mathbb{Z})}=1$ for
any $q\in[1,+\infty]$, the operator $\mathrm{Id}-\beta\mathbf{S}^{-1}$ is
invertible on $\ell^{q}(\mathbb{Z})$ for any $q\in[1,+\infty]$ with inverse

$\left(\mathrm{Id}-\beta\mathbf{S}^{-1}\right)^{-1}=\sum_{\ell=0}^{\infty}\beta^{\ell}\mathbf{S}^{-\ell}.$

As a consequence, the recurrence equation can be recast as a convolution
operator across the network layers with infinite support, namely

$e^{n+1}_{j}=\alpha\sum_{\ell=0}^{\infty}\beta^{\ell}e^{n}_{j-\ell-1}+(1-\beta-\lambda-\alpha)\sum_{\ell=0}^{\infty}\beta^{\ell}e^{n}_{j-\ell}+\lambda\sum_{\ell=0}^{\infty}\beta^{\ell}e^{n}_{j-\ell+1},\quad
j\in\mathbb{Z},\quad n\geq 0.$

From the above expression, we readily deduce that the positivity of the
solution is preserved whenever $0<\beta<1-\lambda-\alpha$. Furthermore, for
the fundamental solution starting from the Dirac delta solution which solves

$\mathcal{G}_{j}^{n+1}-\beta\mathcal{G}_{j-1}^{n+1}=\alpha\mathcal{G}_{j-1}^{n}+(1-\beta-\lambda-\alpha)\mathcal{G}_{j}^{n}+\lambda\mathcal{G}_{j+1}^{n},\quad
j\in\mathbb{Z},\quad n\geq 0,$

we only have that

$\mathcal{G}_{j}^{n}=0,\quad j<-n,$

which implies that $-1<c_{0},c_{\pi}<+\infty$. Indeed, from the formula of
$c_{0}$ we get that

$c_{0}\sim\frac{1-\alpha-\beta}{1-\beta}\longrightarrow+\infty,\quad\text{ as
}\beta\rightarrow 1^{-}.$

Once again, as in the case with $\beta=0$, we can precise the behavior of the
fundamental solution by using the combined results of [10, 12].

(a) $\beta+\alpha<\lambda$. (b) $\beta+\alpha=\lambda$. (c)
$\lambda<\beta+\alpha$.

Figure 5: Effects of turning on $\beta>0$ when $\alpha+\lambda<1$. When
$\beta+\alpha<\lambda$ we have $c_{0}<0$ and observe backward propagation in
panel (a) while when $\lambda<\beta+\alpha$ we have $c_{0}>0$ and have forward
propagation as seen in panel (c). At the transition $\beta+\alpha=\lambda$,
the wave speed vanishes $c_{0}=0$ and there is no propagation as illustrated
in panel (b).Intuitively, $\alpha$ and $\beta$, both propagating signals in
the forward (rightward) direction, compete with $\lambda$ carrying the
feedback (leftward) prediction signals; this competition determines the main
direction of propagation of neural activity in the system.

  * •

Case: $0\leq\lambda+\alpha<1$. There exist some universal constants
$C,\kappa>0$ and $L>0$ such that

$\mathcal{G}_{j}^{n}=\frac{1}{\sqrt{4\pi\sigma_{0}n}}\exp\left(-\frac{|j-c_{0}n|^{2}}{4\sigma_{0}n}\right)+\mathcal{N}_{j}^{n},\quad-n\leq
j\leq Ln,$

where the remainder term satisfies a Gaussian estimate

$\left|\mathcal{N}_{j}^{n}\right|\leq\frac{C}{n}\exp\left(-\kappa\frac{|j-c_{0}n|^{2}}{n}\right),\quad-n\leq
j\leq Ln.$

While for $j>nL$ we simply get a pure exponential bound

$\left|\mathcal{G}_{j}^{n}\right|\leq Ce^{-\kappa n-\kappa j},\quad nL<j.$

Inspecting the formula for $c_{0}$, we notice that when
$\alpha+\beta\lessgtr\lambda$ we have $c_{0}\lessgtr 0$ and the wave speed
vanishes precisely when $\alpha+\beta=\lambda$. This is illustrated in Figure
5, where we see that $\alpha$ and $\beta$, both propagating signals in the
forward (rightward) direction, compete with $\lambda$ carrying the feedback
(leftward) prediction signals; this competition determines the main direction
of propagation of neural activity in the system.

(a) $\alpha<\beta+\lambda$ and $\beta+\alpha<\lambda$. (b)
$\alpha<\beta+\lambda$ and $\beta+\alpha=\lambda$. (c) $\alpha<\beta+\lambda$
and $\lambda<\beta+\alpha$. (d) $\beta+\lambda<\alpha$. (e)
$\beta+\lambda=\alpha$.

Figure 6: Effects of turning on $\beta>0$ when $\alpha+\lambda=1$. We now
observe a secondary wave with associated wave speed $c_{\pi}$ whose sign
depends on the competition between $\alpha$ and $\beta+\lambda$. (a)-(b)-(c)
When $\alpha<\beta+\lambda$, the wave speed of the secondary wave always
verifies $c_{\pi}<0$, and the competition between $\lambda$ and $\beta+\alpha$
gives the direction of the primary wave as previously reported in Figure 5.
(d) When $\beta+\lambda<\alpha$ which always implies that
$\alpha+\beta>\lambda$, we have $0<c_{\pi}<c_{0}$ traducing forward
propagation for both waves. We remark that the secondary wave is slower. (e)
When $\beta+\lambda=\alpha$ which also implies that $\alpha+\beta>\lambda$, we
get $0=c_{\pi}<c_{0}$ such that the secondary wave is blocked. (f) Summary of
the sign of the wave speeds $c_{0}$ and $c_{\pi}$ when $\beta>0$ when
$\alpha+\lambda\leq 1$.

  * •

Case: $\lambda+\alpha=1$. What changes in that case is the existence of a
secondary wave with associated wave speed $c_{\pi}$ whose sign depends on the
competition between $\alpha$ and $\beta+\lambda$. When $\alpha<\beta+\lambda$
then we have $c_{\pi}<0$, and the competition between $\lambda$ and
$\beta+\alpha$ will determine the sign of $c_{0}$, as illustrated in panels
(a)-(b)-(c) of Figure 6. On the other hand, when $\beta+\lambda<\alpha$
implying that $c_{\pi}>0$, we note that $\alpha+\beta>\lambda$ and thus
$c_{0}>0$. In that case, the explicit formula for $c_{\pi}$ and $c_{0}$ shows
that $0<c_{\pi}<c_{0}$ and the secondary wave associated to $c_{\pi}$ is
slower to propagate into the network, see Figure 6(d). Finally, when
$\beta+\lambda=\alpha$ we have $0=c_{\pi}<c_{0}$ and the secondary wave is
blocked, see Figure 6(e).

We have summarized in the diagram of Figure 6(f) all possible configurations
for the sign of the wave speeds $c_{0}$ and $c_{\pi}$ when $\beta\in(0,1)$
when $\alpha+\lambda\leq 1$. We notably observe that when $\beta$ is increased
the region of parameter space where $c_{0}<0$ diminishes while the region of
parameter space where $c_{\pi}<0$ increases, indicating that for high values
of $\beta$ the primary wave is most likely to be forward while the secondary
wave is most likely to be backward.

###  3.2 Wave propagation on a semi-infinite network with a forcing source
term

Now that we have understood the intrinsic underlying mechanisms of wave
propagation for our model (3.1) set on an infinite domain, we turn to the case
where the network is semi-infinite. That is, the network admits an _input_
layer that is only connected to the layer above. The problem now reads

$\left\\{\begin{split}e_{j}^{n+1}-\beta e_{j-1}^{n+1}&=\alpha
e_{j-1}^{n}+(1-\beta-\lambda-\alpha)e_{j}^{n}+\lambda e_{j+1}^{n},\quad j\geq
1,\quad n\geq 0,\\\ e_{0}^{n}&=s_{0}^{n},\quad n\geq 0,\\\
e_{j}^{0}&=h_{j},\quad j\geq 1.\end{split}\right.$ (3.9)

We see that the system depends on the source term $s_{0}^{n}$ applied to its
input layer at each time step, also called a _boundary value_ , and on the
starting activation value $(h_{j})$ applied to each layer at the initial time
point, also called the _initial value_. In fact, the linearity principle tells
us that the solutions of the above problem can be obtained as the linear
superposition of the solutions to the following two problems, the boundary
value problem, where all layers except the input layer are initialized at
zero:

$\left\\{\begin{split}g_{j}^{n+1}-\beta g_{j-1}^{n+1}&=\alpha
g_{j-1}^{n}+(1-\beta-\lambda-\alpha)g_{j}^{n}+\lambda g_{j+1}^{n},\quad j\geq
1,\quad n\geq 0,\\\ g_{0}^{n}&=s_{0}^{n},\quad n\geq 0,\\\ g_{j}^{0}&=0,\quad
j\geq 1,\end{split}\right.$ (3.10)

and the initial value problem, where the input layer source term is set to
zero for all time steps:

$\left\\{\begin{split}f_{j}^{n+1}-\beta f_{j-1}^{n+1}&=\alpha
f_{j-1}^{n}+(1-\beta-\lambda-\alpha)f_{j}^{n}+\lambda f_{j+1}^{n},\quad j\geq
1,\quad n\geq 0,\\\ f_{0}^{n}&=0,\quad n\geq 0,\\\ f_{j}^{0}&=h_{j},\quad
j\geq 1.\end{split}\right.$ (3.11)

Subsequently, the generic solution sequence $(e_{j}^{n})_{j\geq 1}$ can be
obtained as

$e_{j}^{n}=f_{j}^{n}+g_{j}^{n},\quad j\geq 1,\quad n\geq 1.$

####  3.2.1 The initial value problem (3.11)

It is first natural to investigate the initial value problem (3.11) since it
is really close to the infinite network case of the previous section. Here, we
consider the effect of the initial value assigned to each layer $j>0$ at the
first time step ($n=0$), except the input layer ($j=0$) which is set to zero.
The dynamics of (3.11) is still read out from the amplification factor
function $\rho$ defined in (3.6) and once again the solutions to (3.11) can be
obtained as the convolution of the initial sequence with the fundamental
solution associated to the problem. For $j_{0}\geq 1$, we denote by
$\boldsymbol{\delta}^{j_{0}}$ the Dirac delta sequence defined as
$\boldsymbol{\delta}_{j_{0}}^{j_{0}}=1$ and
$\boldsymbol{\delta}_{j}^{j_{0}}=0$ for all $j\geq 1$ and $j\neq j_{0}$.
Correspondingly, we denote by
$\mathcal{G}^{n}_{\mathrm{ivp}}(\cdot,j_{0})=(\mathcal{G}^{n}_{\mathrm{ivp}}(j,j_{0}))_{j\geq
1}$ the solution to (3.11) starting from $\boldsymbol{\delta}^{j_{0}}$, and
let us remark that the solutions to (3.11) starting from any initial condition
$(h_{j})_{j\geq 1}$ can be represented as

$f_{j}^{n}=\sum_{j_{0}=1}^{+\infty}\mathcal{G}^{n}_{\mathrm{ivp}}(j,j_{0})h_{j_{0}},\quad
j\geq 1,\quad n\geq 1.$

Combining the results of [12, 10] together with those of [13, 16, 17] which
precisely deal with recurrence equations with boundary conditions, one can
obtain very similar results as in the previous case. The very first obvious
remark that we can make is that for all $j,j_{0}\geq 1$ and $1\leq n<j_{0}$ we
have

$\mathcal{G}^{n}_{\mathrm{ivp}}(j,j_{0})=\mathcal{G}^{n}_{j-j_{0}},$

meaning that it takes $n=j_{0}$ iterations before the solution arrives at the
boundary $j=0$ and for $1\leq n<j_{0}$ the problem is similar to the one set
on the infinite network. This behavior is illustrated in Figure 7 for several
values of the hyper-parameters where we represent the spatio-temporal
evolution of the rescaled solution sequence
$(\sqrt{n}\,\mathcal{G}^{n}_{\mathrm{ivp}}(j,j_{0}))_{j\geq 1}$. We clearly
observe a Gaussian behavior before the solution reaches the boundary. And for
all $n\geq j_{0}$, we can write

$\mathcal{G}^{n}_{\mathrm{ivp}}(j,j_{0})=\mathcal{G}^{n}_{j-j_{0}}+\mathcal{G}_{\mathrm{bl}}^{n}(j,j_{0}),$

where $\mathcal{G}_{\mathrm{bl}}^{n}(j,j_{0})$ is a remainder term generated
by the boundary condition at $j=0$. It is actually possible to bound
$\mathcal{G}_{\mathrm{bl}}^{n}(j,j_{0})$ in each of the cases treated above.

Figure 7: Space-time evolution of the rescaled solution sequence
$(\sqrt{n}\,\mathcal{G}^{n}_{\mathrm{ivp}}(j,j_{0}))_{j\geq 1}$ to (3.11)
starting with a Dirac delta sequence at $j_{0}=25$ in different cases with
leftward propagation. (a) $\beta=0$ & $\alpha+\lambda<1$ with $c_{0}<0$. (b)
$\beta=0$ & $\alpha+\lambda=1$ with $c_{0}<0$. (c) $0<\beta<1$ &
$\alpha+\lambda=1$ with $-1<c_{\pi}<c_{0}<0$.

When $\beta=0$ and $\alpha+\lambda<1$ with $\alpha<\lambda$ such that
$c_{0}<0$, then $\mathcal{G}_{\mathrm{bl}}^{n}(j,j_{0})$ is well approximated
by

$\mathcal{G}_{\mathrm{bl}}^{n}(j,j_{0})\approx\left\\{\begin{array}[]{lc}-\frac{1}{\sqrt{4\pi\sigma_{0}n}}\exp\left(-\frac{|-j_{0}-c_{0}n|^{2}}{4\sigma_{0}n}\right)\left(\frac{\alpha}{\lambda}\right)^{j},&\quad
1\leq j\leq j_{0},\\\ e^{-\kappa n-\kappa(j-j_{0})},&\quad
j>j_{0},\end{array}\right.$

while when $\lambda<\alpha$ with $c_{0}>0$, then
$\mathcal{G}_{\mathrm{bl}}^{n}(j,j_{0})$ is well approximated by

$\mathcal{G}_{\mathrm{bl}}^{n}(j,j_{0})\approx\left\\{\begin{array}[]{lc}e^{-\kappa
n-\kappa(j_{0}-j)},&\quad 1\leq j\leq j_{0},\\\
-\frac{1}{\sqrt{4\pi\sigma_{0}n}}\exp\left(-\frac{|j-c_{0}n|^{2}}{4\sigma_{0}n}\right)\left(\frac{\lambda}{\alpha}\right)^{j_{0}},&\quad
j_{0}<j,\end{array}\right.$

this is illustrated in Figure 8 in the case $c_{0}<0$.

(a) Gaussian approximation. (b) Visualization of the boundary layer. (c)
Gaussian & boundary layer approximation.

Figure 8: Illustration of the solution
$\mathcal{G}^{n}_{\mathrm{ivp}}(j,j_{0})$ to (3.11) in the case where
$\beta=0$ and $\alpha+\lambda<1$ with $\alpha<\beta$. (a) Visualizations of
the solution $\mathcal{G}^{n}_{\mathrm{ivp}}(j,j_{0})$ (circles) at different
time iterations. The plain lines correspond to the Gaussian approximation
$\frac{1}{\sqrt{4\pi\sigma_{0}n}}\exp\left(-\frac{|j-j_{0}-c_{0}n|^{2}}{4\sigma_{0}n}\right)$
and remark the presence of a boundary layer (seen as a mismatch between the
circles and the Gaussian lines approximation). (b) We represent the boundary
layer by plotting (circles)
$\mathcal{G}^{n}_{\mathrm{ivp}}(j,j_{0})-\frac{1}{\sqrt{4\pi\sigma_{0}n}}\exp\left(-\frac{|j-j_{0}-c_{0}n|^{2}}{4\sigma_{0}n}\right)$
and we compare it to our boundary layer approximation
$-\frac{1}{\sqrt{4\pi\sigma_{0}n}}\exp\left(-\frac{|-j_{0}-c_{0}n|^{2}}{4\sigma_{0}n}\right)\left(\frac{\alpha}{\lambda}\right)^{j}$
(plain lines). (c) Finally we compare the solution
$\mathcal{G}^{n}_{\mathrm{ivp}}(j,j_{0})$ (circles) to its first order
approximation
$\frac{1}{\sqrt{4\pi\sigma_{0}n}}\exp\left(-\frac{|j-j_{0}-c_{0}n|^{2}}{4\sigma_{0}n}\right)-\frac{1}{\sqrt{4\pi\sigma_{0}n}}\exp\left(-\frac{|-j_{0}-c_{0}n|^{2}}{4\sigma_{0}n}\right)\left(\frac{\alpha}{\lambda}\right)^{j}$
(plain lines).

On the other hand for $\alpha+\lambda=1$ with $\alpha<\lambda$ such that
$c_{0}<0$, then $\mathcal{G}_{\mathrm{bl}}^{n}(j,j_{0})$ is well approximated
by

$\mathcal{G}_{\mathrm{bl}}^{n}(j,j,_{0})\approx\left\\{\begin{array}[]{lc}-\frac{1+(-1)^{n}}{\sqrt{4\pi\sigma_{0}n}}\exp\left(-\frac{|-j_{0}-c_{0}n|^{2}}{4\sigma_{0}n}\right)\left(\frac{\alpha}{\lambda}\right)^{j},&\quad
1\leq j\leq j_{0},\\\ e^{-\kappa n-\kappa(j-j_{0})},&\quad
j>j_{0},\end{array}\right.$

while when $\lambda<\alpha$ with $c_{0}>0$, then
$\mathcal{G}_{\mathrm{bl}}^{n}(j,j_{0})$ is well approximated by

$\mathcal{G}_{\mathrm{bl}}^{n}(j,j_{0})\approx\left\\{\begin{array}[]{lc}e^{-\kappa
n-\kappa(j_{0}-j)},&\quad 1\leq j\leq j_{0},\\\
-\frac{1+(-1)^{n}}{\sqrt{4\pi\sigma_{0}n}}\exp\left(-\frac{|j-c_{0}n|^{2}}{4\sigma_{0}n}\right)\left(\frac{\lambda}{\alpha}\right)^{j_{0}},&\quad
j_{0}<j.\end{array}\right.$

When $0<\beta<1$ and $\alpha+\lambda<1$ the approximations are similar as for
the case with $\beta=0$. We thus need to discuss three cases.

  * •

Case $-1<c_{\pi}<c_{0}<0$. In that case, we have for $1\leq j\leq j_{0}$ that

$\mathcal{G}_{\mathrm{bl}}^{n}(j,j_{0})\approx-\frac{1}{\sqrt{4\pi\sigma_{0}n}}\exp\left(-\frac{|-j_{0}-c_{0}n|^{2}}{4\sigma_{0}n}\right)\left(\frac{\alpha+\beta}{\lambda}\right)^{j}-\frac{(-1)^{n}}{\sqrt{4\pi\sigma_{\pi}n}}\exp\left(-\frac{|-j_{0}-c_{\pi}n|^{2}}{4\sigma_{\pi}n}\right)\left(\frac{\alpha-\beta}{\lambda}\right)^{j},$

with an exponential bound for $j>j_{0}$. This situation is presented in Figure
7(c)

  * •

Case $-1<c_{\pi}<0<c_{0}$. In this case we have

$\mathcal{G}_{\mathrm{bl}}^{n}(j,j_{0})\approx\left\\{\begin{array}[]{lc}-\frac{(-1)^{n}}{\sqrt{4\pi\sigma_{\pi}n}}\exp\left(-\frac{|-j_{0}-c_{\pi}n|^{2}}{4\sigma_{\pi}n}\right)\left(\frac{\alpha-\beta}{\lambda}\right)^{j},&\quad
1\leq j\leq j_{0},\\\
-\frac{1}{\sqrt{4\pi\sigma_{0}n}}\exp\left(-\frac{|j-c_{0}n|^{2}}{4\sigma_{0}n}\right)\left(\frac{\lambda}{\alpha+\beta}\right)^{j_{0}},&\quad
j_{0}<j<Ln.\end{array}\right.$

  * •

Case $-1<0<c_{\pi}<c_{0}$. In this case we have

$\mathcal{G}_{\mathrm{bl}}^{n}(j,j_{0})\approx-\frac{1}{\sqrt{4\pi\sigma_{0}n}}\exp\left(-\frac{|j-c_{0}n|^{2}}{4\sigma_{0}n}\right)\left(\frac{\lambda}{\alpha+\beta}\right)^{j_{0}}-\frac{(-1)^{n}}{\sqrt{4\pi\sigma_{\pi}n}}\exp\left(-\frac{|j-c_{\pi}n|^{2}}{4\sigma_{\pi}n}\right)\left(\frac{\lambda}{\alpha-\beta}\right)^{j_{0}}$

for $j_{0}<j<Ln$.

####  3.2.2 The boundary value problem (3.10)

We now turn our attention to the boundary value problem (3.10) where the
network is initialized with zero activity, for all layers except the input.
Motivated by applications, we will only focus on the case where
$s_{0}^{n}=s_{0}\in\mathbb{R}$ for all $n\geq 0$ (i.e., a constant sensory
input) and thus study:

$\left\\{\begin{split}g_{j}^{n+1}-\beta g_{j-1}^{n+1}&=\alpha
g_{j-1}^{n}+(1-\beta-\lambda-\alpha)g_{j}^{n}+\lambda g_{j+1}^{n},\quad j\geq
1,\quad n\geq 0,\\\ g_{0}^{n}&=s_{0},\quad n\geq 0,\\\ g_{j}^{0}&=0,\quad
j\geq 1.\end{split}\right.$ (3.12)

##### Case $\beta=0$.

Here, the stimulus information $s_{o}$ does not directly propagate through the
network via its feedforward connections (since $\beta=0$), but may still
propagate towards higher layers $j>0$ via the feedforward prediction error
correction mechanism, governed by parameter $\alpha$. When $\alpha+\lambda\leq
1$, we distinguish between three cases. Here and throughout, we denote by
$\mathrm{erf}$ the error function defined by

$\mathrm{erf}(x):=\frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{-z^{2}}\mathrm{d}z,\quad
x\in\mathbb{R}.$

  * •

Case $\alpha<\lambda$. In this case we have

$g_{j}^{n}=s_{0}\left(\frac{\alpha}{\lambda}\right)^{j}\left(1+\omega_{j}^{n}\right),\quad\text{
with }\quad\left|\omega_{j}^{n}\right|\leq Ce^{-\kappa n-\kappa j},\quad j\geq
1,\quad n\geq 1.$

It is interesting to note that the sequence
$\left(s_{0}\left(\frac{\alpha}{\lambda}\right)^{j}\right)_{j\leq 1}$ is a
stationary solution to (3.12) and we have uniform convergence at exponential
rate toward this stationary solution, that is

$\underset{j\geq
1}{\sup}\left|g_{j}^{n}-s_{0}\left(\frac{\alpha}{\lambda}\right)^{j}\right|\leq
Ce^{-\kappa n}\underset{n\rightarrow+\infty}{\longrightarrow}0.$

We illustrate this uniform convergence in Figure 9 (a)-(d).

  * •

Case $\alpha=\lambda$. We have

$\left|g_{j}^{n}-s_{0}\left(1-\mathrm{erf}\left(\frac{j}{\sqrt{4\sigma_{0}n}}\right)\right)\right|\leq\frac{C}{n}\exp\left(-\kappa\frac{j^{2}}{n}\right),\quad
j\geq 1,\quad n\geq 1.$

In this case, we observe a slow convergence to the steady state $s_{0}$.
Indeed, for each $\delta\in(0,1/2)$ we have

$\underset{n\rightarrow+\infty}{\lim}\,\underset{1\leq j\leq
n^{\delta}}{\sup}\,\left|g_{j}^{n}-s_{0}\right|=0,$

while for any $\delta>1/2$ we get

$\underset{n\rightarrow+\infty}{\lim}\,\underset{j\geq
n^{\delta}}{\sup}\,\left|g_{j}^{n}\right|=0.$

The propagation is thus diffusive along $j\sim\sqrt{n}$. This can be seen in
Figure 9 (b)-(e).

  * •

Case $\lambda<\alpha$. In this case we have

$\left|g_{j}^{n}-\frac{s_{0}}{2}\left(1-\mathrm{erf}\left(\frac{j-c_{0}n}{\sqrt{4\sigma_{0}n}}\right)\right)\right|\leq\frac{C}{\sqrt{n}}\exp\left(-\kappa\frac{(j-c_{0}n)^{2}}{n}\right),\quad
j\geq 1,\quad n\geq 1.$

In this case, we deduce that we have local uniform convergence towards the
steady state $s_{0}$, actually we have spreading at speed $c_{0}$. More
precisely, for any $c\in(0,c_{0})$ we have

$\underset{n\rightarrow+\infty}{\lim}\,\underset{1\leq j\leq
cn}{\sup}\,\left|g_{j}^{n}-s_{0}\right|=0,$

while for any $c>c_{0}$, we get

$\underset{n\rightarrow+\infty}{\lim}\,\underset{j\geq
cn}{\sup}\,\left|g_{j}^{n}\right|=0.$

We refer to Figure 9 (c)-(f) for an illustration. The figure clearly shows the
competition between hyperparameters $\alpha$ and $\lambda$, with forward
propagation of the sensory input only when $\alpha\geq\lambda$.

(a) $\alpha<\lambda$. (b) $\alpha=\lambda$. (c) $\lambda<\alpha$. (d)
$\alpha<\lambda$. (e) $\alpha=\lambda$. (f) $\lambda<\alpha$.

Figure 9: Case $\beta=0$ and $\alpha+\lambda\leq 1$. Visualization of the
solution $(g_{j}^{n})_{j\geq 1}$ of (3.12) when $s_{0}^{n}=1$ for all $n\geq
0$. Top row: space-time plots of the solution depending on $\alpha$ and
$\lambda$. Bottom row: solution profiles at different time iterations
(circles) compared with the leading order approximations. In the case
$\alpha<\lambda$, we observe uniform convergence (in space) at exponential
rate (in time) toward the stationary solution
$s_{0}\left(\left(\frac{\alpha}{\lambda}\right)^{j}\right)_{j\geq 1}$ while
when $\lambda<\alpha$ we observe the propagation of a wavefront where the
uniform steady state $s_{0}$ is propagating at speed $c_{0}>0$. In the
intermediate case where $\alpha=\lambda$ we get a diffusive invasion as
illustrated by the curve $j=\sqrt{4\sigma_{0}n}$.

##### Case $0<\beta<1$.

Here, the stimulus information $s_{o}$ propagates through the network not only
via its feedforward connections (governed by $\beta>0$) but also via the
feedforward prediction error correction mechanism, governed by parameter
$\alpha$. In the case where $\alpha+\lambda\leq 1$, the results from the case
$\beta=0$ remain valid, the only differences coming from the fact that the
above approximations in the case $\lambda\leq\alpha$ are only valid for $1\leq
j\leq Ln$ for some large constant $L>0$ with exponential localized bounds for
$j\geq Ln$ and that the steady state is now
$\left(s_{0}\left(\frac{\alpha+\beta}{\lambda}\right)^{j}\right)_{j\geq 1}$
whenever $\alpha+\beta<\lambda$. This confirms that the feedforward
propagation of the input $s_{0}$ is now dependent on both terms $\alpha$ and
$\beta$, jointly competing against the feedback term $\lambda$.

Let us remark that when $0<\beta<\alpha-\lambda$ and in the special case
$\alpha+\lambda=1$, where a second stable point exists for the amplification
factor function at $\rho(\pi)$, we can get a slightly more accurate
description of the solution in the form

$g_{j}^{n}=\frac{s_{0}}{2}\left(1-\mathrm{erf}\left(\frac{j-c_{0}n}{\sqrt{4\sigma_{0}n}}\right)\right)-\frac{s_{0}}{2(1+\beta)}\frac{(-1)^{j}}{\sqrt{4\pi\sigma_{\pi}n}}\exp\left(-\frac{(j-c_{\pi}n)^{2}}{4\sigma_{\pi}n}\right)+r_{j}^{n},\quad
n\geq 1,\quad 1\leq j\leq Ln,$

where the remainder term satisfies an estimate of the form

$\left|r_{j}^{n}\right|\leq\frac{C}{\sqrt{n}}\exp\left(-\kappa\frac{(j-c_{0}n)^{2}}{n}\right)+\frac{C}{n}\exp\left(-\kappa\frac{(j-c_{\pi}n)^{2}}{n}\right).$

This is illustrated in Figure 10. It should be noted here that, while the main
wavefront reflecting solution $c_{0}$ is a generic property of our network in
the entire range of validity of parameters $0\leq\beta<1$ and
$\alpha+\lambda\leq 1$, the second oscillatory pattern reflecting $c_{\pi}$
only appears in the special case of $\beta\neq 0$ and $\alpha+\lambda=1$. This
oscillation is, in fact, an artifact from the discrete formulation of our
problem, as will become evident in the next section, where we investigate
continuous formulations of the problem.

(a) Space-time plot. (b) Profile of the solution at $n=70$. (c) Zoom of the
profile near the boundary.

Figure 10: Case $0<\beta<1$ and $\alpha+\lambda=1$ with
$0<\beta<\alpha-\lambda$. Visualization of the solution $(g_{j}^{n})_{j\geq
1}$ of (3.12) when $s_{0}^{n}=1$ for all $n\geq 0$. The solution is the super-
position of a leading rightward front spreading at speed $c_{0}>0$ and an
oscillatory Gaussian profile propagating at speed $c_{\pi}>0$ with
$c_{\pi}<c_{0}$.

###  3.3 Towards continuous predictive models

Starting from a discrete approximation of our system made sense, not only for
mathematical convenience but also because artificial neural networks and deep
learning systems implementing similar predictive coding principles are
intrinsically discrete. Nonetheless, it can be useful to discard this discrete
approximation and investigate our system in the continuous limit. Note that in
the following, we will explore continuous extensions of our model in both time
and space. Biological neural networks, like any physical system, operate in
continuous time and thus it is more biologically accurate to relax the
temporal discretization assumption. This is what we do in the first part of
this section. In the spatial domain, however, the discretization of our system
into successive processing layers was not just an approximation, but also a
reflection of the hierarchical anatomy of the brain. Nonetheless, we can still
represent neuronal network depth continuously, even if only as a mathematical
abstraction. This is what we will do in the subsequent part of this section.
Understanding such continuous limits can allow us to test the robustness of
our framework, and to relate it to canonical models whose dynamics have been
more exhaustively characterized.

####  3.3.1 Continuous in time interpretation

As a first step, we present a continuous in time interpretation of the model
(3.9). We let $\Delta t>0$ be some parameter which will represent a time step
and reformulate the recurrence equation as

$(1-\beta)\frac{e_{j}^{n+1}-e_{j}^{n}}{\Delta t}=\frac{\beta}{\Delta
t}\left(e_{j-1}^{n+1}-e_{j}^{n+1}\right)+\frac{\lambda}{\Delta
t}\left(e_{j+1}^{n}-e_{j}^{n}\right)-\frac{\alpha}{\Delta
t}\left(e_{j}^{n}-e_{j-1}^{n}\right).$

We now interpret $e_{j}^{n}$ as the approximation of some smooth function of
time ${\bf e}_{j}(t)$ evaluated at $t_{n}:=n\Delta t$, that is
$e_{j}^{n}\sim{\bf e}_{j}(t_{n})$. As a consequence, we get that

$e_{j}^{n+1}\sim{\bf e}_{j}(t_{n+1})={\bf e}_{j}(t_{n}+\Delta t)={\bf
e}_{j}(t_{n})+\Delta t\frac{\mathrm{d}}{\mathrm{d}t}{\bf
e}_{j}(t_{n})+\mathcal{O}(|\Delta t|^{2}),\text{ as }\Delta t\rightarrow 0,$

such that

$\frac{e_{j}^{n+1}-e_{j}^{n}}{\Delta t}\sim\frac{\mathrm{d}}{\mathrm{d}t}{\bf
e}_{j}(t_{n}),\text{ as }\Delta t\rightarrow 0.$

Now, introducing the scaled parameters

$\widetilde{\beta}:=\frac{\beta}{\Delta
t},\quad\widetilde{\lambda}:=\frac{\lambda}{\Delta t},\text{ and
}\widetilde{\alpha}:=\frac{\alpha}{\Delta t},$

we get at the limit $\Delta t\rightarrow 0$ the following lattice ordinary
differential equation

$\frac{\mathrm{d}}{\mathrm{d}t}{\bf
e}_{j}(t)=(\widetilde{\beta}+\widetilde{\alpha}){\bf
e}_{j-1}(t)-(\widetilde{\beta}+\widetilde{\alpha}+\widetilde{\lambda}){\bf
e}_{j}(t)+\widetilde{\lambda}{\bf e}_{j+1}(t),\quad t>0.$ (3.13)

When defined on the infinite lattice $\mathbb{Z}$, one can represent the
solutions as

${\bf e}_{j}(t)=\sum_{k\in\mathbb{Z}}\mathbf{G}_{j-k}(t)h_{k},\quad
j\in\mathbb{Z},\quad t>0,$

starting from the initial sequence ${\bf e}_{j}(t=0)=(h_{j})_{j\in\mathbb{Z}}$
where $(\mathbf{G}_{j}(t))_{j\in\mathbb{Z}}$ is the fundamental solution to
(3.13) starting from the Dirac delta sequence $\boldsymbol{\delta}$. Once
again, each $\mathbf{G}_{j}(t)$ can be represented by the inverse Fourier
transform and reads

$\mathbf{G}_{j}(t)=\frac{1}{2\pi}\int_{-\pi}^{\pi}e^{\nu(\theta)t}e^{\mathbf{i}j\theta}\mathrm{d}\theta,\quad
t>0,\quad j\in\mathbb{Z},$

where the function $\nu(\theta)$ is defined as

$\nu(\theta):=(\widetilde{\beta}+\widetilde{\alpha})e^{-\mathbf{i}\theta}-(\widetilde{\beta}+\widetilde{\alpha}+\widetilde{\lambda})+\widetilde{\lambda}e^{\mathbf{i}\theta},\quad\theta\in[-\pi,\pi].$

The function $\nu(\theta)$ serves as an amplification factor function for the
time continuous equation (3.13). To ensure stability222Note that the notions
of stability/unstability and marginal stability introduced in the fully
discrete setting naturally extend to the semi-continuous setting., one needs
to impose that $\mathrm{Re}(\nu(\theta))\leq 0$ for each
$\theta\in[-\pi,\pi]$. From its formula, we obtain that

$\mathrm{Re}(\nu(\theta))=(\widetilde{\beta}+\widetilde{\alpha}+\widetilde{\lambda})(\cos(\theta)-1),\quad\theta\in[-\pi,\pi],$

such that we deduce that $\mathrm{Re}(\nu(0))=0$ and
$\mathrm{Re}(\nu(\theta))<0$ for all $\theta\in[-\pi,\pi]\backslash\\{0\\}$.
In particular, it is now evident that, contrary to the discrete case,
$\nu(\pi)$ cannot be a stable solution for the continuous system (except in
the trivial case where all hyperparameters
$\widetilde{\alpha},\widetilde{\beta},\widetilde{\lambda}$ are zero). This
confirms that the previously observed oscillations associated with $\rho(\pi)$
in specific cases were merely an artifact of the temporal discretization.

We note that, near the tangency point at $\theta=0$, the function
$\nu(\theta)$ has the following asymptotic expansion

$\nu(\theta)=-(\widetilde{\beta}+\widetilde{\alpha}+\widetilde{\lambda})\mathbf{i}\theta-\frac{\widetilde{\beta}+\widetilde{\alpha}+\widetilde{\lambda}}{2}\theta^{2}+\mathcal{O}(|\theta|^{3}),\text{
as }\theta\rightarrow 0.$

It is also possible to prove a Gaussian approximation in that case, and
following for example [6], we have

$\mathbf{G}_{j}(t)=\frac{1}{\sqrt{4\pi\widetilde{\sigma_{0}}t}}\exp\left(-\frac{|j-\widetilde{c_{0}}t|^{2}}{4\widetilde{\sigma_{0}}t}\right)+\mathbf{R}_{j}(t),\quad
j\in\mathbb{Z},\quad t>0,$

with

$\left|\mathbf{R}_{j}(t)\right|\leq\frac{C}{\sqrt{t}}\exp\left(-\kappa\frac{|j-\widetilde{c_{0}}t|^{2}}{t}\right),\quad
j\in\mathbb{Z},\quad t>0,$

for some universal constants $C>0$ and $\kappa>0$. Here, $\widetilde{c_{0}}$
and $\widetilde{\sigma_{0}}$ are given by

$\widetilde{c_{0}}=\widetilde{\beta}+\widetilde{\alpha}-\widetilde{\lambda},\quad\text{
and
}\quad\widetilde{\sigma_{0}}=\frac{\widetilde{\beta}+\widetilde{\alpha}+\widetilde{\lambda}}{2}>0.$

We remark that both $\widetilde{c_{0}}$ and $\widetilde{\sigma_{0}}$ are
linked to $c_{0}$ and $\sigma_{0}$ (the propagation speed and spread of the
solution in the case of the discrete model) in the following sense

$\frac{c_{0}}{\Delta t}\rightarrow\widetilde{c_{0}},\quad\text{ and
}\quad\frac{\sigma_{0}}{\Delta t}\rightarrow\widetilde{\sigma_{0}},\quad\text{
as }\Delta t\rightarrow 0.$

We also note that the spatially homogeneous solutions of (3.13) are trivial in
the sense that if we assume that ${\bf e}_{j}(t)={\bf e}(t)$ for all
$j\in\mathbb{Z}$ then the equation satisfied by ${\bf e}(t)$ is simply

$\frac{\mathrm{d}}{\mathrm{d}t}{\bf e}(t)=0.$

Finally, we conclude by noticing that in this continuous in time regime, there
is no possible oscillations either in space or time, in the sense that the
fundamental solution always resembles a fixed Gaussian profile advected at
wave speed $\widetilde{c_{0}}$. The formula for $\widetilde{c_{0}}$ highlights
the intuitive functional relation between the propagation (or advection)
direction and the “competition” between the feedforward influences
$\widetilde{\alpha}+\widetilde{\beta}$ and the feedback influence
$\widetilde{\lambda}$.

####  3.3.2 Fully continuous interpretation: both in time and depth

In this section, we give a possible physical interpretation of the discrete
model (3.9) via continuous transport equations, in which both time and space
(i.e., neuronal network depth) are made continuous. Let us introduce $\Delta
t>0$, $\Delta x>0$ and set $\nu:=\frac{\Delta x}{\Delta t}$. As before, we can
view ${\Delta t}$ as a time step for our system; additionally, ${\Delta x}$
can be viewed as a spatial step in the (continuous) neuronal depth dimension,
and thus $\nu$ becomes akin to a neural propagation speed or a conduction
velocity. We then rewrite the recurrence equation as

$(1-\beta)\frac{e_{j}^{n+1}-e_{j}^{n}}{\Delta
t}=\beta\nu\frac{e_{j-1}^{n+1}-e_{j}^{n+1}}{\Delta
x}+\lambda\nu\frac{e_{j+1}^{n}-e_{j}^{n}}{\Delta
x}-\alpha\nu\frac{e_{j}^{n}-e_{j-1}^{n}}{\Delta x}.$

The key idea is to now assume that $e_{j}^{n}$ represents an approximation of
some smooth function ${\bf e}(t,x)$ evaluated at $t_{n}:=n\Delta t$ and
$x_{j}:=j\Delta x$, that is $e_{j}^{n}\sim{\bf e}(t_{n},x_{j})$. Then passing
to the limit $\Delta t\to 0$, $\Delta x\to 0$ with $\frac{\Delta x}{\Delta
t}=\nu>0$ fixed and assuming that $\beta+\alpha\neq\lambda$, one gets the
partial differential equation

$\partial_{t}{\bf
e}(t,x)+\frac{\nu(\beta+\alpha-\lambda)}{1-\beta}\partial_{x}{\bf
e}(t,x)=0,\quad t>0,\quad x>0,$ (3.14)

with boundary condition ${\bf e}(t,x=0)=s_{0}(t)$ and initial condition ${\bf
e}(t=0,x)=h(x)$ satisfying the compatibility condition $s_{0}(0)=h(0)$ where
$s_{0}(t)$ is a smooth function such that $s_{0}(t^{n})=s_{0}^{n}$ and
$h(x_{j})=h_{j}$. The above partial differential equation is a transport
equation with associated speed $\frac{\nu(\beta+\alpha-\lambda)}{1-\beta}=\nu
c_{0}$. Depending on the sign of $c_{0}$, we have a different representation
for the solutions of (3.14).

  * •

Case $c_{0}<0$. Solution is given by

${\bf e}(t,x)=h\left(x-\nu c_{0}t\right),\quad t>0,\quad x>0.$

Let us remark that when $c_{0}<0$ the trace of the solution at $x=0$ is
entirely determined by the initial data $h(x)$ since

${\bf e}(t,x=0)=h\left(-\nu c_{0}t\right),\quad t>0.$

Intuitively, this reflects the dominance of backward (leftward) propagation in
this network, with solutions determined entirely by the initial value $h(x)$,
even for $x=0$ (the source term, $s_{0}(t)$, having no influence in this
case).

  * •

Case $c_{0}>0$. Solution is given by

${\bf e}(t,x)=\left\\{\begin{array}[]{lc}s_{0}\left(t-\frac{x}{\nu
c_{0}}\right),&x\leq\nu c_{0}t,\\\ h\left(x-\nu c_{0}t\right),&x>\nu
c_{0}t,\end{array}\right.\quad t>0,\quad x>0.$

Intuitively, this reflects the dominance of forward (rightward) propagation in
this network, with both the source term $s_{0}(t)$ and the initial values
$h(x)$ transported at constant velocity $\nu c_{0}$.

Thanks to the explicit form of the solutions, we readily obtain many
qualitative properties of the solution ${\bf e}(t,x)$. Boundedness and
positivity of the solutions are inherited from the functions $s_{0}(t)$ and
$h(x)$. In the case where $\beta+\alpha=\lambda$ (i.e., with balanced feed-
forward and feedback influences), the limiting equation is slightly different.
Indeed, in this case, introducing $\delta:=\frac{\Delta x^{2}}{\Delta t}$ and
letting $\Delta t\to 0$, $\Delta x\to 0$ with $\delta>0$ fixed, on gets the
partial differential equation

$\partial_{t}{\bf
e}(t,x)=\frac{\delta(\beta+\alpha+\lambda)}{2(1-\beta)}\partial_{x}^{2}{\bf
e}(t,x),\quad t>0,\quad x>0,$ (3.15)

and we readily observe that when $\beta+\alpha=\lambda$, we have that

$\frac{\beta+\alpha+\lambda}{2(1-\beta)}=\frac{\beta(1-\alpha-\lambda)+\alpha+\lambda-(\lambda-\alpha)^{2}}{2(1-\beta)^{2}}=\sigma_{0}>0.$

We obtain a heat equation with a boundary condition ${\bf e}(t,x=0)=s_{0}(t)$
and initial condition ${\bf e}(t=0,x)=h(x)$. Upon denoting

$\mathcal{S}(t,x,y):=\frac{1}{\sqrt{4\pi\delta\sigma_{0}t}}\left(e^{-\frac{(x-y)^{2}}{4\delta\sigma_{0}t}}-e^{-\frac{(x+y)^{2}}{4\delta\sigma_{0}t}}\right),$

the solution of the equation is given by

${\bf
e}(t,x)=s_{0}(t)+\int_{0}^{+\infty}\mathcal{S}(t,x,y)\left(h(y)-s_{0}(0)\right)\mathrm{d}y-\int_{0}^{t}\int_{0}^{+\infty}\mathcal{S}(t-s,x,y)s_{0}^{\prime}(s)\mathrm{d}y\mathrm{d}s,\quad
t>0,\quad x>0,$

Let us remark that when $s_{0}(t)=s_{0}\in\mathbb{R}$ is constant for all
$t\geq 0$, the above expression simplifies to

${\bf
e}(t,x)=s_{0}\left(1-\mathrm{erf}\left(\frac{x}{\sqrt{4\delta\sigma_{0}t}}\right)\right)+\int_{0}^{+\infty}\mathcal{S}(t,x,y)h(y)\mathrm{d}y,\quad
t>0,\quad x>0.$

In conclusion, this section extended our discrete model towards a continuous
limit in both space and time. In the temporal domain, it allowed us to
understand our stable solution as an advection behavior, and alerted us that
the other apparently oscillatory solutions previously observed in specific
cases were mainly due to our discretization approximation. In the spatio-
temporal domain, the continuous limit (3.14) allowed us to realize that our
main equation (3.1) was merely a discrete version of a transport equation.
In the following sections, we will systematically return to discrete
implementations (with gradually increasing functionality), before considering,
again, their continuous formulations.

##  4 Beyond the identity case

In the previous section we have studied in depth the case where
$\mathcal{W}^{f}$ and $\mathcal{W}^{b}$ are both the identity matrix: each
neuron in any given layer directly conveys its activation value to a single
corresponding neuron in the next layer, and to a single neuron in the previous
layer. Motivated by concrete implementations of the model in deep neural
networks [32, 9], we aim to investigate more realistic situations with more
complex connectivity matrices. While the generic unconstrained case (i.e. two
unrelated and dense connection matrices $\mathcal{W}^{f}$ and
$\mathcal{W}^{b}$) does not easily lend itself to analytical study, we will
consider here two situations of practical interest: in the first one, the
forward and backward connection matrices are symmetric and identical; in the
second case, each matrix is symmetric, but the two are not necessarily
identical.

###  4.1 The symmetric Rao & Ballard case

Following the pioneering work of Rao & Ballard [27], we will assume in this
section that $\mathcal{W}^{f}=(\mathcal{W}^{b})^{\bf t}$ and $\mathcal{W}^{f}$
is symmetric, which implies that

$\mathcal{W}^{f}=\mathcal{W}^{b}\in\mathscr{S}_{d}(\mathbb{R}),$

where we denoted $\mathscr{S}_{d}(\mathbb{R})$ the set of symmetric matrices
on $\mathbb{R}^{d}$.

The underlying interpretation is that, if a strong synaptic connection exists
from neuron $a$ to neuron $b$, then there is also a strong connection from $b$
to $a$. This assumption, which follows from Hebbian plasticity rules (“neurons
that fire together wire together”) does not capture all of the diversity of
brain connectivity patterns, but can be considered a good first approximation
[27].

####  4.1.1 Neural basis change and neural assemblies

The spectral theorem ensures that $\mathcal{W}^{f}$ (and thus
$\mathcal{W}^{b}$) is diagonalizable in an orthonormal basis. Namely, there
exists an orthogonal invertible matrix $P\in\mathscr{M}_{d}(\mathbb{R})$ such
that $PP^{\bf t}=P^{\bf t}P=\mathbf{I}_{d}$ and there exists a diagonal matrix
denoted $\mathcal{D}\in\mathscr{M}_{d}(\mathbb{R})$ such that

$\mathcal{W}^{f}=\mathcal{W}^{b}=P\mathcal{D}P^{\bf t}.$

We denote by $\gamma_{p}\in\mathbb{R}^{d}$, $p=1,\cdots,d$ the diagonal
elements of $\mathcal{D}$ and without loss of generality we may assume that

$\gamma_{1}\leq\cdots\leq\gamma_{d}.$

Thanks to this diagonalization, we can now perform a change of basis for our
neuronal space. We set $\mathcal{U}_{j}^{n}:=P^{t}\mathcal{E}_{j}^{n}$ as the
new basis, with $P\mathcal{U}_{j}^{n}:=\mathcal{E}_{j}^{n}$. Each
$\mathcal{U}_{j}^{n}$ can now be understood as a neural assembly, reflecting
one of the principal components of the weight matrix
$\mathcal{W}^{f}=\mathcal{W}^{b}$. Importantly, although assemblies may
overlap, activity updates induced by feedforward or feedback connections to
one given assembly do not affect the other assemblies, since the matrix $P$ is
orthogonal. Therefore, our problem is much simplified when considering
activity update equations at the level of these neural assemblies
$\mathcal{U}_{j}^{n}$ rather than across individual neurons
$\mathcal{E}_{j}^{n}$. Our model (2.1) becomes

$\mathcal{U}_{j}^{n+1}=\beta\mathcal{D}\mathcal{U}_{j-1}^{n+1}+\alpha\mathcal{D}\mathcal{U}_{j-1}^{n}+\left[(1-\beta-\lambda)\mathbf{I}_{d}-\alpha\mathcal{D}^{2}\right]\mathcal{U}_{j}^{n}+\lambda\mathcal{D}\mathcal{U}_{j+1}^{n}.$

Note that, because all matrices in the above equation are diagonal, we have
totally decoupled the $d$ components of the vector $\mathcal{U}_{j}^{n}$. More
precisely, by denoting $u_{j,p}^{n}$ the $p$th component of
$\mathcal{U}_{j}^{n}$, that is
$\mathcal{U}_{j}^{n}=(u_{j,1}^{n},\cdots,u_{j,d}^{n})^{\mathbf{t}}$, we obtain

$u_{j,p}^{n+1}-\beta\gamma_{p}u_{j-1,p}^{n+1}=\alpha\gamma_{p}u_{j-1,p}^{n}+(1-\beta-\lambda-\alpha\gamma_{p}^{2})u_{j,p}^{n}+\lambda\gamma_{p}u_{j+1,p}^{n},\quad
p=1,\cdots,d.$

This indicates that one needs to study

$u_{j}^{n+1}-\beta\gamma u_{j-1}^{n+1}=\alpha\gamma
u_{j-1}^{n}+(1-\beta-\lambda-\alpha\gamma^{2})u_{j}^{n}+\lambda\gamma
u_{j+1}^{n},$ (4.1)

where $\gamma\in\mathbb{R}$ is a given parameter. Here, $\gamma$ can be
thought of as the connection strength across layers (both feedforward and
feedback, since we assumed here symmetric connectivity) of the neural assembly
under consideration. By construction, each assembly in a given layer is only
connected to the corresponding assembly in the layer above, and similarly in
the layer below. Note that when $\gamma=1$, we encounter again the exact
situation that we studied in the previous section (3.9), but now with neural
assemblies in lieu of individual neurons.

####  4.1.2 Study of the amplification factor function

Based on our previous analysis, the behavior of the solutions to (4.1) are
intrinsically linked to the properties of the amplification factor function:

$\rho_{\gamma}(\theta)=\frac{\alpha\gamma\left(e^{-\mathbf{i}\theta}-\gamma\right)+1-\beta+\lambda\left(\gamma
e^{\mathbf{i}\theta}-1\right)}{1-\beta\gamma
e^{-\mathbf{i}\theta}},\quad\theta\in[-\pi,\pi],$

where one needs to ensure that $|\rho_{\gamma}(\theta)|\leq 1$ for all
$\theta\in[-\pi,\pi]$. The very first condition is to ensure that
$1\neq\beta|\gamma|$ (to avoid division by zero). Next, we investigate the
behavior of $\rho_{\gamma}(\theta)$ at $\theta=0$ and check under which
condition on $\gamma$ we can ensure that $-1\leq\rho_{\gamma}(0)\leq 1$. We
have

$\rho_{\gamma}(0)=1+\frac{(1-\gamma)(\alpha\gamma-\lambda-\beta)}{1-\beta\gamma},$

which readily tells us that $\rho_{\gamma}(0)=1$ if and only if $\gamma=1$ or
$\gamma=\frac{\lambda+\beta}{\alpha}$. And on the other hand
$\rho_{\gamma}(0)=-1$ if and only if $\gamma=\gamma_{\pm}^{0}$ with

$\gamma_{\pm}^{0}=\frac{\lambda+\alpha-\beta\pm\sqrt{(\lambda+\alpha-\beta)^{2}+4\alpha(2-\lambda-\beta)}}{2\alpha},$

with $\gamma_{-}^{0}<0<\gamma_{+}^{0}$ since $\lambda+\beta<2$ by assumption.
One also notices that

$(\lambda+\alpha-\beta)^{2}+4\alpha(2-\lambda-\beta)=(\lambda+3\alpha-\beta)^{2}+8\alpha(1-\alpha-\lambda),$

such that either $\alpha+\lambda=1$ and in that case $\gamma_{-}^{0}=-1$ and
$\gamma_{+}^{0}=1+\frac{1-\beta}{\alpha}>1$, or $\alpha+\lambda<1$ and in that
case $\gamma_{-}^{0}<-1$ and $\gamma_{+}^{0}>1+\frac{1-\beta}{\alpha}>1$.
Next, we remark that

$\rho_{\gamma}(\pi)=1-\frac{(1+\gamma)(\alpha\gamma+\lambda+\beta)}{1+\beta\gamma}=\rho_{-\gamma}(0),$

which then implies that $\rho_{\gamma}(\pi)=1$ if and only if $\gamma=-1$ or
$\gamma=-\frac{\lambda+\beta}{\alpha}$ and $\rho_{\gamma}(\pi)=-1$ if and only
if $\gamma=-\gamma_{\pm}^{0}$.

As explained in the beginning, our aim is to completely characterize under
which conditions on $\gamma\in\mathbb{R}$, $0\leq\beta<1$ and
$0<\alpha,\lambda<1$ with $\alpha+\lambda\leq 1$, one can ensure that
$|\rho_{\gamma}(\theta)|\leq 1$ for all $\theta\in[-\pi,\pi]$.

Potential regions of marginal stability are thus given by those values of the
parameters satisfying $\gamma=\pm 1$,
$\gamma=\pm\frac{\lambda+\beta}{\alpha}$, $\gamma=\pm\gamma_{+}^{0}$ and
$\gamma=\pm\gamma_{-}^{0}$, and it is important to determine the intersections
among the above regions. We have already proved that $\gamma_{-}^{0}=-1$
whenever $\alpha=1-\lambda$. Next, we compute that
$\gamma_{-}^{0}=-\frac{\lambda+\beta}{\alpha}$ whenever
$\alpha=\frac{\lambda(\lambda+\beta)}{1-\lambda-\beta}=:\Lambda$, while
$\gamma_{-}^{0}=-\gamma_{+}^{0}$ whenever $\alpha=\beta-\lambda$ and
$\gamma_{+}^{0}=\frac{\lambda+\beta}{\alpha}$ when
$\alpha=\beta(\lambda+\beta)$. Let us already point out that $\Lambda$ is only
defined if $\lambda+\beta<1$ and in that case $\Lambda>0$.

(a) $(\lambda,\beta)$ in Region $\mathbf{(I)}$. (b) $(\lambda,\beta)$ in
Region $\mathbf{(II)}$. (c) $(\lambda,\beta)$ in Region $\mathbf{(III)}$. (d)
$(\lambda,\beta)$ in Region $\mathbf{(IV)}$. (e) $(\lambda,\beta)$ in Region
$\mathbf{(V)}$.

Figure 11: Stability/instability regions and their boundaries as a function of
$(\alpha,\gamma)$ for (4.1) while $(\lambda,\beta)$ being fixed in one the
five regions given in panel (a). Shaded orange regions correspond to an
instability for (4.1) while purple regions correspond to a stability for
(4.1). The boundaries of the stability/instability regions are given by the
intersections of the parametrized curves $\gamma=\pm 1$ (dark red curves),
$\gamma=\pm\gamma_{-}^{0}$ (dark blue curves), $\gamma=\pm\gamma_{+}^{0}$
(pink curves) and $\gamma=\pm\frac{\lambda+\beta}{\alpha}$ (magenta curves).
Note that the region of interest is $0<\alpha\leq 1-\lambda$. Along each such
parametrized curves equation (4.1) is marginally stable.

We now introduce five regions in the quadrant
$(\beta,\lambda)\in[0,1)\times(0,1)$ which are depicted in Figure 11(a).
First, since $1-\lambda=\Lambda=\lambda+\beta$ if and only if
$2\lambda+\beta=1$ (which corresponds to the blue line in Figure 11(a)), we
deduce that when $2\lambda+\beta\geq 1$ we have
$1-\lambda\leq\min(\Lambda,\lambda+\beta)$ which leads us to define the
following two regions:

$\displaystyle\mathbf{(I)}$
$\displaystyle:=\left\\{(\beta,\lambda)\in[0,1)\times(0,1)~{}|~{}2\lambda+\beta\geq
1\text{ and }\beta\leq\lambda\right\\},$ $\displaystyle\mathbf{(II)}$
$\displaystyle:=\left\\{(\beta,\lambda)\in[0,1)\times(0,1)~{}|~{}2\lambda+\beta\geq
1\text{ and }\beta>\lambda\right\\}.$

Now, when $2\lambda+\beta<1$, we have the strict ordering
$0<\Lambda<\lambda+\beta<1-\lambda$ and when $\beta>\lambda$ it is thus
necessary to compare $\Lambda$ to $\beta-\lambda$. We remark that
$\Lambda=\beta-\lambda=\beta(\lambda+\beta)$ if and only if
$\beta^{2}-\beta(1-\lambda)+\lambda=0$, which corresponds to the yellow
parabola in Figure 11(a). We thus define the following three regions:

$\displaystyle\mathbf{(III)}$
$\displaystyle:=\left\\{(\beta,\lambda)\in[0,1)\times(0,1)~{}|~{}2\lambda+\beta<1\text{
and }\beta\leq\lambda\right\\},$ $\displaystyle\mathbf{(IV)}$
$\displaystyle:=\left\\{(\beta,\lambda)\in[0,1)\times(0,1)~{}|~{}2\lambda+\beta<1,\beta>\lambda\text{
and }\beta^{2}-\beta(1-\lambda)+\lambda\geq 0\right\\},$
$\displaystyle\mathbf{(V)}$
$\displaystyle:=\left\\{(\beta,\lambda)\in[0,1)\times(0,1)~{}|~{}2\lambda+\beta<1,\beta>\lambda\text{
and }\beta^{2}-\beta(1-\lambda)+\lambda<0\right\\}.$

Note that when $(\beta,\lambda)$ is in region $\mathbf{(IV)}$, we have the
ordering

$0<\beta-\alpha\leq\Lambda<\lambda+\beta<1-\lambda,$

while for $(\beta,\lambda)$ in region $\mathbf{(V)}$, we have

$0<\Lambda<\beta(\lambda+\beta)<\beta-\alpha<\lambda+\beta<1-\lambda.$

We can characterize the stability of our equation separately for each of the
five regions defined in Figure 11(a). Since the region already determines the
value of the parameters $\beta$ and $\lambda$, the stability will be expressed
as a function of the two remaining parameters $\alpha$ and $\gamma$ (Figure
11(b-f)). We refer to Figures 11(b-f) for a comprehensive representation of
the stability regions. Note that the boundaries of the stability/instability
regions are precisely given by the intersections of the parametrized curves
$\gamma=\pm 1$ (dark red curves), $\gamma=\pm\gamma_{-}^{0}$ (dark blue
curves), $\gamma=\pm\gamma_{+}^{0}$ (pink curves) and
$\gamma=\pm\frac{\lambda+\beta}{\alpha}$ (magenta curves). Along each such
parametrized curves equation (4.1) is marginally stable. We comment below the
case $(\lambda,\beta)$ in Region $\mathbf{(III)}$. The other cases can be
described in the same way, but we leave this out for conciseness.

| $\gamma=1$ | $\gamma=-1$ | $\gamma=\frac{\lambda+\beta}{\alpha}$ | $\gamma=-\frac{\lambda+\beta}{\alpha}$ | $\gamma=\gamma_{\pm}^{0}$ | $\gamma=-\gamma_{\pm}^{0}$
---|---|---|---|---|---|---
$c_{0}^{\gamma}$ | $\frac{\beta+\alpha-\lambda}{1-\beta}$ | $\frac{\alpha-\lambda-\beta}{1+\beta}$ | $\frac{(\beta+\lambda)(\beta+\alpha-\lambda)}{\alpha-\beta(\lambda+\beta)}$ |  | $\frac{\gamma_{\pm}^{0}(\beta+\lambda-\alpha)}{1-\beta\gamma_{\pm}^{0}}$ |
$c_{\pi}^{\gamma}$ | $\frac{\alpha-\lambda-\beta}{1+\beta}$ | $\frac{\beta+\alpha-\lambda}{1-\beta}$ |  | $\frac{(\beta+\lambda)(\beta+\alpha-\lambda)}{\alpha-\beta(\lambda+\beta)}$ |  | $\frac{\gamma_{\pm}^{0}(\beta+\lambda-\alpha)}{1-\beta\gamma_{\pm}^{0}}$
Table 1: Expressions of the wave speed $c_{0}^{\gamma}$ and $c_{\pi}^{\gamma}$
for values of $\gamma$ corresponding the boundaries of the stability regions
from Figure 11. Let us note that $c_{\pi}^{1}$ and $c_{0}^{-1}$ only exist
when $\alpha+\lambda=1$.

Suppose that $(\lambda,\beta)$ belongs to Region $\mathbf{(III)}$. We present
the results of Figure 11(d) by letting $\alpha$ vary between $0$ and
$1-\lambda$ and $\gamma\in\mathbb{R}$. More precisely, for each fixed
$\alpha\in(0,1-\lambda)$ we investigate the stability properties as a function
of $\gamma$. We have to distinguish between several subcases.

  1. (i)

If $0<\alpha<\Lambda$. Then, equation (4.1) is stable for each
$\gamma\in(-1,1)$, unstable for $|\gamma|>1$ and marginally stable at
$\gamma=\pm 1$ with $\rho_{1}(0)=1$ and $\rho_{-1}(\pm\pi)=1$.

  2. (ii)

If $\alpha=\Lambda$. Then $\gamma_{-}^{0}=-\frac{\lambda+\beta}{\alpha}$ and
equation (4.1) is stable for each $\gamma\in(-1,1)$, unstable for
$|\gamma|>|\gamma_{-}^{0}|$ and $|\gamma_{-}^{0}|>|\gamma|>1$, whereas it is
marginally stable at $\gamma=\pm 1$ and at $\gamma=\pm\gamma_{-}^{0}$ with
$\rho_{1}(0)=1$, $\rho_{-1}(\pm\pi)=1$, $\rho_{\gamma_{-}^{0}}(0)=-1$ and
$\rho_{-\gamma_{-}^{0}}(\pm\pi)=-1$ together with
$\rho_{-\gamma_{-}^{0}}(0)=1$, $\rho_{\gamma_{-}^{0}}(\pm\pi)=1$.

  3. (iii)

If $\Lambda<\alpha<\lambda+\beta$. Then, equation (4.1) is stable for each
$\gamma\in(-1,1)$ and
$\frac{\lambda+\beta}{\alpha}<|\gamma|<|\gamma_{-}^{0}|$, unstable for
$|\gamma|>|\gamma_{-}^{0}|$ and $\frac{\lambda+\beta}{\alpha}>|\gamma|>1$ and
marginally stable at $\gamma\in\left\\{\pm
1,\pm\frac{\lambda+\beta}{\alpha},\pm\gamma_{-}^{0}\right\\}$ with
$\rho_{1}(0)=1$, $\rho_{-1}(\pm\pi)=1$,
$\rho_{\frac{\lambda+\beta}{\alpha}}(0)=1$,
$\rho_{-\frac{\lambda+\beta}{\alpha}}(\pm\pi)=1$,
$\rho_{\gamma_{-}^{0}}(0)=-1$ and $\rho_{-\gamma_{-}^{0}}(\pm\pi)=-1$.

  4. (iv)

If $\alpha=\lambda+\beta$. Then, equation (4.1) is stable for each
$\gamma\in(\gamma_{-}^{0},-\gamma_{-}^{0})\backslash\\{\pm 1\\}$, unstable for
$|\gamma|>|\gamma_{-}^{0}|$ and marginally stable at $\gamma\in\left\\{\pm
1,\pm\gamma_{-}^{0}\right\\}$ with $\rho_{1}(0)=1$, $\rho_{-1}(\pm\pi)=1$,
$\rho_{\gamma_{-}^{0}}(0)=-1$ and $\rho_{-\gamma_{-}^{0}}(\pm\pi)=-1$. Remark
that in this case we have $\frac{\lambda+\beta}{\alpha}=1$.

  5. (v)

If $\lambda+\beta<\alpha<1-\lambda$. Then, equation (4.1) is stable for each
$\gamma\in\left(-\frac{\lambda+\beta}{\alpha},\frac{\lambda+\beta}{\alpha}\right)$
and $1<|\gamma|<|\gamma_{-}^{0}|$, unstable for $|\gamma|>|\gamma_{-}^{0}|$
and $1>|\gamma|>\frac{\lambda+\beta}{\alpha}$ and marginally stable at
$\gamma\in\left\\{\pm
1,\pm\frac{\lambda+\beta}{\alpha},\pm\gamma_{-}^{0}\right\\}$ with
$\rho_{1}(0)=1$, $\rho_{-1}(\pm\pi)=1$,
$\rho_{\frac{\lambda+\beta}{\alpha}}(0)=1$,
$\rho_{-\frac{\lambda+\beta}{\alpha}}(\pm\pi)=1$,
$\rho_{\gamma_{-}^{0}}(0)=-1$ and $\rho_{-\gamma_{-}^{0}}(\pm\pi)=-1$.

  6. (vi)

If $\alpha=1-\lambda$. Then, equation (4.1) is stable for each
$\gamma\in\left(-\frac{\lambda+\beta}{\alpha},\frac{\lambda+\beta}{\alpha}\right)$,
unstable for $|\gamma|>|\gamma_{-}^{0}|=1$ and
$1>|\gamma|>\frac{\lambda+\beta}{\alpha}$ and marginally stable at
$\gamma\in\left\\{\pm 1,\pm\frac{\lambda+\beta}{\alpha}\right\\}$ with
$\rho_{1}(0)=1$ and $\rho_{1}(\pm\pi)=-1$, $\rho_{-1}(\pm\pi)=1$ and
$\rho_{-1}(0)=-1$, with $\rho_{\frac{\lambda+\beta}{\alpha}}(0)=1$,
$\rho_{-\frac{\lambda+\beta}{\alpha}}(\pm\pi)=1$. Remark that in this case we
have $\gamma_{-}^{0}=-1$.

##### Summary.

In summary, we see that stability is nearly guaranteed whenever $-1<\gamma<1$,
regardless of the values of other parameters (as long as
$0<\alpha<1-\lambda$). This makes intuitive sense, as $\gamma$ represents the
connection strength across layers of a particular neural assembly, a
connection weight $|\gamma|<1$ implies that activity of this assembly will
remain bounded across layers. Additionally, and perhaps more interestingly, in
some but not all regions (e.g. Regions II and V) stability can be obtained for
much larger values of $|\gamma|$; this, however, appears to coincide with low
values of the $\alpha$ parameter. In other words, for high connection
strengths $|\gamma|$, the feedforward error correction term $\alpha$ makes the
system unstable.

####  4.1.3 Wave speed characterization

In the previous section (The Identity Case), we have proved that the direction
of propagation was given by the sign of $c_{0}$ and $c_{\pi}$ whenever they
exist which could be read off from the behavior of $\rho_{\gamma}(\theta)$
near $\theta=0$ or $\theta=\pm\pi$. We have reported the values of
$c_{0}^{\gamma}$ and $c_{\pi}^{\gamma}$ for different values of $\gamma$ in
Table 1. For example, in Figure 12 we illustrate the changes in propagation
speed and direction for $c_{0}^{\gamma}$ for the case $(\lambda,\beta)$ in
Region $\mathbf{(III)}$ (as defined in Figure 11a), but the calculations
remain valid for the other regions.

Figure 12: Plot of the wave speeds $c_{0}^{1}$,
$c_{0}^{\frac{\lambda+\beta}{\alpha}}$ and $c_{0}^{\gamma_{-}^{0}}$ as a
function of $\alpha$ in the case $(\lambda,\beta)$ in Region $\mathbf{III}$
associated to Figure 11(d).

It is worth emphasizing that for fixed values of the hyper-parameters
$\alpha$, $\beta$ and $\lambda$, we see here that varying $\gamma$ can give
rise to different propagation speeds or even different directions. As each
neuronal assembly $u_{j,p}$ in a given layer $j$ is associated with its own
connection strength $\gamma_{p}$, it follows that different speeds and even
different directions of propagation can concurrently be obtained in a single
network, one for each assembly. For instance, in a given network with
hyperparameters $\alpha=0.2$, $\beta=0.2$ and $\lambda=0.3$ (region
$\mathbf{III}$), a neural assembly with a connection strength of $\gamma=1$
would propagate forward at a relatively slow speed, while another with
$\gamma=2.5$ would propagate in the same direction at a much faster speed, and
yet another assembly with $\gamma=\gamma^{0}_{-}\approx-2.09$ would
simultaneously propagate in the opposite backward direction.

####  4.1.4 Continuous in time interpretation

We can repeat the “continuous system” analysis conducted in the previous
section (The Identity Case), which has lead to (3.13), but this time with Rao-
Ballard connection matrices between layers. With the same scaling on the
hyperparameters

$\widetilde{\beta}:=\frac{\beta}{\Delta
t},\quad\widetilde{\lambda}:=\frac{\lambda}{\Delta t},\text{ and
}\widetilde{\alpha}:=\frac{\alpha}{\Delta t},$

we get that, at the limit $\Delta t\rightarrow 0$, the equation (4.1) becomes
the following lattice ordinary differential equation

$\frac{\mathrm{d}}{\mathrm{d}t}{\bf
u}_{j}(t)=(\widetilde{\beta}+\widetilde{\alpha})\gamma{\bf
u}_{j-1}(t)-(\widetilde{\beta}+\widetilde{\lambda}+\widetilde{\alpha}\gamma^{2}){\bf
u}_{j}(t)+\widetilde{\lambda}\gamma{\bf u}_{j+1}(t),\quad t>0.$ (4.2)

Note that the neuronal layer activity is now expressed in terms of neural
assemblies ${\bf u}_{j}$ rather than individual neurons ${\bf e}_{j}$.

The amplification factor function in this case is given by

$\nu_{\gamma}(\theta)=(\widetilde{\beta}+\widetilde{\alpha})\gamma
e^{-\mathbf{i}\theta}-(\widetilde{\beta}+\widetilde{\lambda}+\widetilde{\alpha}\gamma^{2})+\widetilde{\lambda}\gamma
e^{\mathbf{i}\theta},\quad\theta\in[-\pi,\pi],$

whose real part is given by

$\mathrm{Re}(\nu_{\gamma}(\theta))=(\widetilde{\beta}+\widetilde{\alpha}+\widetilde{\lambda})\gamma\cos(\theta)-(\widetilde{\beta}+\widetilde{\lambda}+\widetilde{\alpha}\gamma^{2}),\quad\theta\in[-\pi,\pi].$

When $\gamma>0$, we observe that

$\underset{\theta\in[-\pi,\pi]}{\max}\mathrm{Re}(\nu_{\gamma}(\theta))=\mathrm{Re}(\nu_{\gamma}(0))=(\widetilde{\lambda}+\widetilde{\beta}-\widetilde{\alpha}\gamma)(\gamma-1),$

whereas when $\gamma<0$, we have

$\underset{\theta\in[-\pi,\pi]}{\max}\mathrm{Re}(\nu_{\gamma}(\theta))=\mathrm{Re}(\nu_{\gamma}(\pm\pi))=-(\widetilde{\lambda}+\widetilde{\beta}+\widetilde{\alpha}\gamma)(\gamma+1).$

As a consequence, the stability analysis in this case is very simple and
depends only on the relative position of $\gamma$ with respect to $\pm 1$ and
$\pm\frac{\widetilde{\lambda}+\widetilde{\beta}}{\widetilde{\alpha}}$. It is
summarized in Figure 13.

Figure 13: Stability/instability regions and their boundaries as a function of
$(\widetilde{\alpha},\gamma)$ for (4.2) for any
$(\widetilde{\lambda},\widetilde{\beta})$ fixed. The shaded orange region
corresponds to an instability for (4.2) while the purple region corresponds to
a stability for (4.2). The boundaries of the stability/instability regions are
given by the intersections of the parametrized curves $\gamma=\pm 1$ (dark red
curves) and $\gamma=\pm\frac{\lambda+\beta}{\alpha}$ (magenta curves) where
equation (4.2) is marginally stable.

The simple behavior illustrated in Figure 13 for our continuous system
contrasts with the number and diversity of behaviors obtained for the discrete
version of the same system (Figure 11). A number of points are worth
highlighting. For instance, although the values of $\beta$ and $\lambda$ were
critical for the discrete system (to define the region (I) to (V)), they do
not affect the qualitative behavior of the continuous system. Furthermore,
some observations in the continuous system appear to contradict the
conclusions made previously in the discrete case. We see that stability can
still be obtained with high values of the connection weight $\gamma>>1$, but
this time the stable regions coincide with high $\alpha$ values, whereas it
was the opposite in Figure 11 panels (b),(f). This qualitative difference in
behavior can be taken as a point of caution, to remind us that a discrete
approximation of the system can be associated with important errors in
interpretation.
Finally we note that, while stability regions are qualitatively different in
the continuous case compared to the discrete approximation, the speed and
direction of propagation of neural signals (reflected in the variables $c_{0}$
and $c_{\pi}$ when they exist) remains comparable.

####  4.1.5 A class of examples

In this section, we provide a class of examples of $\mathcal{W}^{f}$ amenable
to a complete analysis. Namely we consider $\mathcal{W}^{f}$ as the following
linear combination

$\mathcal{W}^{f}=\zeta\mathbf{I}_{d}+\xi\mathbf{A},$ (4.3)

for some $\zeta,\xi\in\mathbb{R}$ where
$\mathbf{A}\in\mathscr{M}_{d}(\mathbb{R})$ is given by

$\mathbf{A}=\left(\begin{matrix}-2&1&0&\cdots&\cdots&0\\\
1&-2&1&\ddots&\ddots&\vdots\\\ 0&\ddots&\ddots&\ddots&\ddots&\vdots\\\
\vdots&\ddots&\ddots&\ddots&\ddots&\vdots\\\ \vdots&\ddots&\ddots&1&-2&1\\\
0&\cdots&\cdots&0&1&-2\end{matrix}\right).$

The matrix $\mathbf{A}$ is nothing but the discrete laplacian and
$\mathcal{W}^{f}$ acts as a convolution operator on $\mathbb{R}^{d}$. More
precisely, $\mathcal{W}^{f}$ combines a convolution term with a residual
connection term, as in the well-known ResNet architecture [18]. Let us also
note that the spectrum of $\mathbf{A}$ is well known and given by

$\mathrm{Spec}(\mathbf{A})=\left\\{-4\sin^{2}\left(\frac{p\pi}{2(d+1)}\right),\quad
p=1,\cdots,d\right\\}.$

As a consequence, the spectrum of $\mathcal{W}^{f}$ is simply given by

$\mathrm{Spec}(\mathbf{\mathcal{W}^{f}})=\left\\{\zeta-4\xi\sin^{2}\left(\frac{p\pi}{2(d+1)}\right),\quad
p=1,\cdots,d\right\\}.$

One can for example set

$\zeta=\frac{\sin^{2}\left(\frac{d\pi}{2(d+1)}\right)+\sin^{2}\left(\frac{\pi}{2(d+1)}\right)}{\sin^{2}\left(\frac{d\pi}{2(d+1)}\right)-\sin^{2}\left(\frac{\pi}{2(d+1)}\right)}\quad\text{
and
}\quad\xi=\frac{1}{2\left(\sin^{2}\left(\frac{d\pi}{2(d+1)}\right)-\sin^{2}\left(\frac{\pi}{2(d+1)}\right)\right)},$

such that

$\zeta-4\xi\sin^{2}\left(\frac{d\pi}{2(d+1)}\right)=-1,\quad\zeta-4\xi\sin^{2}\left(\frac{\pi}{2(d+1)}\right)=1,$

and for all $p=2,\cdots,d-1$

$\zeta-4\xi\sin^{2}\left(\frac{p\pi}{2(d+1)}\right)\in(-1,1).$

Next, for any $p=1,\cdots,d$ the eigenvector corresponding to the eigenvalue
$-4\sin^{2}\left(\frac{p\pi}{2(d+1)}\right)$ is

$U_{p}=\left(\cos\left(\frac{p\pi}{d+1}\right),\cdots,\cos\left(\frac{pk\pi}{d+1}\right),\cdots,\cos\left(\frac{pd\pi}{d+1}\right)\right)^{\bf
t}\in\mathbb{R}^{d}.$

$U_{p}$ is the projection vector that corresponds to the $p^{th}$ neural
assembly $u_{j,p}$ as defined above.
Along $U_{1}$, the recurrence equation reduces to (4.1) with $\gamma=1$, while
along $U_{d}$, the recurrence equation reduces to (4.1) with $\gamma=-1$, and
we can apply the results of the previous section (the Identity case). In
between (for all $1\leq p\leq d$) we see that the eigenvalues of our
connection matrix $\mathbf{\mathcal{W}^{f}}$ span the entire range between
$-1$ and $1$, that they can be explicitly computed, and thus that the
stability, propagation speed and direction of activity in the corresponding
neural assembly can be determined.

####  4.1.6 Fully continuous interpretation in time, depth and width.

For the same class of example (connection matrix composed of a convolution and
residual terms), we now wish to provide a fully continuous interpretation for
model (2.1) in the special case $\zeta=1$ and $\xi$ adjusted as follows. By
fully continuous, we mean that we explore the limit of our model when not only
time $t$, but also network depth $x$ and neuronal layer width $y$ are
considered as continuous variables. Although we already presented a model that
was continuous in both time and depth in subsection 3.3.2, the layers in that
model only comprised a single neuron, and had no intrinsic spatial dimension.
We now introduce this third continuous dimension. The starting point is to see
$\mathcal{E}_{j,k}^{n}$, the $k$th element of $\mathcal{E}_{j}^{n}$, as an
approximations of some continuous function $\mathcal{E}(t,x,y)$ evaluated at
$t_{n}=n\Delta t$, $x_{j}=j\Delta x$ and $y_{k}=k\Delta y$ for some $\Delta
t>0$, $\Delta x>0$ and $\Delta y>0$. Let us first remark that the action of
$\mathbf{A}$ on $\mathcal{E}_{j}^{n}$ is given by

$(\mathbf{A}\mathcal{E}_{j}^{n})_{k}=\mathcal{E}_{j,k-1}^{n}-2\mathcal{E}_{j,k}^{n}+\mathcal{E}_{j,k+1}^{n},$

which can be seen at a discrete approximation of
$\partial_{y}^{2}\mathcal{E}(t_{n},x_{j},y_{k})$ up to a scaling factor of
order $\Delta y^{2}$. Once again, setting $\nu=\frac{\Delta x}{\Delta t}$ and
introducing $\kappa=\frac{\Delta y^{2}}{\Delta t}$, we may rewrite (2.1) with
$\mathcal{W}^{f}=\mathcal{W}^{b}=\mathbf{I}_{d}+\xi\mathbf{A}$ as

$\displaystyle(1-\beta)\frac{\mathcal{E}_{j}^{n+1}-\mathcal{E}_{j}^{n}}{\Delta
t}$
$\displaystyle=\beta\nu\frac{\mathcal{E}_{j-1}^{n+1}-\mathcal{E}_{j}^{n+1}}{\Delta
x}+\lambda\nu\frac{\mathcal{E}_{j+1}^{n}-\mathcal{E}_{j}^{n}}{\Delta
x}-\alpha\nu\frac{\mathcal{E}_{j}^{n}-\mathcal{E}_{j-1}^{n}}{\Delta x}$
$\displaystyle~{}~{}~{}+\beta\xi\kappa\frac{\mathbf{A}}{\Delta
y^{2}}\mathcal{E}_{j-1}^{n+1}+\alpha\xi\kappa\frac{\mathbf{A}}{\Delta
y^{2}}\mathcal{E}_{j-1}^{n}-2\alpha\xi\kappa\frac{\mathbf{A}}{\Delta
y^{2}}\mathcal{E}_{j}^{n}-\alpha\xi^{2}\kappa\Delta
y^{2}\frac{\mathbf{A}}{\Delta y^{2}}\frac{\mathbf{A}}{\Delta
y^{2}}\mathcal{E}_{j}^{n}$
$\displaystyle~{}~{}~{}+\lambda\xi\kappa\frac{\mathbf{A}}{\Delta
y^{2}}\mathcal{E}_{j+1}^{n}.$

Now letting $\Delta t\to 0$, $\Delta x\to 0$ and $\Delta y\to 0$ with $\nu$
and $\kappa$ fixed, we obtain the following partial differential equation

$\partial_{t}\mathcal{E}(t,x,y)+\frac{\nu(\beta+\alpha-\lambda)}{1-\beta}\partial_{x}\mathcal{E}(t,x,y)=\frac{\xi\kappa(\beta+\lambda-\alpha)}{1-\beta}\partial_{y}^{2}\mathcal{E}(t,x,y).$

This is a diffusion equation along the $y$ dimension while being a transport
in the $x$ direction. As such, it is only well defined (or stable) when the
sign of the diffusion coefficient in front of
$\partial_{y}^{2}\mathcal{E}(t,x,y)$ is positive. This depends on the sign of
$\xi$ and $\beta+\lambda-\alpha$, which need to verify
$\xi(\beta+\lambda-\alpha)>0$. In that case, the system diffuses neural
activity along the dimension $y$ such that the entire neuronal layer converges
to a single, uniform activation value when $t\to\infty$.

###  4.2 The general symmetric case

Finally, we now wish to relax some of the assumptions made in the previous
Rao-Ballard case. Thus, the last case that we present is one where we assume
that

  * (i)

$\mathcal{W}^{f}$ and $\mathcal{W}^{b}$ are symmetric matrices, that is
$\mathcal{W}^{f},\mathcal{W}^{b}\in\mathscr{S}_{d}(\mathbb{R})$,

  * (ii)

$\mathcal{W}^{f}$ and $\mathcal{W}^{b}$ commute, that is
$\mathcal{W}^{f}\mathcal{W}^{b}=\mathcal{W}^{b}\mathcal{W}^{f}$.

But we do not necessarily impose that $\mathcal{W}^{f}=(\mathcal{W}^{b})^{\bf
t}$ as in the Rao & Ballard’s previous case. Let us already note that examples
of matrices verifying the above conditions are residual convolution matrices
introduced in (4.3), that is
$\mathcal{W}^{f}=\zeta_{f}\mathbf{I}_{d}+\xi_{f}\mathbf{A}$ and
$\mathcal{W}^{b}=\zeta_{b}\mathbf{I}_{d}+\xi_{b}\mathbf{A}$ for some
$\zeta_{b,f},\xi_{b,f}\in\mathbb{R}$. Under assumptions (i) and (ii),
$\mathcal{W}^{f}$ and $\mathcal{W}^{b}$ can be diagonalized in the same
orthonormal basis, meaning that there exist an invertible orthogonal matrix
$P\in\mathscr{M}_{d}(\mathbb{R})$ such that
$PP^{\mathbf{t}}=P^{\mathbf{t}}P=I_{d}$, and two diagonal matrices
$\mathcal{D}^{f}\in\mathscr{M}_{d}(\mathbb{R})$ and
$\mathcal{D}^{b}\in\mathscr{M}_{d}(\mathbb{R})$ with the properties that

$P^{\mathbf{t}}\mathcal{W}^{f}P=\mathcal{D}^{f},\quad\text{ and }\quad
P^{\mathbf{t}}\mathcal{W}^{b}P=\mathcal{D}^{b}.$

For future reference, we denote by $\gamma_{p}^{f,b}$ for each $1\leq p\leq d$
the diagonal elements of $\mathcal{D}^{f,b}$. Once again, we can use the
matrix $P$ to apply an orthonormal basis change and create neural asssemblies
$\mathcal{U}_{j}^{n}:=P^{t}\mathcal{E}_{j}^{n}$. With
$P\mathcal{U}_{j}^{n}:=\mathcal{E}_{j}^{n}$, the recurrence equation becomes

$\mathcal{U}_{j}^{n+1}-\beta\mathcal{D}^{f}\mathcal{U}_{j-1}^{n+1}=\alpha\mathcal{D}^{b}\mathcal{U}_{j-1}^{n}+\left[(1-\beta-\lambda)\mathbf{I}_{d}-\alpha\mathcal{D}^{b}\mathcal{D}^{b}\right]\mathcal{U}_{j}^{n}+\lambda\mathcal{D}^{b}\mathcal{U}_{j+1}^{n}.$

Note that, because all matrices in the above equation are diagonal, we have
also totally decoupled the $d$ components of the vector $\mathcal{U}_{j}^{n}$.
More precisely, by denoting $u_{j,p}^{n}$ the $p$th component of
$\mathcal{U}_{j}^{n}$, that is
$\mathcal{U}_{j}^{n}=(u_{j,1}^{n},\cdots,u_{j,d}^{n})^{\mathbf{t}}$, we obtain

$u_{j,p}^{n+1}-\beta\gamma_{p}^{f}u_{j-1,p}^{n+1}=\alpha\gamma_{p}^{b}u_{j-1,p}^{n}+(1-\beta-\lambda-\alpha\left(\gamma_{p}^{b}\right)^{2})u_{j,p}^{n}+\lambda\gamma_{p}^{b}u_{j+1,p}^{n},\quad
p=1,\cdots,d.$

This indicates that one needs to study

$u_{j}^{n+1}-\beta\gamma_{1}u_{j-1}^{n+1}=\alpha\gamma_{2}u_{j-1}^{n}+(1-\beta-\lambda-\alpha\gamma_{2}^{2})u_{j}^{n}+\lambda\gamma_{2}u_{j+1}^{n},$
(4.4)

where $\gamma_{1,2}\in\mathbb{R}$ are now two given parameters. As before,
$\gamma_{1,2}$ can be thought of as the connection strength across layers of
the neural assembly under consideration. By construction, each assembly in a
given layer is only connected to the corresponding assembly in the layer
above, and similarly in the layer below, with $\gamma_{1}$ for the feedforward
direction and $\gamma_{2}$ for the feedback direction. Note that
$\gamma_{1}=\gamma_{2}$ would then correspond to the Rao-Ballard situation
studied previously.

####  4.2.1 Study of the amplification factor function

Repeating the previous analysis, one needs to understand the amplification
factor

$\rho_{\gamma_{1},\gamma_{2}}(\theta)=\frac{\alpha\gamma_{2}\left(e^{-\mathbf{i}\theta}-\gamma_{2}\right)+1-\beta+\lambda\left(\gamma_{2}e^{\mathbf{i}\theta}-1\right)}{1-\beta\gamma_{1}e^{-\mathbf{i}\theta}},\quad\theta\in[-\pi,\pi].$

We already note a symmetry property of the amplification factor function which
reads

$\rho_{\gamma_{1},\gamma_{2}}(\theta)=\rho_{-\gamma_{1},-\gamma_{2}}(\theta\pm\pi),\quad\theta\in[-\pi,\pi].$

As a consequence, whenever $\rho_{\gamma_{1},\gamma_{2}}(0)=\pm 1$ one has
$\rho_{-\gamma_{1},-\gamma_{2}}(\pm\pi)=\pm 1$ for the same values of the
parameters. Then, we note that

$\rho_{\gamma_{1},\gamma_{2}}(0)=1\Longleftrightarrow\gamma_{1}=\chi(\gamma_{2}),$

where the function $\chi(x)$, depending only on the hyper-parameters, is given
by

$\chi(x):=\frac{\alpha
x^{2}-(\alpha+\lambda)x+\lambda+\beta}{\beta},x\in\mathbb{R}.$ (4.5)

Thus, using the above symmetry, we readily deduce that

$\rho_{\gamma_{1},\gamma_{2}}(\pm\pi)=1\Longleftrightarrow\gamma_{1}=-\chi(-\gamma_{2}).$

Finally, we compute that

$\rho_{\gamma_{1},\gamma_{2}}(0)=-1\Longleftrightarrow\gamma_{1}=\zeta(\gamma_{2}),$

where the function $\zeta(x)$, depending only on the hyper-parameters, is
given by

$\zeta(x):=\frac{-\alpha
x^{2}+(\alpha+\lambda)x+2-\lambda-\beta}{\beta},x\in\mathbb{R}.$ (4.6)

Using the above symmetry, we readily deduce that

$\rho_{\gamma_{1},\gamma_{2}}(\pm\pi)=-1\Longleftrightarrow\gamma_{1}=-\zeta(-\gamma_{2}).$

A complete and exhaustive characterization of all possible cases as a function
of $\gamma_{1,2}$ and the hyper-parameters is beyond the scope of this paper.
Nevertheless, we can make some few further remarks. The four above curves
$\gamma_{1}=\chi(\gamma_{2})$, $\gamma_{1}=-\chi(-\gamma_{2})$,
$\gamma_{1}=\zeta(\gamma_{2})$ and $\gamma_{1}=-\zeta(-\gamma_{2})$ form
parabolas in the plane $(\gamma_{1},\gamma_{2})$ that can intersect and
provide the boundaries of the stability regions. For example, we can notice
that $\gamma_{1}=\zeta(\gamma_{2})$ and $\gamma_{1}=-\zeta(-\gamma_{2})$
intersect if and only if $\gamma_{2}=\pm\sqrt{\frac{2-\lambda-\beta}{\alpha}}$
whereas $\gamma_{1}=\chi(\gamma_{2})$ and $\gamma_{1}=-\chi(-\gamma_{2})$ can
never intersect. We refer to Figure 14 for an illustration of the stability
regions and their boundaries in the case
$(\alpha,\beta,\lambda)=(0.4,0.2,0.3)$. Here, we see that stability can be
obtained with large values of the feedforward connection strength
$\gamma_{1}$, but this requires the feedback connections strength $\gamma_{2}$
to remain low. Of course, different qualitative behaviors and stability
regions may be obtained for different choices of the hyperparameters
$(\alpha,\beta,\lambda)$; while it is beyond the scope of the present study to
characterize them all, it is important to point out that such a
characterization is feasible using the present method, for any choice of the
hyperparameters.

More interestingly, we can investigate the dependence of the wave speed as a
function of the parameters $\gamma_{1}$ and $\gamma_{2}$. For example, when
$\gamma_{1}=\chi(\gamma_{2})$, we have that

$\rho_{\chi(\gamma_{2}),\gamma_{2}}(\theta)=\exp\left(-\mathbf{i}\frac{(\alpha-\lambda)\gamma_{2}+\beta\chi(\gamma_{2})}{1-\beta\chi(\gamma_{2})}\theta+\mathcal{O}(|\theta|^{2})\right),\text{
as }\theta\rightarrow 0,$

such that the associated wave speed is given by

$c_{0}^{\chi}=\frac{(\alpha-\lambda)\gamma_{2}+\beta\chi(\gamma_{2})}{1-\beta\chi(\gamma_{2})},$

whose sign may vary as $\gamma_{2}$ is varied. We refer to the forthcoming
section 4.2.4 below for a practical example (see Figure 18).

Figure 14: Stability/instability regions and their boundaries as a function of
$(\gamma_{1},\gamma_{2})$ for (4.4) for fixed values of the hyperparameters
$(\alpha,\beta,\lambda)$. The shaded orange region corresponds to an
instability for (4.4) while the purple region corresponds to a stability for
(4.4). The boundaries of the stability/instability regions are given by the
intersections of the parametrized curves $\gamma_{1}=\chi(\gamma_{2})$ (blue
curve), $\gamma_{1}=-\chi(-\gamma_{2})$ (light blue curve),
$\gamma_{1}=\zeta(\gamma_{2})$ (dark green curves) and
$\gamma_{1}=-\zeta(-\gamma_{2})$ (light green curves) where equation (4.4) is
marginally stable. We represented the line $\gamma_{1}=\gamma_{2}$ (black
curve) which corresponds to the case studied in Figure 11(d) with
$(\beta,\lambda)$ in Region $\mathbf{(III)}$ and $\alpha$ fixed in
$(\Lambda,\lambda+\beta)$.

####  4.2.2 Continuous in time interpretation

As done in previous sections, we now perform a continuous in time limit of the
model (4.4). With the same scaling on the hyperparameters

$\widetilde{\beta}:=\frac{\beta}{\Delta
t},\quad\widetilde{\lambda}:=\frac{\lambda}{\Delta t},\text{ and
}\widetilde{\alpha}:=\frac{\alpha}{\Delta t},$

we get that, at the limit $\Delta t\rightarrow 0$, the equation (4.4) becomes
the following lattice ordinary differential equation

$\frac{\mathrm{d}}{\mathrm{d}t}{\bf
u}_{j}(t)=(\widetilde{\beta}\gamma_{1}+\widetilde{\alpha}\gamma_{2}){\bf
u}_{j-1}(t)-(\widetilde{\beta}+\widetilde{\lambda}+\widetilde{\alpha}\gamma_{2}^{2}){\bf
u}_{j}(t)+\widetilde{\lambda}\gamma_{2}{\bf u}_{j+1}(t),\quad t>0.$ (4.7)

The amplification factor function in this case is given by

$\nu_{\gamma_{1},\gamma_{2}}(\theta)=(\widetilde{\beta}\gamma_{1}+\widetilde{\alpha}\gamma_{2})e^{-\mathbf{i}\theta}-(\widetilde{\beta}+\widetilde{\lambda}+\widetilde{\alpha}\gamma_{2}^{2})+\widetilde{\lambda}\gamma_{2}e^{\mathbf{i}\theta},\quad\theta\in[-\pi,\pi],$

whose real part is given by

$\mathrm{Re}(\nu_{\gamma_{1},\gamma_{2}}(\theta))=(\widetilde{\beta}\gamma_{1}+(\widetilde{\alpha}+\widetilde{\lambda})\gamma_{2})\cos(\theta)-(\widetilde{\beta}+\widetilde{\lambda}+\widetilde{\alpha}\gamma_{2}^{2}),\quad\theta\in[-\pi,\pi].$

When
$\widetilde{\beta}\gamma_{1}+(\widetilde{\alpha}+\widetilde{\lambda})\gamma_{2}>0$,
we observe that

$\underset{\theta\in[-\pi,\pi]}{\max}\mathrm{Re}(\nu_{\gamma_{1},\gamma_{2}}(\theta))=\mathrm{Re}(\nu_{\gamma_{1},\gamma_{2}}(0))=\widetilde{\beta}\gamma_{1}+(\widetilde{\alpha}+\widetilde{\lambda})\gamma_{2}-(\widetilde{\beta}+\widetilde{\lambda}+\widetilde{\alpha}\gamma_{2}^{2}),$

such that

$\underset{\theta\in[-\pi,\pi]}{\max}\mathrm{Re}(\nu_{\gamma_{1},\gamma_{2}}(\theta))=0\Longleftrightarrow\gamma_{1}=\frac{\widetilde{\alpha}\gamma_{2}^{2}-(\widetilde{\alpha}+\widetilde{\lambda})\gamma_{2}+\widetilde{\beta}+\widetilde{\lambda}}{\widetilde{\beta}}.$

Whereas, when
$\widetilde{\beta}\gamma_{1}+(\widetilde{\alpha}+\widetilde{\lambda})\gamma_{2}<0$,
we observe that

$\underset{\theta\in[-\pi,\pi]}{\max}\mathrm{Re}(\nu_{\gamma_{1},\gamma_{2}}(\theta))=\mathrm{Re}(\nu_{\gamma_{1},\gamma_{2}}(\pm\pi))=-\widetilde{\beta}\gamma_{1}-(\widetilde{\alpha}+\widetilde{\lambda})\gamma_{2}-(\widetilde{\beta}+\widetilde{\lambda}+\widetilde{\alpha}\gamma_{2}^{2}),$

such that

$\underset{\theta\in[-\pi,\pi]}{\max}\mathrm{Re}(\nu_{\gamma_{1},\gamma_{2}}(\theta))=0\Longleftrightarrow\gamma_{1}=\frac{-\widetilde{\alpha}\gamma_{2}^{2}-(\widetilde{\alpha}+\widetilde{\lambda})\gamma_{2}-\widetilde{\beta}-\widetilde{\lambda}}{\widetilde{\beta}}.$

As a consequence, the stability regions are determined by the locations of the
parabolas
$\gamma_{2}\mapsto\frac{\widetilde{\alpha}\gamma_{2}^{2}-(\widetilde{\alpha}+\widetilde{\lambda})\gamma_{2}+\widetilde{\beta}+\widetilde{\lambda}}{\widetilde{\beta}}$
and
$\gamma_{2}\mapsto\frac{-\widetilde{\alpha}\gamma_{2}^{2}-(\widetilde{\alpha}+\widetilde{\lambda})\gamma_{2}-\widetilde{\beta}-\widetilde{\lambda}}{\widetilde{\beta}}$
in the plane $(\gamma_{1},\gamma_{2})$. We observe that they never intersect
and are oriented in the opposite directions and refer to Figure 15 for a
typical configuration. Here, we see that the system is stable for a very large
range of values of both $\gamma_{1}$ and $\gamma_{2}$. In particular, for
large enough values of the feedback connection weight (e.g. $|\gamma_{2}|>3$),
stability is guaranteed regardless of the value of the feedforward connection
weight $\gamma_{1}$ (within a reasonable range, e.g. $\gamma_{1}\in(-10,10)$).
This is the opposite behavior as that obtained for the discrete system in
Figure 14, where stability was impossible under the same conditions for
$\gamma_{1,2}$. This highlights again the errors of interpretation that can
potentially be caused by discrete approximation of a continuous system.

Figure 15: Stability/instability regions and their boundaries as a function of
$(\gamma_{1},\gamma_{2})$ for (4.7) for any
$(\widetilde{\alpha},\widetilde{\lambda},\widetilde{\beta})$ fixed with
$\widetilde{\alpha}=\widetilde{\lambda}+\widetilde{\beta}$. The shaded orange
region corresponds to an instability for (4.7) while the purple region
corresponds to a stability for (4.7). The boundaries of the
stability/instability regions are given by the intersections of the parabolas
$\gamma_{1}=\frac{\widetilde{\alpha}\gamma_{2}^{2}-(\widetilde{\alpha}+\widetilde{\lambda})\gamma_{2}+\widetilde{\beta}+\widetilde{\lambda}}{\widetilde{\beta}}$
and
$\gamma_{1}=\frac{-\widetilde{\alpha}\gamma_{2}^{2}-(\widetilde{\alpha}+\widetilde{\lambda})\gamma_{2}-\widetilde{\beta}-\widetilde{\lambda}}{\widetilde{\beta}}$
(magenta curves) where equation (4.7) is marginally stable. We represented the
line $\gamma_{1}=\gamma_{2}$ (blue curve) which corresponds to the case
studied in Figure 13.

####  4.2.3 Fully continuous interpretation when
$\mathcal{W}^{f}=\mathbf{I}_{d}+\xi_{f}\mathbf{A}$ and
$\mathcal{W}^{b}=\mathbf{I}_{d}+\xi_{b}\mathbf{A}$.

When $\mathcal{W}^{f}=\mathbf{I}_{d}+\xi_{f}\mathbf{A}$ and
$\mathcal{W}^{b}=\mathbf{I}_{d}+\xi_{b}\mathbf{A}$, one can once again
identify $\mathcal{E}_{j,k}^{t}$ as the approximation of some smooth function
$\mathcal{E}(t,x,y)$ at $t_{n}=n\Delta t$, $x_{j}=j\Delta x$ and
$y_{k}=k\Delta y$, along the three dimensions of time, network depth and
neuronal layer width. We may rewrite (2.1) in this case as

$\displaystyle(1-\beta)\frac{\mathcal{E}_{j}^{n+1}-\mathcal{E}_{j}^{n}}{\Delta
t}$
$\displaystyle=\beta\nu\frac{\mathcal{E}_{j-1}^{n+1}-\mathcal{E}_{j}^{n+1}}{\Delta
x}+\lambda\nu\frac{\mathcal{E}_{j+1}^{n}-\mathcal{E}_{j}^{n}}{\Delta
x}-\alpha\nu\frac{\mathcal{E}_{j}^{n}-\mathcal{E}_{j-1}^{n}}{\Delta x}$
$\displaystyle~{}~{}~{}+\beta\xi_{f}\kappa\frac{\mathbf{A}}{\Delta
y^{2}}\mathcal{E}_{j-1}^{n+1}+\alpha\xi_{b}\kappa\frac{\mathbf{A}}{\Delta
y^{2}}\mathcal{E}_{j-1}^{n}-2\alpha\xi_{b}\kappa\frac{\mathbf{A}}{\Delta
y^{2}}\mathcal{E}_{j}^{n}-\alpha\xi_{b}^{2}\kappa\Delta
y^{2}\frac{\mathbf{A}}{\Delta y^{2}}\frac{\mathbf{A}}{\Delta
y^{2}}\mathcal{E}_{j}^{n}$
$\displaystyle~{}~{}~{}+\lambda\xi_{b}\kappa\frac{\mathbf{A}}{\Delta
y^{2}}\mathcal{E}_{j+1}^{n},$

such that in the limit $\Delta t\to 0$, $\Delta x\to 0$ and $\Delta y\to 0$
with $\nu$ and $\kappa$ fixed, we obtain the following partial differential
equation

$\partial_{t}\mathcal{E}(t,x,y)+\frac{\nu(\beta+\alpha-\lambda)}{1-\beta}\partial_{x}\mathcal{E}(t,x,y)=\kappa\frac{\beta\xi_{f}+(\lambda-\alpha)\xi_{b}}{1-\beta}\partial_{y}^{2}\mathcal{E}(t,x,y).$

As before, this is a diffusion equation along the $y$ dimension, whose
stability depends on the positivity of the diffusion coefficient, i.e.
$\beta\xi_{f}+(\lambda-\alpha)\xi_{b}\geq 0$ .

####  4.2.4 Application to a ring model of orientations

Going back to our discrete system, in this section we consider the case where
neurons within each layer encode for a given orientation in $[0,\pi]$. Here,
we have in mind visual stimuli which are made of a fixed elongated black bar
on a white background with a prescribed orientation. We introduce the
following matrix $\mathbf{A}_{\mathrm{per}}\in\mathscr{M}_{d}(\mathbb{R})$
given by

$\mathbf{A}_{\mathrm{per}}=\left(\begin{matrix}-2&1&0&\cdots&0&1\\\
1&-2&1&\ddots&\ddots&0\\\ 0&\ddots&\ddots&\ddots&\ddots&\vdots\\\
\vdots&\ddots&\ddots&\ddots&\ddots&\vdots\\\ 0&\ddots&\ddots&1&-2&1\\\
1&0&\cdots&0&1&-2\end{matrix}\right),$

which is nothing but the discretizing of the Laplacian with boundary
condition. Indeed, for each $\mathcal{E}_{j}^{n}\in\mathbb{R}^{d}$, we assume
that neuron $\mathcal{E}_{j,k}^{n}$ encodes for orientation $\frac{k}{d}\pi$
for $k=1,\cdots,d$. We readily remark that
$0\in\mathrm{Spec}(\mathbf{A}_{\mathrm{per}})$ with corresponding eigenvector
$U_{1}=(1,\cdots,1)^{\bf t}\in\mathbb{R}^{d}$. Furthermore, we have:

  * •

if $d=2m+1$ is odd, then $\lambda_{p}=-4\sin\left(\frac{p\pi}{d}\right)^{2}$
with $p=1,\cdots,m$ is an eigenvalue of $\mathbf{A}_{\mathrm{per}}$ of
multiplicity $2$ with associated eigenvectors

$\displaystyle U_{2p}$
$\displaystyle=\left(\cos\left(\frac{2p\pi}{d}\right),\cdots,\cos\left(\frac{2kp\pi}{d}\right),\cdots,1\right)^{\bf
t}\in\mathbb{R}^{d},$ $\displaystyle U_{2p+1}$
$\displaystyle=\left(\sin\left(\frac{2p\pi}{d}\right),\cdots,\sin\left(\frac{2kp\pi}{d}\right),\cdots,0\right)^{\bf
t}\in\mathbb{R}^{d};$

  * •

if $d=2m$ is even, then $\lambda_{p}=-4\sin\left(\frac{p\pi}{d}\right)^{2}$
with $p=1,\cdots,m-1$ is an eigenvalue of $\mathbf{A}_{\mathrm{per}}$ of
multiplicity $2$ with associated eigenvectors $U_{2p}$ and $U_{2p+1}$ as
above. And $\lambda=-4$ is a simple eigenvalue of $\mathbf{A}_{\mathrm{per}}$
with associated eigenvector $U_{d}=(-1,1,-1,1,\cdots,-1,1)\in\mathbb{R}^{d}$.

(a) Eigenvectors $U_{2}$ and $U_{3}$. (b) Eigenvectors $U_{4}$ and $U_{5}$.

Figure 16: We plot the eigenvectors $U_{2p}$ and $U_{2p+1}$ for $p=1$ and
$p=2$ as a function of $\frac{k}{d}\pi$ for $k=1,\cdots,d$. We note that
$U_{2p}$ and $U_{2p+1}$ encode the first Fourier modes. Here we have set
$d=2^{5}$.

It may be interesting to note that any linear combinations of $U_{2p}$ and
$U_{2p+1}$ can always be written in the form

$aU_{2p}+bU_{2p+1}=A\left(\cos\left(\frac{2p\pi}{d}+\varphi\right),\cdots,\cos\left(\frac{2kp\pi}{d}+\varphi\right),\cdots,\cos(\varphi)\right)^{\bf
t}\in\mathbb{R}^{d},$

where $A=\sqrt{a^{2}+b^{2}}>0$ and
$\varphi=-\mathrm{arctan}\left(\frac{b}{a}\right)\in(-\pi/2,\pi/2)$ whenever
$a\neq 0$ and $b\neq 0$. This means that $U_{2p}$ and $U_{2p+1}$ span all
possible translations modulo $[0,\pi]$ of a fixed profile. We refer to Figure
16 for a visualization of the first eigenvectors. In short, these eigenvectors
$U_{i}$ implement a Fourier transform of the matrix
$\mathbf{A}_{\mathrm{per}}$.

We now set $\mathcal{W}^{b}$ to be

$\mathcal{W}^{b}=\frac{1}{2}\mathbf{I}_{d}-\frac{1}{4}\mathbf{A}_{\mathrm{per}}=\left(\begin{matrix}1&-\frac{1}{4}&0&\cdots&0&-\frac{1}{4}\\\
-\frac{1}{4}&1&-\frac{1}{4}&\ddots&\ddots&0\\\
0&\ddots&\ddots&\ddots&\ddots&\vdots\\\
\vdots&\ddots&\ddots&\ddots&\ddots&\vdots\\\
0&\ddots&\ddots&-\frac{1}{4}&1&-\frac{1}{4}\\\
-\frac{1}{4}&0&\cdots&0&-\frac{1}{4}&1\end{matrix}\right),$

which means that $\mathcal{W}^{b}$ acts as a convolution with local excitation
and lateral inhibition. From now on, to fix ideas, we will assume that $d=2m$
is even. We define the following matrix

$P=\left(U_{1},U_{2},\cdots,U_{d}\right)\in\mathscr{M}_{d}(\mathbb{R}).$

As a consequence, we have the decomposition

$P^{\mathbf{t}}\mathcal{W}^{b}P=\mathcal{D}^{b},$

with
$\mathcal{D}^{b}=\mathrm{diag}\left(\frac{1}{2},\frac{1}{2}-\frac{1}{4}\lambda_{1},\frac{1}{2}-\frac{1}{4}\lambda_{1},\cdots,\frac{1}{2}-\frac{1}{4}\lambda_{m-1},\frac{1}{2}-\frac{1}{4}\lambda_{m-1},\frac{3}{2}\right)\in\mathscr{M}_{d}(\mathbb{R})$.
Now, for given values of the hyper-parameters $(\alpha,\beta,\lambda)$ with
$\beta>0$, we set $\mathcal{D}^{f}:=\chi(\mathcal{D}^{b})$ where the map
$\chi$, defined in (4.5), is applied to the diagonal elements of
$\mathcal{D}^{b}$, that is

$\mathcal{D}^{f}=\mathrm{diag}\left(\chi\left(\frac{1}{2}\right),\chi\left(\frac{1}{2}-\frac{1}{4}\lambda_{1}\right),\cdots,\chi\left(\frac{1}{2}-\frac{1}{4}\lambda_{m-1}\right),\chi\left(\frac{3}{2}\right)\right)\in\mathscr{M}_{d}(\mathbb{R}).$

And then we set $\mathcal{W}^{f}:=P\mathcal{D}^{f}P^{\bf t}$. We refer to
Figure 17 for an illustration of the structures of matrices $\mathcal{W}^{f}$
and $\mathcal{W}^{b}$. For large set of values of the hyper-parameters,
$\mathcal{W}^{f}$ still present a band structure with positive elements on the
diagonals indicating that $\mathcal{W}^{f}$ can also be interpreted as a
convolution with local excitation. For the values of the hyper-parameters
fixed in Figure 17, the feedforward matrix $\mathcal{W}^{f}$ is purely
excitatory.

(a) $\mathcal{W}^{f}$. (b) $\mathcal{W}^{b}$.

Figure 17: Illustration of the matrices $\mathcal{W}^{f}$ and
$\mathcal{W}^{b}$ for $d=2^{5}$ neurons and values of the hyper-parameters
fixed to $(\alpha,\beta,\lambda)=(0.1,0.1,0.5)$. Note the band structure of
$\mathcal{W}^{f}$ with local excitation.

Reproducing the analysis developed in the previous Subsection 4.2, we perform
a change of orthonormal basis to express neural activities in terms of the
relevant assemblies $\mathcal{U}_{j}^{n}:=P^{t}\mathcal{E}_{j}^{n}$. With
$P\mathcal{U}_{j}^{n}:=\mathcal{E}_{j}^{n}$, the recurrence equation becomes

$\mathcal{U}_{j}^{n+1}-\beta\mathcal{D}^{f}\mathcal{U}_{j-1}^{n+1}=\alpha\mathcal{D}^{b}\mathcal{U}_{j-1}^{n}+\left[(1-\beta-\lambda)\mathbf{I}_{d}-\alpha\mathcal{D}^{b}\mathcal{D}^{b}\right]\mathcal{U}_{j}^{n}+\lambda\mathcal{D}^{b}\mathcal{U}_{j+1}^{n}.$

Then, if we denote by $\gamma_{p}$ the $p$th diagonal element of
$\mathcal{D}^{b}$, then for each $p=1,\cdots,d$ the above recurrence writes

$u_{j,p}^{n+1}-\beta\chi(\gamma_{p})u_{j-1,p}^{n+1}=\alpha\gamma_{p}u_{j-1,p}^{n}+(1-\beta-\lambda-\alpha\gamma_{p}^{2})u_{j,p}^{n}+\lambda\gamma_{p}u_{j+1,p}^{n},$

where $u_{j,p}^{n}$ is the $p$th component (or neural assembly) of
$\mathcal{U}_{j}^{n}$. For each $p=1,\cdots,d$, the associated amplification
factor function reads

$\rho_{p}(\theta)=\frac{\alpha\gamma_{p}\left(e^{-\mathbf{i}\theta}-\gamma_{p}\right)+1-\beta+\lambda\left(\gamma_{p}e^{\mathbf{i}\theta}-1\right)}{1-\beta\chi(\gamma_{p})e^{-\mathbf{i}\theta}},\quad\theta\in[-\pi,\pi],$

and with our specific choice of function $\chi$, we have that $\rho_{p}(0)=1$
with

$\rho_{p}(\theta)=\exp\left(-\mathbf{i}\frac{(\alpha-\lambda)\gamma_{p}+\beta\chi(\gamma_{p})}{1-\beta\chi(\gamma_{p})}\theta-\sigma_{0}^{p}\theta^{2}+\mathcal{O}(|\theta|^{3})\right),\text{
as }\theta\rightarrow 0,$

such that the associated wave speed is given by

$c_{0}^{p}=\frac{(\alpha-\lambda)\gamma_{p}+\beta\chi(\gamma_{p})}{1-\beta\chi(\gamma_{p})},$

and where we have set

$\sigma_{0}^{p}=\frac{\alpha(1+4\lambda)\gamma_{p}^{2}+\beta+\lambda-(\alpha+\lambda)\gamma_{p}(\alpha\gamma_{p}^{2}+\beta+\lambda)}{2(1-\beta\chi(\gamma_{p}))^{2}}.$

From now on, we assume that we have tuned the hyper-parameters such that
$|\rho_{p}(\theta)|<1$ for all
$\theta\in[-\pi,\pi]\backslash\left\\{0\right\\}$ and each $p=1,\cdots,d$.
This can in fact be systematically checked numerically for a given set of
hyper-parameters. We report in Figure 18 the shape of $p\mapsto c_{0}^{p}$ for
the same values of the hyper-parameters as the ones in Figure 17 and
$d=2^{5}$. We first remark that $p\mapsto c_{0}^{p}$ is a monotone decreasing
map, and in our specific case we have

$c^{d}_{0}<c^{d-1}_{0}=c^{d-2}_{0}<\cdots<c^{9}_{0}=c^{8}_{0}<0<c^{7}_{0}=c^{6}_{0}<c^{5}_{0}=c^{4}_{0}<c^{3}_{0}=c^{2}_{0}<c^{1}_{0}.$

Figure 18: Plot of the wave speed $c_{0}^{p}$ for $p=1,\cdots,d$ (colored
dots). The color code (blue/red) refers to the sign of $c_{0}^{p}$: blue when
positive and dark red when negative. Note that only the elements associated to
the eigenvalues $\frac{1}{2}$, $\frac{1}{2}-\frac{1}{4}\lambda_{1}$,
$\frac{1}{2}-\frac{1}{4}\lambda_{2}$ and $\frac{1}{2}-\frac{1}{4}\lambda_{3}$
are positive. We also remark that $c_{0}^{p}$ is a monotone decreasing
function. Here $d=2^{5}$ and values of the hyper-parameters are fixed to
$(\alpha,\beta,\lambda)=(0.1,0.1,0.5)$.

Given a fixed input entry $\mathcal{E}_{0}\in\mathbb{R}^{d}$ presented at
$j=0$ to the network continually at each time step, we can deduce which
components of $\mathcal{E}_{0}\in\mathbb{R}^{d}$ will be able to propagate
forward through the network. More precisely, we can decompose
$\mathcal{E}_{0}$ along the basis $\left(U_{1},\cdots U_{d}\right)$ of
eigenvectors, that is

$\mathcal{E}_{0}=\sum_{p=1}^{d}a_{p}U_{p},$

for some real coefficients $a_{p}$ for $p=1,\cdots,d$. Assuming that the
network was at rest initially, we get that the dynamics along each eigenvector
(or neural assembly) is given by

$\left\\{\begin{split}u_{j,p}^{n+1}-\beta\chi(\gamma_{p})u_{j-1,p}^{n+1}&=\alpha\gamma_{p}u_{j-1,p}^{n}+(1-\beta-\lambda-\alpha\gamma_{p}^{2})u_{j,p}^{n}+\lambda\gamma_{p}u_{j+1,p}^{n},\quad
j\geq 1,\quad n\geq 0,\\\ u_{0,p}^{n}&=a_{p},\quad n\geq 0,\\\
u_{j,p}^{0}&=0,\quad j\geq 1.\end{split}\right.$ (4.8)

Thus, we readily obtain that

$\mathcal{E}_{j}^{n}=\sum_{p=1}^{d}u_{j,p}^{n}U_{p},\quad j\geq 1,\quad n\geq
1,$

where $u_{j,p}^{n}$ is a solution to (4.8).

As a consequence, the monotonicity property of the map $p\mapsto c_{0}^{p}$
indicates that the homogeneous constant mode $U_{1}$ is the fastest to
propagate forward into the network with associated spreading speed
$c^{1}_{0}$, it is then followed by the modes $(U_{2},U_{3})$ propagating at
speed $c^{2}_{0}=c^{3}_{0}$. In our numerics, we have set the parameters such
that $c^{1}_{0}\approx c^{2}_{0}=c^{3}_{0}$ with a significant gap with the
other wave speeds. Lets us remark, that all modes $U_{p}$ with $p\geq 8$ are
not able to propagate into the network (see Figure 19). Thus our architecture
acts as a mode filter.

(a) Case $p=2$ with $\gamma_{2}=\frac{1}{2}-\frac{1}{4}\lambda_{1}\sim
0.5096$. (b) Case $p=10$ with
$\gamma_{10}=\frac{1}{2}-\frac{1}{4}\lambda_{5}\sim 0.7222$.

Figure 19: Space-time plot of the solution of the recurrence equation for
(4.8) for $p=2$ and $p=10$ associated to respectively positive wave speed
$c_{0}^{2}>0$ and negative wave speed $c_{0}^{10}<0$. The neural assembly
associated with the 2nd eigenvector of the connectivity matrix propagates its
input signal into the network at constant speed; but the neural assembly
associated with the 10th eigenvector does not propagate the signals it
receives on the input layer.

Even more precisely, let us remark that the sequence
$\left(a_{p}\left(\frac{\alpha\gamma_{p}+\beta\chi(\gamma_{p})}{\lambda\gamma_{p}}\right)^{j}\right)_{j\geq
0}$ is a stationary solutions of (4.8) which remains bounded whenever $p$ is
such that the associated wave speed is negative, that is $c_{0}^{p}<0$, since
in that case, one has
$\alpha\gamma_{p}+\beta\chi(\gamma_{p})<\lambda\gamma_{p}$. The solution
$\mathcal{E}_{j}^{n}$ can then be approximated as

$\mathcal{E}_{j}^{n}\simeq\sum_{p~{}:~{}c_{0}^{p}>0}\frac{a_{p}}{2}\left(1-\mathrm{erf}\left(\frac{j-c_{0}^{p}n}{\sqrt{4\sigma_{0}^{p}n}}\right)\right)U_{p}+\sum_{p~{}:~{}c_{0}^{p}<0}a_{p}\left(\frac{\alpha\gamma_{p}+\beta\chi(\gamma_{p})}{\lambda\gamma_{p}}\right)^{j}U_{p},\quad
j\geq 1,\quad n\geq 1,$

(a) $j=0$. (b) $j=1$. (c) $j=2$. (d) $j=3$. (e) $j=4$. (f) $j=5$.

Figure 20:  A fixed input $\mathcal{E}_{0}$ (red) which is the superposition
of a tuned curve at $\vartheta=0$ (blue) with some fixed random noise is
presented at layer $j=0$. Profile (yellow) of $\mathcal{E}_{j}^{n}$ at time
iteration $n=200$ along the first layers of the network
$j\in\left\\{1,2,3,4,5\right\\}$.

This is illustrated by a first example simulation in Figures 20 and 21. We
present at $j=0$ a fixed input $\mathcal{E}_{0}$ which is generated as the
superposition of a tuned curve at $\vartheta=0$ (blue) with some fixed random
noise: namely we select $a_{1}=0$, $a_{2}=1$, $a_{3}=0$ and all other
coefficients $a_{p}$ for $p=4,\cdots,d$ are drawn from a normal law with an
amplitude pre-factor of magnitude $\varepsilon$ set to $\varepsilon=0.1$. The
shape of the input $\mathcal{E}_{0}$ is shown in Figure 20(a). The profile of
$\mathcal{E}_{j}^{n}$ at time iteration $n=200$ along the first layers of the
network $j\in\left\\{1,2,3,4,5\right\\}$ is given in Figure
20(b)-(c)-(d)-(e)-(f) respectively. We first observe that the network indeed
acts as a filter since across the layers of the network the solution profile
$\mathcal{E}_{j}^{n}$ has been denoised and get closer to the tuned curve at
$\vartheta=0$. Let us also remark that the filtering is more efficient for
layers away from the boundary and is less efficient for those layers near the
boundary. This is rather natural since the impact of the input
$\mathcal{E}_{0}$ is stronger on the first layers. We see that already at
layer $j=5$, we have almost fully recovered the tuned curve at $\vartheta=0$