MathMods Erasmus Mundus Program
University of L’Aquila
University of Hamburg
University of Barcelona
Master Thesis

Neural network models

_Author:_

Plamen Dimitrov

_Advisors:_
Prof. Leonardo Guidoni Prof. Debora Amadori

###### Abstract

This work presents the current collection of mathematical models related to
neural networks and proposes a new family of such with extended structure and
dynamics in order to attain a selection of cognitive capabilities. It starts
by providing a basic background to the morphology and physiology of the
biological and the foundations and advances of the artificial neural networks.
The first part then continues with a survey of all current mathematical models
and some of their derived properties. In the second part, a new family of
models is formulated, compared with the rest, and developed analytically and
numerically. Finally, important additional aspects and any limitations to deal
with in the future are discussed.

###### Contents

  1. 1 Background
    1. 1.1 Introduction
    2. 1.2 Biological neural networks
    3. 1.3 Artificial neural networks
  2. 2 Mathematical models of neural networks
    1. 2.1 Realistic neuron models
      1. 2.1.1 Conductance-based models
      2. 2.1.2 Multi-compartment models
      3. 2.1.3 Reduced conductance-based models
    2. 2.2 Minimal neuron models
      1. 2.2.1 Spiking (pulsed) models
      2. 2.2.2 Firing-rate (rate-based) models
      3. 2.2.3 Population models
      4. 2.2.4 Binary models
    3. 2.3 Extension to networks
      1. 2.3.1 Activity - study of the potential dynamics
      2. 2.3.2 Connectivity - study of the synapse plasticity
  3. 3 Proposed model
    1. 3.1 Formulation
      1. 3.1.1 Construction of a family of models
      2. 3.1.2 Motivation for the family of models
    2. 3.2 Mathematical analysis
      1. 3.2.1 Universal approximation
      2. 3.2.2 Expanded graph interpretation
      3. 3.2.3 Phase space decomposition
      4. 3.2.4 Network motif study
    3. 3.3 Numerical implementations
      1. 3.3.1 Feature detector
      2. 3.3.2 Object tracker
      3. 3.3.3 Pattern classifier
  4. 4 Discussion
    1. 4.1 Additional considerations
    2. 4.2 Current limitations
    3. 4.3 Future developments
  5. References

##  1 Background

###  1.1 Introduction

The brain is the only human organ that has spawned a wide variety of
disciplines originating in and focused on completely independent approaches
that don’t truly intersect each other for a common theory. Scientific
disciplines like

  * •

psychology

  * •

behavioral science

  * •

cognitive science

  * •

psychiatry

  * •

artificial intelligence

  * •

machine learning

  * •

neurology

  * •

neuroscience

  * •

neuroanatomy

are all devoted to studying different aspects of this organ with the
psychology type sciences taking a higher level (symbolic) approach and the
neuroanatomy type sciences - a lower level (circuit) approach. The hope is
that they will intersect somewhere in the middle, fully explaining higher
level cognitive phenomena through lower level physiological and morphological
such. However, this hasn’t happened as of yet and when it does, it will be
analogical to finally creating a comprehensive coherent theory demystifying
questions about the human mind which are currently only asked in the field of
philosophy.

More engineering oriented fields like artificial intelligence and machine
learning linger in the middle between the two ends with attempts to reproduce
high level phenomena by both more abstract methods (Lisp, A* search) as well
as methods resembling the natural structure of the brain (artificial neural
networks, neocognitron). The ultimate goal of such fields is achieving some
practical use of the implemented algorithms - be it human-dependent tasks like
image/voice/text recognition and classification or self-driving cars and real
time translation. Their approaches usually stem in developing an algorithm for
the job which strives to achieve maximal accuracy rate on specific commonly-
accepted data set. For this purpose, obtaining good empirical results on the
given data set is often sufficient for adoption of the developed techniques
and results in narrow specialization of the implemented algorithms. However,
many results in neuroscience suggest the existence of a universal approach
that could be the potential solution in all different areas of AI application
- from voice and image recognition to language comprehension and translation
[73, 5, 72]. Even though a well-specialized implementation will perform better
on the specific area of application, one must then advocate for modelling and
reproducing more general cognitive behavior and use more mathematical rigor in
order to draw useful conclusions from such generality.

To construct such a universal algorithm however is not an easy task. This is
one of the reasons for multiple disappointments and periods of loss of
interest and high criticism in the area of artificial intelligence that are
commonly referred to as the AI winters. In order to build something, one
should completely understand it, and the lack of complete or at least
sufficient understanding of the brain is among the primary reasons for the
occurrence of so many attempts to explain it.

Figure 1: The complexity of the brain studied from many different scientific
angles gets all participants separated and lost in the end.

Figure 1 is just one of many depictions of the same problem - each discipline
is like a flash of light illuminating different angle of the same giant object
and none of them sees the whole of it or how large it really is. And the
realization about all the limitations that have to be faced is by far not
restricted only to the engineering field. In psychiatry, using psychiatric
drugs has so many side effects because using them to treat complex psychiatric
disorder is a bit like trying to change a car’s engine oil by opening a can
and pouring it all over the engine block - some of it will dribble into the
right place but a lot of it will do more harm then good [3]. This is how
neurological issues related with the brain are solved at the present day after
all the progress in medicine from the last centuries. General developments in
neuroscience like fMRI have been really helpful in pinpointing active areas of
the brain to certain stimuli but the map is very crude and uses increased
blood density rather than actual neuron spikes. Even if we were able to record
the complete history of every spike of every neuron and store the entire
neural circuitry or connectome as proposed by [1], the data generated just
from the brain of the smallest organisms with a minimal diversity of possible
behaviors like a fruit fly would include approximately $10^{4}$ neurons to
analyze.

It is therefore vital to be under no illusion about the difficulties involved
here. Nevertheless, the possibility of existence of such a universal method
alone is attractive enough so that we would like to make a step in this
direction. In particular, we will investigate the mathematical modelling
setting of the scientific fields above, concentrating on neural networks as
the least common denominator among all of them. We will study the reasons and
derivation and include further references for the major neuron models and the
resulting network models in our investigation. Using these as a starting
point, we will then propose a family of models that could be best motivated
with a set of cognitive behaviors they must exhibit. We have compiled these
cognitive phenomena with the expectation that they are general enough to be
shaped into solutions of different areas of application like the ones we
mentioned before. The validity of the models then depends on the validity of
these requirements which is studied in depth both analytically and numerically
through the formulated models. We will discuss it in a final section but in
order to be able to talk about the summary of existing models, we need some
additional background for the relevant life and computational sciences in this
first section.

###  1.2 Biological neural networks

The human brain comprises 100 billion neurons as its principal cellular
elements, each one connected with around 10000 others. Consequently, the
potential complexity of the resulting network is vast in terms of both
possibilities and interpretation. It is called by many the most complex
machine on Earth and perhaps the universe [46, 18]. In addition to the
combinatorial size of possibilities for connectivity, the neurons
differentiate in many different types with different functions and by far
don’t constitute the entire brain. Other cellular types are also present and
even estimated to be five times more than the neurons (90% of the brain),
namely several types of glial cells that play a crucial role in the
maintenance of the neural network, neural development and in the generation of
myelin which is used for isolation and is one of the main reasons for fast
impulse conductance. Despite their essential supportive nature, we will not
discuss them in more detail and instead focus on the main communication
infrastructure. This section will begin with an overview of the morphology and
physiology of the neuron, the main elementary node or element of the brain,
then provide a short description of the ways in which individual neurons
communicate with each other in pairs, networks, and systems, and finish with
plasticity as one of the most characteristic neural network properties.

The 100 billion neurons in the brain share a number of common features. The
anatomical variation of these neurons is large, but the general morphology and
their electrical and ligand dependant responsiveness allows these cells to be
classified as _neurons_ (coined in 1891 by Wilhelm von Waldeyer) [56]. Neurons
are different from most other cells in the fact that they are polarized and
have distinct morphological regions, each with specific functions.

Figure 2: Neuron morphology and physiology.

Neurons are generically characterized by a central cell body or _soma_ that
comes in different shapes. The soma contains the cell nucleus and most of the
genomic expression and synthetic machinery producing the proteins, lipids, and
sugars that constitute neuron’s cytoplasm and membranes. For a general neuron,
the soma extends into inputs and outputs as in figure 2.

The region where a neuron receives connections from other neurons or the input
pole consists of extensively branching tree-like extensions of the soma
membrane known as _dendrites_ (coined in 1889 by William His from dendros
meaning ”tree” in Greek). The dendrites arise directly from the cell body in
vertebrate neurons similarly the ones in the figure. However they arise from
the axon in invertebrate neurons. Furthermore, the body is also a receiving
site in most neurons.

The output pole, called the _axon_ (coined in 1896 by Rudolph Albert von
Kolliker) arises as a single structure from the soma (and occasionally from a
dendrite). The axon conducts propagating electrochemical signals termed
_action potentials_ (also spikes or impulses) that are usually initiated at
the axon base or _hillock_ and move away from the soma to the _terminal_
regions of the neuron. Axons can be rather long extending up to about a meter
in some human sensory and motor nerve cells. A _synapse_ in the terminal
region of the axon is the place where one neuron forms a connection with
another (called _postsynaptic_ neuron, the first one being respectively
_presynaptic_) and conveys information through the process of synaptic
transmission.

However, it is important to note that there are many exceptions to these
general rules of neuronal organization. There are some dendrites that also
serve as output systems [69]. In some neurons, e.g. peripheral sensory
neurons, the input occurs via axons. Neurons can have many types of branching
or no branches at all, examples being receptors cells in the carotid glomus,
in gustatory system in vertebrate tongue or as photoreceptors in the retina
[56].

Last important morphological feature is that the presynaptic cell is not
directly connected to the postsynaptic cell - the two are rather separated by
a gap known as the _synaptic cleft_. Therefore, the communication between a
presynaptic and a postsynaptic neuron (_synaptic transmission_) is realized
through a _neurotransmitter_ , i.e. a chemical messenger released through an
action potential at the presynaptic terminal (a process called _exocytosis_)
that diffuses at the synaptic cleft and binds to selective _receptors_. The
binding to the receptors leads to a change in the permeability of ion channels
in the membrane and in turn a change in the membrane potential of the
postsynaptic neuron known as a _postsynaptic synaptic potential (PSP)_ which
may then lead to an action potential in the postsynaptic neuron. This brings
to the conclusion that signaling among neurons is associated with changes in
their electrophysiological properties.

Now let us provide more details and important terminology regarding the
eletrophysiological properties of a neuron. The change in the PSP caused by a
presynaptic action potential may be a decrease in the polarized state of the
membrane called _depolarization_ or an increase called _hyperpolarization_.
The membrane potential in the absence of any presynaptic stimulation is about
-60 mV inside with respect to the outside of the postsynaptic cell and is
called _the resting potential_. Hyperpolarization then makes the potential
inside the cell even more negative while sufficiently large depolarization is
the one to trigger an action potential. The action potential is associated
with a very rapid depolarization to achieve a peak value of about +40 mV in
the brief period of 0.5 msec [6]. The peak is followed by an equally rapid
period of _refractoriness_ when the neuron cannot be excited which associated
also with a _repolarization_ phase that completes a depolarization-
repolarization cycle.

The voltage at which the depolarization becomes sufficient to trigger an
action potential is called a _threshold_. It can be reached through multiple
presynaptic spikes leading to smaller changes in the membrane potential which
accumulate in time (_temporal summation_ which together with spatial summation
over the dendritic area is also called _synaptic integration_) or through a
single suprathreshold spike but the resulting action potential is identical in
amplitude, shape, and duration. The action potential will also not change if
the input stimulus leads to a depolarization much larger than the threshold.
The presynaptic spikes can also cause hyperpolarization which reduces the
change of action potential on the postsynaptic side. The PSP can therefore be
divided in to _excitatory postsynaptic potential (EPSP)_ which is the one
responsible for the depolarization and _inhibitory postsynaptic potential
(IPSP)_ which is the one responsible for the hyperpolarization. An IPSP is
called inhibitory because it tends to prevent the postsynaptic neuron from
firing an action potential thus regulating the ability of excitatory signal to
bring about this same event. The conclusion from this is that synaptic
transmission has two basic forms, namely _excitation_ and _inhibition_.

Generally, a single action potential in a presynaptic cell does not produce an
EPSP large enough to reach threshold and necessarily trigger an action
potential. But longer-duration of even subthreshold stimulus can lead to
multiple action potentials (spike sequences or _spike trains_) with frequency
depending on the intensity of the stimulus. The same is true for receptor
cells from pressure of touch to intensity of light as well as motor cells
where the output frequency determines the contraction level of muscle [56].
This, together with the identical properties of action potentials suggests
that the nervous system encodes information in frequency of isotropic pulses
instead of their amplitude and shape.

Figure 3: Ion gates comprise ion pumps that maintain nontrivial equilibrium
potential of the membrane and ion channels which are responsible for the
active electric properties of a neuron [63].

The generation of action potentials is an example of the active electrical
properties of neurons because they are brought about by active, voltage-
dependent means [56]. More specifically, the membrane potential might be
affected by activation of ion channels which is caused by the neurotransmitter
release on the presynaptic side. The passive electrical properties of neurons
are also very important for describing their overall dynamics in mathematical
models presented in later sections. The cell membrane as a bilipid layer is a
nearly perfect insulator but has _ion gates_ as proteins embedded into it.
These ion gates can be _ion pumps_ or _ion channels_. There are also _non-
gated ’leakage’_ channels (constantly open channels permeable to potassium)
which we mention for completeness but later on we will only consider the ion
gates.

Ion pumps permanently (for an energy cost) transport ions from one side of the
membrane to the other (depending on the ions transported and therefore the ion
pumps transporting them, e.g. sodium is less while potassium is more
concentrated inside) maintaining a voltage difference through a difference in
the ion concentrations on the inside and on the outside of the neuron. An ion
channel will then reach an equilibrium when there is no more flow of ions
between the two sides. The membrane potential at equilibrium is the
_reversal/Nernst potential_ of that particular ion assuming a single ion
dominated system. Among the main ion species found in the central nervous
system are sodium, potassium, calcium, and chloride ions.

Similarly the neuron’s morphology, there are many exceptions when it comes to
its physiological properties. Although the conductance of most ion (voltage-
gated) channels is increased by membrane depolarization, the conductance of
some channels is increased when the membrane is hyperpolarized. In addition to
axonal spikes, action potentials can also travel along dendrites either in the
direction of the soma (centripetal spikes) or away of the soma and reach the
most distal dendritic branches (centrifugal spike) depending mostly on the
distribution of ion channels over the dendritic tree. Axons can also branch
out into parallel pathways early in the spike propagation or later on into a
distal dendrite, in both cases sending very similar spike trains but with
different resulting conduction due to different axonal diameters and possible
spike failure.

Now let’s look at the network level and more specifically at the general types
of network connectivity that is experimentally observable. There are three
major classes of connections in all areas of the brain - _feedforward_
(bottom-up) connections transfer activation to neurons which are further along
the processing path that started with external input, _lateral_ (recurrent)
connections are used for communication among neurons that are considered to be
within the same stage along the processing path, and _feedback_ (top-down)
connections which project activation back. The feedforward and feedback
connections are typically comparable in terms of numbers while the lateral
connections outnumber each. The _feedforward networks_ are built only from
feedforward connections and exclude the possibility of cycles, the _recurrent
networks_ include lateral connections as well, and the _fully recurrent
networks_ allow for connections in any direction among all neurons. This is
also the reason why when talking about fully recurrent networks, one talks
simply about connections and and forgets the overall categorization.

To complete the various anatomical organizations in the brain, we must include
a description of the system level structure. The emphasis of any network model
described later on is mainly on _neocortical circuits_ so we will briefly
outline some principal organizations of the _neocortex_ , the outer and
developed later in evolution layer of the _cerebral cortex_. Even though the
neocortex or simply _cortex_ can be horizontally split into four lobes
responsible for different cognitive functions and modalities, its detailed
structure is similar among all of them which distinguishes it significantly
from older parts of the brain like the _brainstem_ [83]. The structure
observed in the neocortex comprises of six (to ten depending on enumeration)
_layers_ where the fourth layer is generally viewed as an input layer (due to
predominant number of white matter afferents) while the fifth layer projects
mostly outwards and is therefore viewed as an output layer. The _white matter_
consists mostly of axons and connects distant areas which implies that it must
be used for global communication. The second and third layers are thought to
be responsible for long range lateral (within layer) communication. The sixth
layer neurons mostly project into the first layer. The neurons are also
arranged in _columns_ which are thought of to be the functional units (with
neurons within the column storing redundant information) since similar
connectivity patterns are observed within and among these columns.

Finally, perhaps the most significant property of neural networks which makes
them an evolutionary advantage in the long term is their _synaptic
plasticity_. A persistent (ten minutes or longer) increase in synaptic
transmission efficacy is called _long-term potentiation (LTP)_ and a decrease
of the same kind is called _long-term depression (LTD)_. There could be
multiple physiological and biophysical mechanisms for realizing such synaptic
plasticity, among them changes of the number of release sites, the probability
of neurotransmitter release, the number of transmitter receptors, and the
conductance and kinetics of ion channels all of which have been demonstrated
for specific cell types in multiple papers [83]. Additionally, a spike timing
as well as calcium dependent plasticity was observed empirically in some cells
and is examined in further detail in some of the models and learning rules to
follow.

###  1.3 Artificial neural networks

No doubt a large portion of the models developed concerning neural networks
are algorithms without a rigorous mathematical treatment but achieving a very
good measurable accuracy at a specific task. These are computational models
within the field of machine learning that deserve at least a brief overview
here together with any computational terminology involved in the following
sections. With the increasing advanced and complexity of such algorithms or
_artificial neural networks (ANNs)_ , we should also make an effort to
understand them a bit more rigorously and not just their biological
counterparts. We will start with an overview of the basic terminology and
concepts used in machine learning and therefore important also for the ANNs
and our numerical implementations later on. We will then move to the classical
ideas and developments and into their latest extensions with deep learning. We
will complete the section with some existing criticisms about all these
computational models.

All machine learning algorithms are based on induction as opposed to (e.g.
mathematical) deduction. As the name of the discipline implies, they are
learning to perform a specific task from empirical data which cannot be
algorithmically specified in any straightforward way. If the learning they
perform is _supervised_ , the task is to predict an output value or _label_
from an input value given a _training set_ of examples as input-output value
pairs. It is called supervised due to the presence of the training labels that
teach the network the right answers as opposed to _unsupervised learning_
where these are not provided. There are also paradigms in between like
_reinforcement learning_ where the right answers are not provided but the
algorithm now termed _agent_ still receives feedback in the form of a _reward
or punishment_ from the environment as well as _semi-supervised learning_
where both supervised and unsupervised approaches are combined during the
_training phase_ with both _labelled_ and _unlabelled data_.

The two major tasks in supervised learning are _classification_ where the
algorithm must predict a discrete label value and _regression_ where the value
is continuous. As such, they define supervised learning as a task of
approximating a function over all possible inputs from finitely many known
values. While this is also studied in much depth in statistics, the point here
is to find algorithms that could gradually learn (approximate and use to
predict labels for new inputs) the function that generalizes best over all
data. For this purpose, a _test set_ is often used in addition to the training
set in order to avoid _overfitting_ , i.e. performing too well on the training
set and a lot worse on a ’fresh new’ test set). _Underfitting_ is another
although generally smaller danger. Overfitting is often caused by too many
free parameters or too small training set so it could also be reduced by a
technique called _cross-validation_ which helps in selecting these parameters
by partitioning the training set and obtaining validation errors through
testing on each partition and training on the rest. Comparison in performance
of the algorithms in terms of _accuracy rate_ or _test error_ is then possible
on standardized _datasets_ that are shared and known in the machine learning
community. Major tasks in unsupervised learning are _clustering_ (while there
are no explicit labels for the inputs, similarities among them are still
inferable), _dimensionality reduction_ (removing unnecessary information for
easier interpretation or in the case of ANNs - compression performed by
_autoencoders_), and _outlier detection_.

It is surprising how many practical applications (from navigating a road to
speech recognition) can be translated into the tasks described above and that
all of the tasks can also be performed well by ANNs among many other machine
learning algorithms. This is at least the case after multiple periods of
abandonment due to their realized limitations and earlier overestimation of
their capabilities typical for the artificial intelligence field.

The very first major work in the direction of pure computation with regard to
the nervous system was done by a neurophysiologist Warren McCulloch and a
mathematician Walter Pitts who showed that a simple model of a neuron (simply
a threshold sum of inputs described also in later sections) can reproduce the
basic AND/OR/NOT logical functions [62]. It was very influential especially
because of the belief that logical reasoning could solve AI [50] but lacked
the crucial ability to learn so was later on extended by a neurobiologist
Frank Rosenblatt in his _perceptron_ model to include weights (and a bias term
similar to a current injected straight to the neuron’s soma) that could make
one input more influential than another [75].

The guiding principle for learning that he used was based on another
fundamental work by the psychologist Donald Hebb who conjectured that synaptic
changes are driven by correlated activity of pre- and postsynaptic neurons,
i.e. neurons that fire together also wire together [34]. More formally, the
strength of a synapse from neuron A to neuron B will increase if the firing of
neuron A often contributes to the firing of neuron B through that synapse.
_Hebbian learning_ is the way this conjecture is usually referred to and can
be generalized in various ways, one in particular to also reduce the strength
of the synapse should the neuron A repeatedly fail to contribute to the firing
of neuron B.

The perceptron then learns by readjusting the weights in the direction of
lower error if its output does not match the desired label and thus correlates
better the desired and actual outputs. Due to the presence of sum thresholding
or _activation function_ for the sum of inputs, the perceptron performs binary
classification which can be extended to multi-category classification by
adding a _unit_ or neuron for each category to form a _layer_. This layer is
an _output layer_ fully connected to the inputs and the resulting two layer
network is one of the most canonical examples of an ANN - a perceptron ANN
approximating a vector-valued function. One way to classify the vector of
input states is thus to choose the maximum weighted sum among the outputs.

An even more daring deviation from biological plausibility is Widrow and
Hoff’s _ADALINE_ [87]. Even though its main contribution is in implementing
neurons in electrical circuits through resistors with memory (memistors), it
explores the option of simplifying the perceptron altogether by dropping the
activation function and simply outputting the weighted sum of each neuron.
This linear form is easier to study mathematically because it removes jump
discontinuities and nonlinearities from the model neuron. In addition, the
learning/approximation problem can be restated as optimization problem where
the derivative of the error can be used in a stochastic or other gradient
descent.

Although _connectionism_ as the idea that complex computations can be
performed by many simpler connected units comes with these first computational
models, the models themselves are fundamentally limited in their computational
power. A major critic of the perceptron for instance is that it cannot
reproduce the XOR logical function since it violates the linear separability
condition [66]. A solution to this, extension to the perceptron, and an
eventual exit of the first AI winter is the addition of _hidden layers_ and
_backpropagation_ learning. If there are multiple layers between the input and
output layer responsible for multiple levels of abstraction and feature
extraction, much more complex computations including approximating the XOR
functions can be performed. Since the first weight optimization approach
described above can no longer be applied (we only know the error at the output
layer), backpropagation of the error from the output layer to the weights of
all hidden layers using chain rule splits and redistributes the correction in
the entire network. This approach together with the greater computational
power of the multilayer ANNs was shown in a sequence of papers most notable
from which is the description by Rumelhart, Hinton, and Williams [76]. A
result of major mathematical importance to once and for all solve any
approximation issues is the proof that multilayer ANNs are _universal
approximators_ , i.e. with sufficiently many hidden units in a middle layer
they can approximate any Borel-measurable function between finite dimensional
spaces to an arbitrary degree of accuracy [42].

The multilayer ANNs with backpropagation learning are the state of the art
ANNs until the present which after a second AI winter returned with a few
small modifications as _deep neural networks or DNNs_. These modifications
involve more careful weight initialization and choice of activation function
[28]. The paper which is among the first to talk about ”Deep learning” is
again due to Hinton [36] and essentially considers a variant of DNN, namely
’deep belief network’ described below, which is optimally trained layer-by-
layer through a semi-supervised learning (unsupervised learning to initialize
the weights). Such variations of ANNs/DNNs are what remains to be discussed
until the end of this section. The major ANN classification is based on the
previous connectivity classification of feedforward, recurrent and fully
recurrent networks. The DNN variations are then based on modality
specialization (CNNs for visual input, LSTMs for sequential input) and
connectivity (DBNs as fully recurrent networks).

As a variation of an ANN/DNN which is very successful in computer vision, a
_convolutional neural network or CNN_ was first applied to handwritten zip
code recognition by LeCun [52]. The main practical improvement that a CNN
introduced involves the addition of a _convolution layer_ where instead of
each neuron having a unique weight for each of its inputs, we have neurons
with weights that are identical (and identically modified) at multiple
disjoint subregions of the same size within the previous layer. In this way we
don’t have to learn the same feature like specific edge orientation by
multiple neurons and the neuron’s selectivity becomes invariant with respect
to the position of the feature in the image. This draws inspiration from
Fukushima’s neocognitron [24] and in return from Hubel and Wiesel’s emprical
results on the hierarchical organization of the visual nervous system [45]
since the complex cell observed and modelled there pertain a certain degree of
spatial invariance. A second practical improvement that resembles also the
hierarchical arrangement of edge-sensitive simple, complex to hypercomplex
cells there is the addition of a _pooling layer_ as a form of down-sampling so
that a layer is projected in a smaller size successive layer. In this way the
number of free parameters and therefore chance of overfitting as well as the
training time is reduced while a neuron from a later layer ’sees’ the features
from the previous layer and thus a larger region in the original input image.

Another somewhat successful variation of an ANN/DNN is the _deep belief
network or DBN_ which as their name speaks is capable of internal
representations of the perceived inputs. The groundwork for a DBN is a
_Boltzmann machine_ which is a _graphical (graph) model_ performing _energy-
based learning_ with nodes as hidden and visible stochastic variables where
the visible variables determine the probabilities of the hidden variables. The
essential symmetry feature of Bolzmann machines is that the probability of a
visible variable (input) can be calculated from a hidden variable (output)
which allows for the original input to be generated from an internal
representation. This models are called _generative graphical models_. DBNs are
then multilayer networks of restricted Boltzmann machines (only with
feedforward/feedback connections) which were introduced and later on sped up
in performance by Hinton, Neal, and Dayan. In successive papers, they included
recognition weights (feedforward) and generative weights (feedback) and a
’wake-sleep’ algorithm for unsupervised training [35].

A variation of a recurrent neural network, which excels in processing
sequential information such as text or speech recognition, is the _Long Short
Term Memory or LSTM_. It has significant advantages over its preceding
competitors like time-delay neural networks (separate weights for multiple
discrete times within a sliding window) and recurrent networks with
backpropogation through time in the fact that the latter two suffer from
inability to learn long-term information - the first due to the fixed window
length and the second due to limitations of backpropagation for deep learning
[50]. The LSTM solution by [38] contains special types of units called
Constant Error Carousels (CECs) that have an activation function, an identity
function, and a connection to itself with weight of $1.0$. In this way errors
as derivatives of identity (chain rule) don’t vanish in the long term and the
network also discovers and learns correlations that are very distant in time.
CECs are then connected to several nonlinear adaptive units with possibly
multiplicative activation functions for learning nonlinear behavior.

A recurring and often criticized issue in all of these models is their
biological plausibility [15]. Since they are purely computational, they are
optimized for accuracy on comparable datasets like MNIST, Imagenet, TIMIT, to
name a few. Since they don’t necessarily model natural phenomena and only draw
inspiration from such, they get the engineering freedom of deviating away as
much as they need as long as they achieve their original goal. At the same
time biological resemblance is what distinguishes neural networks from other
machine learning techniques like random forests and support vector machines
which also perform with very high accuracy ratios but don’t assume such
resemblance. In this sense, ANNs are a middle ground between features taken
from nature and features taken from engineering ingenuity. The problem of
balancing the two makes a large difference in the final performance of each
computation model. It is therefore possible that engineering constructions
like

  * •

the backpropagation learning of an ANN

  * •

the identical weights on a region of a fixed size of a CNN

  * •

the probabilistic (although not necessarily) feedback connections of a DBN

  * •

the CECs of a LSTM

could all have more natural and universal explanations. If a computational
deviation is necessary, it would still be preferable to try to find the most
natural solution possible and explicitly state any unavoidable computational
reductions.

Other sources of significant criticism concern the lack of feedback
connections in most of these models as well as the phenomenon of _catastrophic
interference/forgetting_ whereby a completely new input pattern may lead to
very fast forgetting of the already learned one. Considering feedback
connections is probably important at least since such connections are
empirically observed to be approximately the same in number as the feedforward
connections. However, most of the ANN models ignore them which implies
ignoring a significant part of the overall BNN (biological neural network)
dynamics. A group of fully recurrent networks that take a very different
direction from DNNs and which we only mention here for completeness are self-
organizing or Kohonen maps (SOMs) [49] as well as models from the Active
Resonance Theory (ART) [8]. The ART proposal also solves the catastrophic
interference problem but is statistically inconsistent and uses centralized
(instead of parallel and purely local) policies. By local policies we mean
policies that make use of entirely local information within a synapse or
neuron and not global information about the network state of any kind. A more
natural solution would therefore resemble swarm intelligence (SI) where
decentralized decisions lead to complex self-organized systems.

##  2 Mathematical models of neural networks

###  2.1 Realistic neuron models

In our review of neuron models we will mainly include their mathematical
formulation and sometimes their construction. For any additional mathematical
study we will include the main references where this study was performed. Some
predictions and matched observations in the cognitive sciences are also
mentioned with further directions to their main sources for a deeper
investigation.

One of the most important tasks in modelling neurons together with any
resulting neural networks is to capture neuron’s essential features for a
desired cognitive level behavior while choosing among many different levels of
biological plausibility. While the belief that the more we pertain to the
complex empirical reality the greater the chance that we will reach our goal
is not entirely false, in most cases this will fail in practice due to the
computational resources for and multiple possible interpretations of the
complexity that we have added. However, essential features of the neuron might
still reside in the more realistic models even though there is already a wide
variety of reduced to computationally minimized models that have made their
selections of what to retain and what to lose. For this reason, it is
important that we include these neuron models as well.

####  2.1.1 Conductance-based models

The most biologically detailed yet still sufficiently general neuron models
are the ones including dynamics arising from the currents that pass through
the ion channels in the cell membrane. The membrane conductance then results
in action potentials and becomes the main focus of this type of models. The
Hodgkin-Huxley neuron introduced in [39] is a milestone in computational
neuroscience and a starting point for many easy extensions and therefore
cannot be skipped in any study of conductance-based models of a neuron.

Figure 4: Electrical circuit analogue to the Hodgkin-Huxley neuron model as
presented in [39].

While studying the giant axon of the squid, Hodgkin and Huxley found three
types of current: sodium, potassium, and a leak current of Cl- ions which also
takes care of other channels. Since the membrane lipid bilayer isolates the
ions inside and outside the cell body it plays the role of a capacitor. Since
the ion pumps lead to difference in the concentration of ions and therefore a
nonzero equilibrium potential (Nernst potential), they play the role of a
battery. The three currents from the ion channels and the current charging the
capacitor comprise the total applied current.

$I(t)=I_{C}(t)+\sum_{k=1}^{3}I_{k}(t)=I_{C}(t)+I_{Na}(t)+I_{K}(t)+I_{L}(t)$
(1)

We can easily construct an equation for a membrane conductance of a passive
neuron, i.e. including only the $I_{C}$ and $I_{L}$ terms using the above
analogies with electrical circuits. From the definition of current $I=dQ/dt$
and the definition of capacity with $Q=VC$ where $Q$ is the charge, $V$ is the
voltage, and $C$ is the capacitance, we get that $I_{C}=dQ/dt=CdV/dt$. From
Ohm’s law $V-E=\Delta V=IR$ where $E$ is the resting potential (steady state
voltage for a passive membrane) and $R$ is the resistance of the leak channel,
we get that $I_{L}=\Delta V/R=(V-E)/R$. The passive membrane equation would
then follow from the conservation of current $I_{C}+I_{L}=0$. For each of the
ion channels we have a similar equation to the leak channel and can use
conductance as the reciprocal of resistance $g_{L}=1/R$ to simplify notation
with the only difference that the leak channel conductance is a constant while
the the other channels can open and close implying conductance that depends on
time and voltage. This means that we can simply use their maximum conductances
$g_{Na}$ and $g_{K}$ multiplied by a number between 0 and 1 which will
represent the probability of a ion channel to be open. We can then expand the
currents and rearrange (1) like

$C\frac{dV}{dt}=I_{C}(t)=I(t)-(g_{Na}m^{3}h(V-E_{Na})+g_{K}n^{4}(V-E_{K})+g_{L}(V-E_{L}))$

where the sodium channels have two coefficients (gating variables) $m$ and $h$
since they have activation and inactivation gates and all such coefficients
are raised to powers to match experimental data. $E_{Na}$ and $E_{K}$ are now
the reversal potentials and the quantities $V-E_{K}$ are the driving forces
for each ion. To complete the model, we add differential equations to govern
the three gating variables which are based on protein kinetics with activation
rate $\alpha$ and inactivation rate $\beta$:

$\displaystyle\begin{dcases}&C\frac{dV}{dt}=I(t)-(g_{Na}m^{3}h(V-E_{Na})+g_{K}n^{4}(V-E_{K})+g_{L}(V-E_{L}))\\\
&\frac{dm}{dt}=\alpha_{m}(V)(1-m)-\beta_{m}m\\\
&\frac{dn}{dt}=\alpha_{n}(V)(1-n)-\beta_{n}n\\\
&\frac{dh}{dt}=\alpha_{h}(V)(1-h)-\beta_{h}h\end{dcases}$ (2)

This model by Hodgkin and Huxley captures the essence of spike generation by
sodium and potassium ion channels but it could easily be extended to other
types of ion channels. All a protagonist of biological realism has to do is
add an equation for the dynamic of the new ion channel and a term in the main
equation sum. A general form system can then be

$\displaystyle\begin{dcases}&C\frac{dV}{dt}=I(t)-\sum_{k=1}I_{k}(t)=I(t)-\sum_{k=1}g_{k}m_{k}^{p_{k}}h_{k}^{q_{k}}(V-E_{k})\\\
&\frac{dm_{k}}{dt}=\alpha_{m,k}(V)(1-m_{k})-\beta_{m,k}m_{k}\\\
&\frac{dh_{k}}{dt}=\alpha_{h,k}(V)(1-h_{k})-\beta_{h,k}h_{k}\end{dcases}$ (3)

The ion channels that have been added in various models include:

  1. 1.

Sodium channels

    1. (a)

Fast sodium current $I_{Na}$ is similar to the one in the Hodgkin-Huxley
model.

    2. (b)

Persistent sodium current $I_{NaP}$ lack the inactivation gate and therefore
the variable $h=1$.

  2. 2.

Potassium channels with current $I_{K}$

    1. (a)

Rapidly inactivating potassium current with constant 1 $I_{A1}$

    2. (b)

Rapidly inactivating potassium current with constant 2 $I_{A2}$

    3. (c)

Slowly inactivating potassium current with constant 1 $I_{K2a}$

    4. (d)

Slowly inactivating potassium current with constant 2 $I_{K2b}$

  3. 3.

Calcium channels with current $I_{Ca^{2+}}$

    1. (a)

Low threshold calcium current $I_{T}$ \- no inactivation gate

    2. (b)

High threshold calcium current $I_{L}$ \- inactivating current where the
channel shuts down after depolarization

Different combinations of these give rise to different dynamics which in some
cases may only be quantitatively different but in others also qualitatively
[27]. A good example for this is the case of low threshold calcium current
which predicts and explains a phenomenon of post-inhibitory rebound, i.e. an
action potential arising from an inhibitory input. Another model, this time
with the addition of rapidly inactivating potassium current, is known as the
Connor-Stevens neuron model [14]. disable,inlinedisable,inlinetodo:
disable,inlinePDE version: derivation, bistability, periodicity - no
inhibition

####  2.1.2 Multi-compartment models

A second direction towards greater biological plausibility of a neuron model
one could take is to add spatial variance within the cell. Complex neuron
geometries and natural differences in the membrane potential across different
compartments can bring about different dynamical behavior which can then be
accounted for by such models. Such differences in current are especially
noticeable along long and thin dendrites. Considering the above models as
single-compartment models, the step to make must now be clear - simply
separate the neuron into different regions, each governed by the same single-
compartment model but with different ion channel concentrations (potential and
gating variables) and add longitudinal resistance among the compartments in
addition to the transversal resistance of each compartment’s ion channels. Of
course, this assumes that the variation in potential within a compartment is
either negligible or produces a small enough error that can be ignored for
interpretation of the results.

An example for a non-branching cable would be

$C\frac{dV_{m}}{dt}=I_{m}(t)-\sum_{k=1}I_{k}^{m}+g_{m,m+1}(V_{m+1}-V_{m})-g_{m,m-1}(V_{m-1}-V_{m})$
(4)

where $I_{m}$ is the current flowing into compartment $m$ and the last two
terms are responsible for the interaction with the previous and respectively
next compartment (with one of the two terms missing in the case of the cable
ends). The constant $g_{m,m^{\prime}}$ determines the conductance (and
therefore resistance) from compartment $m$ to compartment $m^{\prime}$. The
sum of currents can be replaced with the ion currents similarly to the way we
did it in (1) but with different channel conductance for each compartment.
This can then be extended to various branching-based coupling by adding
removing terms similar to the last two terms [17].

Multi-comparment models can also be used with simpler single-compartment
models like the LIF neuron introduced in later sections which could leverage
some of the computational challenges they pose when increasing the
comparmental resolution to simulate a complex dendritic structure. A further
major disadvantage of such multi-compartment models is the large number of
free parameters where multiple combinations of very different values can
produce the same plausible biological behavior [81].

The best approximation of the potential dynamics on a long and narrow
dendritic cable can be reached by making a step beyond such multi-compartment
models and into cable theory using a nonlinear diffusion equation called the
cable equation [70]. The best approximation of the potential dynamics along an
entire dendritic tree can then be achieved with multi-compartment models that
use the cable equation instead of a number of ODEs (simpler approach above),
thus accounting for spatial variations within a single-compartment model as
well. Similar models are sometimes referred to as Rall-type models and some of
them made interesting and true predictions like the dendro-dendritic spikes
[71] but they are beyond the scope of our review.

####  2.1.3 Reduced conductance-based models

A major advantage of the very biologically detailed models proposed above is
the easy interpretation of all quantities and comparability with experimental
data. However, a major drawback for the analysis of networks of such neurons
is the computational challenges they pose for simulation as well as the
complexity of their interpretation on a network-wide scale. One major attempt
to amend this is the reduction of the dimensionality of the dynamical system
representing a full conductance-based model into usually 2-dimensional such
with one variable that could be loosely related with the membrane potential
and one variable which is usually called a ”relaxation” variable.

The first and currently classical such reduction is the Fitzhugh-Nagumo neuron

$\displaystyle\begin{dcases}&\dot{u}=c(v+u-u^{3}/3+z)\\\
&\dot{v}=-(u-a+bv)/c\\\ &1-2b/3<a<1,0<b<1,b<c^{2}\end{dcases}$ (5)

which can be obtained from the damped oscillator

$\ddot{u}+k\dot{u}+u=0$

by replacing the damping constant with a quadratic term of $u$ to obtain a Van
der Pol oscillator

$\ddot{u}+c(u^{2}-1)\dot{u}+u=0$

and then performing a Lienard transformation [22]. The two unknowns are $u$
for membrane potential (or simply voltage) and $v$ for a recovery variable
(also called accommodation or refractoriness), $a$, $b$, and $c$ are constant
coefficients within certain bounds and $z$ represents an external input
stimulus.

Phase plane analysis shows that this model retains the qualitative features of
the Hodking-Huxley model (2) namely threshold-spike dynamics, rafractoriness
period, spike accommodation and other physiological phenomena observed at the
cardiac muscle excitation. When perturbed sufficiently far from its
equilibrium point at the resting potential, the neuron will essentially
perform an excursion through the phase space and return to this equilibrium
point. This trajectory then represents the action potential of the neuron. In
addition to its retained properties in lower dimensionality, (5) also shows
that (2) is a particular member of a large class of nonlinear systems showing
excitable and oscillatory behavior.

Further reductions build on top of (5) in various ways. Most widely used as
such is the Hindmarsch-Rose neuron

$\displaystyle\begin{dcases}&\dot{u}=v-F(u)+z-w\\\ &\dot{v}=G(u)-v\\\
&\tau\dot{w}=H(u)-w\end{dcases}$ (6)

where $z$ retains the meaning if input stimulus, $F(u)$ is a cubic and $G(u)$
is a quadratic polynomial of $u$ both generalized from the Fitzhugh-Nagumo
equations but with addition of a third unknown $w$ asymptotically converging
towards a linear function $H(u)$ [74]. The model (6) could reproduce all
complex neuron behavior presuming we have found correct functions for $F$,
$G$, and $H$ [48].

Yet another model that has become popular due its exact and measurable
parameters is the Morris-Lecar neuron [67] and a recent one that represents
very well all membrane potential dynamics while remaining very computationally
efficient is the Izhikevich neuron namely

$\displaystyle\begin{dcases}&\dot{u}=0.04u^{2}+5u+140-v+z\\\
&\dot{v}=a(bu-v)\end{dcases}$ (7)

with the auxiliary after-spike resetting

$\displaystyle\text{if }u\geq 30\text{mV then
}\left\\{\begin{array}[]{ll}u\leftarrow c\\\ v\leftarrow
v+d\end{array}\right.$

where $u$ and $v$ are the membrane potential and recovery variable as usual
[47, 48]. According to calculations in [48], all reduced conductance-based
(also called Hodgkin-Huxley type) models described here are performing at
least an order of magnitude better than the original Hodgkin-Huxley model (2)
when it comes to simulating a time window of 1ms (calculated in floating point
operations per second or FLOPS). At the same time most of them are capable of
reproducing the larger part of the most prominent features of biological
spiking neurons. The least amount of FLOPS for the simulation is held by the
minimal neuron models that we review in the next part.

###  2.2 Minimal neuron models

The computationally most effective neuron models usually retain only threshold
and sometimes subthreshold dynamics of the membrane potential. They do this in
an abstract fashion excluding entirely the electrophysiological dynamics that
cause the (sub)threshold behavior on the first place. Due to the simple form
of their equations, we usually immediately decompose the soma input into a sum
of multiple presynaptic contributions and specify the nature of these
contributions as well as the neuron output.

There are various assumptions about the nature of the communicated signal or
neural code [17]. Neurons communicate via action potentials as nearly
identical pulses, thus the most realistic variant of communicating neurons and
their resulting networks would be to preserve this feature. If no significant
information is encoded into the particular times and frequencies of such spike
trains, we could simplify them by taking averages in time and deriving spike-
count rates. The resulting neuron models are called firing-rate models. If no
significant information is encoded in a local population of neurons, i.e. they
have a nearly identical response and similar connectivity, we can take
averages within a population and deriving population rates. The resulting
neuron models are called population models. Clearly, the absence of time
coding (mutually independent spikes) and/or population coding (mutually
independent firing probabilities for individual neurons) is a strong
assumption with some experiments that already contradict them [17].

####  2.2.1 Spiking (pulsed) models

Spiking neuron models are believed to comprise the third generation of neural
networks after the rate-based or continuous second generation (reviewed next)
and the binary first generation (reviewed last) [58]. Despite these bold
claims they have had some difficulty in finding their predicted application in
reproducing cognitive phenomena on a sufficiently high level, mostly due to
computational restrictions as well.

One of the major neuron models which are minimized in terms of the biophysical
mechanisms but bring about action potentials is the leaky integrate and fire
model first introduced by [51]. It is very limited in terms of physiological
complexity but reproduces both the threshold phenomena and the subthreshold
dynamics of a neuron fairly well.

The derivation is similar to and simpler than the derivation of (2). From the
conservation of charge we have this time for $k=1$

$I(t)=I_{C}(t)+\sum_{k=1}^{1}I_{k}(t)=I_{C}(t)+I_{L}(t)=C\frac{dV}{dt}+g_{L}(V-E_{L})$

and taking simpler resting potential $E_{L}=0$ together with $g_{L}=1/R$

$I(t)=C\frac{dV}{dt}+\frac{(V-0)}{R}=C\frac{dV}{dt}+\frac{V}{R}$

Finally, multiplying with $R$ and replacing $\tau=RC$ as the time scale
constant of the leaky integrator

$RI(t)=RC\frac{dV}{dt}+V=\tau\frac{dV}{dt}+V$

The final leaky integrate and fire neuron then reads

$\tau\frac{dV}{dt}=-V(t)+RI(t)$ (8)

with the additional condition that when V reaches a threshold value
$V_{threshold}$ it will be reset back to a value $V_{reset}$ to take care of
the threshold dynamic of the neuron. After resetting, the model is governed by
the same ODE until another discrete fire/reset event occurs. The combination
of the fire condition and the leaky integration (decaying to $RI$ steady
state) defines the leaky integrate and fire neuron model.

One major problem with the LIF neuron is oversimplification and the lack of
reproducibility of rather important dynamics like refractoriness (the neuron
can be re-excited for an arbitrarily small interval after the last fire event
for sufficiently large input) and spike-rate adaptation (the neuron will not
reduce its firing frequency given a long constant input which is observed
empirically). On the positive side, the model is simple enough so that it is
possible to find its explicit solution for constant input which is a nice and
rare mathematical property to have. This derivation together with more
elaboration on the model’s upsides and downsides can be found in [17].
disable,inlinedisable,inlinetodo: disable,inlinePDE version: weak form, blow-
up, asymptotics, properties, well-posedness

A straightforward extension of the LIF model could be a general nonlinear
integrate and fire model like

$\tau\frac{dV}{dt}=-F(V(t))+G(V(t))I(t)$ (9)

with the same reset condition [27]. A special case of (9) is the quadratic
integrate and fire model

$\tau\frac{dV}{dt}=-C(V(t)-V_{reset})(V-V_{C})+RI(t)$ (10)

which under a change of variable can also be called the theta model of a
biological neuron (with very interesting stability analysis on the unit
circle, also included in [27]).

The second major neuron model which is minimized in terms of the biophysical
mechanisms that bring about action potentials is the spike response model or
SRM [25]. It is represented by the integral equation

$u_{i}(t)=\eta(t-t_{i})+\sum_{j}w_{ij}\sum_{k}\epsilon_{ij}(t-t_{i},t-t_{j}^{k})+\int_{0}^{\infty}\kappa(t-t_{i},s)I(t-s)ds$
(11)

where a neuron indexed by $i$ has a single quantity $u_{i}$ with resting value
$u_{rest}=0$ responding to a presynaptic spike at a fixed time $t_{i}$ and an
external input current $I(t)$ (the last term). The presynaptic neurons with
respect to neuron $i$ are indexed by $j$, have strengths $w_{ij}$ and each one
has spikes indexed by $k$ at times $t_{j}^{k}<t$. The time evolution of the
response to an incoming spike is represented by $\epsilon$ while the form of
the action potential and the after potential is described by $\eta$. An output
spike is then triggered if after the summation of presynaptic neurons and
their spikes at previous times, the value $u_{i}$ reaches a (not necessarily
fixed) threshold $\alpha=\alpha(t-t_{i})$. Here $t_{i}$ is fixed and is the
time of the last action potential

$t_{i}=\max{t_{i}^{k}<t}$

To add refractory period $\Delta$, we can simply set the threshold to a very
high value until time $t_{i}+\Delta$. The notation can be simplified if we
introduce the total PSP quantity

$h(t|t_{i})=\sum_{j}w_{ij}\sum_{k}\epsilon_{ij}(t-t_{i},t-t_{j}^{k})+\int_{0}^{\infty}\kappa(t-t_{i},s)I(t-s)ds$

in which case the model becomes

$u_{i}(t)=\eta(t-t_{i})+h(t|t_{i})$ (12)

Further interpretation of the SRM and its properties is presented in [27].

The integrate-and-fire (LIF) neuron is a special case of the spike response
(SRM) neuron. To see this, we need to take $V_{i}=u_{i}$ as the potential of
the $i$-th neuron and generalize the input current of (8) to

$\tau\frac{du_{i}}{dt}=-u_{i}(t)+RI(t)=-u_{i}(t)+R(\sum_{j}w_{ij}\sum_{k}\beta(t-t_{j}^{k})+I_{i}(t))$
(13)

where $I(t)$ is external current applied to the soma and $\beta$ are
postsynaptic current pulses at times $t_{j}^{k}$. More details about this
comparison can also be found in [27].

####  2.2.2 Firing-rate (rate-based) models

The most detailed way to represent interneuron communication and therefore a
neural network is by preserving the neuron spikes used for communication among
nodes since there could be timing-related patterns like synchronized and
correlated firing and therefore dynamics which can only be captured by this
type of communication. However, if we assume that this is not the case, we can
greatly reduce computational costs due to number of parameters and time scales
as well as the difficulty of any analytic interpretation, all by replacing the
spikes with firing rates.

Variations in the action potential’s duration, amplitude, and shape are very
small and are typically ignored in models of spiking networks which consider
them as identical events occurring at single time instances $0\leq t_{i}\leq
T$. The neuron response function $\rho(t)$ can then be defined as

$\rho(t)=\sum_{i=1}^{n}\delta(t-t_{i})$ (14)

and can be used to obtain any sequence of responses of amplitude defined by
$h(t)$ as

$\sum_{i=1}^{n}h(t-t_{i})=\int_{-\infty}^{+\infty}\rho(t-\tau)h(\tau)d\tau=\rho\star
h$ (15)

We can derive the firing rate (which in this case means a spike-count rate)
from the neuron response function (14) with a continuum assumption as

$r(t)=\lim_{\Delta t\to 0}\frac{1}{\Delta t}\int_{t}^{t+\Delta
t}\rho(\tau)d\tau$ (16)

or alternatively from $\langle\rho(t)\rangle$ averaged from many experimental
trials where the neuron is subjected to the same stimulus. In this sense,
$r(t)$ can also be interpreted as the probability density of firing.

The most basic firing rate model with current dynamics is constructed by [17]:

$\displaystyle\begin{dcases}&\tau_{s}\frac{dI_{s}}{dt}=-I_{s}+\sum_{i=1}^{N_{u}}w_{i}u_{i}=-I_{s}+\mathbf{w}\cdot\mathbf{u}\\\
&v=F(I_{s})=(I_{s}-\alpha)^{+}\end{dcases}$ (17)

where $u$ is used for input firing rate and $v$ for output firing rate with
vectors $\mathbf{u}$ and $\mathbf{v}$ for collections of input and output
units. The $w_{i}$ coefficients are the weights/strengths of the synapses that
receive input firing rates $u_{i}$ and are constants here but generally could
depend on time to account for neuronal plasticity. The first equation then
relates the total input current at the soma $I_{s}$ to the firing rates
received from a total of $N_{u}$ dendrites. It is a simple dynamic equation
with a decay to a steady state at the total input firing rate. The second
equation relates the output firing rate $v$ with the total input current
through an activation function $F$ which is usually taken to be the positive
part function $w_{+}(x):=\max(x,0)$. This is an instant equation with $\alpha$
playing the role of a threshold constant.

It is important to note that (17) is simplified and somewhat collapsed version
of multiple neurophysiological phenomena:

  1. 1.

The amplitude and sign of the synaptic input $i$ are determined by $w_{i}$
which also absorbs the probability of a neurotransmitter release from a
presynaptic terminal.

  2. 2.

The firing rate $u_{i}$ is the time average of the actual spike train as
described in the definition of $r(t)$ and encompasses all the dynamics of the
synaptic conductance of the presynaptic spike like active and passive
properties of the dendritic cables, linear summation of multiple spikes, etc.

  3. 3.

The activation function $F$ takes care of converting units of $I_{s}$ into
units of firing rate through an extra constant multiplier, rendering the
weights $w_{i}$ dimensionless in the final model.

If the second equation is not assumed to be instant, i.e. the dependence of
the output firing rate on the soma current is dynamic, the firing rate model
extends to

$\displaystyle\begin{dcases}&\tau_{s}\frac{dI_{s}}{dt}=-I_{s}+\sum_{i=1}^{N_{u}}w_{i}u_{i}=-I_{s}+\mathbf{w}\cdot\mathbf{u}\\\
&\tau_{r}\frac{dv}{dt}=-v+(I_{s}-\alpha)^{+}\end{dcases}$ (18)

where $\tau_{r}$ is the time scale determining how rapidly $v(t)$ can approach
its steady state when considering constant input or can follow the changes in
$I_{s}$ otherwise.

Finally, considering large difference in the two time scales $\tau_{s}$ and
$\tau_{r}$ and in particular $\tau_{r}\gg\tau_{s}$, we can perform quasi-
steady state approximation with the assumption of $I_{s}\approx const$ with
respect to the second equation to obtain the single equation firing rate model

$\tau_{r}\frac{dv}{dt}=-v+(\mathbf{w}\cdot\mathbf{u}-\alpha)^{+}$ (19)

####  2.2.3 Population models

The firing rates used in the previous section are spike-count rates which work
well under the assumption that the time averaging window is very small with
respect to the time scale of any change in the input. This is not the case in
many practical situations where the reaction times are short. However, we
could use a firing rate also as an average over population instead of average
over time. We then obtain a generalized neuron which represents a local
population aggregate and interactions with other generalized neurons are
equivalent to interactions with other subpopulations.

First and most important population model to look at is the Wilson-Cowan model
[88]. A local population aggregate is assumed to have similar properties and
connections which is confirmed from observations about nearly identical
responses within relatively small volumes of brain tissue and the presence of
local redundancy. Under the additional assumption of random but dense local
connectivity and close spatial proximity, we can neglect spatial interactions
and deal simply with the temporal dynamics of a population aggregate and more
specifically with the portion of the aggregate which is active per unit time.

Considering excitatory and inhibitory populations then leads to two such
unknowns or respectively $E(t)$ and $I(t)$ with 0 as resting potential state.
To obtain the model then, we have to derive expressions for the proportion of
cells which are sensitive (and not refractory) and the cells which receive at
least a threshold excitation. If refractory cells have refractory period of
$r$ msec their proportion and respectively the proportion of sensitive cells
are

$\int_{t-r}^{t}E(t^{\prime})dt^{\prime}\Rightarrow
1-\int_{t-r}^{t}E(t^{\prime})dt^{\prime}$ (20)

and similarly for the inhibitory population. The cells that receive at least a
threshold excitation are given by the functions $S_{e}$ and $S_{i}$ which are
called subpopulation response functions and represent the expected proportion
of cells in a subpopulation which would respond to a given level of excitation
if none of them were originally refractory. These functions can be defined
from a distribution of per-cell thresholds $D(\theta)$ or per-cell synapse
numbers $C(\omega)$ as

$S(x)=\int_{0}^{x(t)}D(\theta)d\theta\qquad\text{or}\qquad
S(x)=\int_{\theta/x(t)}^{\infty}C(\omega)d\omega$ (21)

where per-cell threshold distribution assumes constant synapse number and vice
versa and $x(t)$ is the average excitation that all cells are subjected to
(reasonable assumption for sufficiently large number of synapses). In the
second case all cells with at least $\theta/x(t)$ synapses will receive
threshold excitation which together with the desire to make $S(x)$ a sigmoidal
function (0 and 1 for no and full subpopulations as asymptotic values) is the
reason for the integral bounds.

As a third multiplier in the final equations, we should add an average level
of excitation generated in an excitatory or resp. inhibitory cell at time $t$

$\displaystyle
A_{E}=\int_{-\infty}^{t}\alpha(t-t^{\prime})(c_{1}E(t^{\prime})-c2I(t^{\prime})+P(t^{\prime}))dt^{\prime}$
(22) $\displaystyle
A_{I}=\int_{-\infty}^{t}\alpha(t-t^{\prime})(c_{3}E(t^{\prime})-c4I(t^{\prime})+Q(t^{\prime}))dt^{\prime}$
(23)

which we can obtain if we assume the the effect of stimulation decays with
rate $\alpha$, each cell sums its input synapses with average number of
excitatory and inhibitory synapses given by the constants $c$, and external
input $P(t)$ or $Q(t)$. From (20), (21), and (22), the Wilson-Cowan model now
follows:

$\displaystyle\begin{dcases}&E(t+\tau)=[1-\int_{t-r}^{t}E(t^{\prime})dt^{\prime}]S_{e}A_{e}\\\
&I(t+\tau^{\prime})=[1-\int_{t-r^{\prime}}^{t}I(t^{\prime})dt^{\prime}]S_{i}A_{i}\end{dcases}$
(24)

Here we take the proportion of cells which are both sensitive and above
threshold at time $t$. Both these events are assumed to be uncorrelated with
mutually independent probabilities of each.

disable,inlinedisable,inlinetodo: disable,inlinePDE version: something general
first, dynamical system, Fokker-Planck equation, numerical (FVM)
approximation, asymptotics, evolution, well-posedness, generalized relative
entropy, consequences, QSSA and reduction to 1D Fokker-Planck

Other population models have been developed with LIF as well as SRM neurons as
comprising units. These result in PDEs tracking the change in neuron state
across time and space which in the first case is a membrane potential density
and in the second case is a rafractory density (due to the state of
refractoriness of an SRM neuron) [27].

Besides Wilson-Cowan, there are other integral formulations of the population
dynamics, one particularly by Gerstner:

$A(t)=\int_{-\infty}^{t}P_{I}(t|t^{\prime})A(t^{\prime})dt^{\prime}$ (25)

where $P_{I}(t|t^{\prime})$ is a kernel representing the probability density
that a neuron fires at time $t$ with an input $I(t^{\prime})$ given that it’s
last spike was at $t^{\prime}$ and $A(t)$ is the population activity at a
given time instance [26]. This equation can be derived under the assumptions
that

  1. 1.

The total number of neurons in the population remains constant.

  2. 2.

The neurons exhibit no adaptation, i.e. the state of neuron $i$ depends only
on the most recent firing at time $t_{i}^{\prime}$.

  3. 3.

The neurons are independent on a sufficiently small time scale, i.e. the
number of spikes in the interval $(t,t+\Delta)$ which is of this time scale
for each neuron is an independent random variable. Then by the law of large
numbers, for a sufficiently large network they converge in probability to
their expectation values which are the ones we will consider.

The second assumption implies that the probability of a neuron which has fired
at $t^{\prime}$ to fire again at $\hat{t}\in[t^{\prime},t]$ is given by
$\int_{t^{\prime}}^{t}P_{I}(s|t^{\prime})ds$. Then define the probability that
a neuron which has fired at $t^{\prime}$ survives without firing until $t$ as
$S_{I}(t|t^{\prime})=1-\int_{t^{\prime}}^{t}P_{I}(s|t^{\prime})ds$.

The third assumption helps us consider the model in the thermodynamic limit of
$N\rightarrow\infty$ where $N$ is the size of the population. Then the average
of the fraction of neurons at time $t$ which have fired their last spike at
$\hat{t}\in[t_{0},t_{1}]$, $t_{0}<t$, $t_{1}<t$ is
$\int_{t_{0}}^{t_{1}}S_{I}(t|t^{\prime})A(t^{\prime})dt^{\prime}$ where
$A(t^{\prime})\Delta t^{\prime}$ is the fraction of neurons that have fired
(not necessarily last) in $\hat{t}\in[t^{\prime},t^{\prime}+\Delta
t^{\prime}]$ and $S_{I}(t|t^{\prime})$ of these are expected to survive from
$t^{\prime}$ to $t$ without firing.

Finally, the first assumption together with an assumption that all neurons
have fired at some point in the past (including $-\infty$ for all that haven’t
fired in a finite time) implies that if we extend the lower bound to $-\infty$
and the upper bound to $t$ we will obtain

$1=\int_{-\infty}^{t}S_{I}(t|t^{\prime})A(t^{\prime})dt^{\prime}$ (26)

which will hold $\forall t$ because the number of neurons is constant and all
of them have fired in the interval $[-\infty,t]$.

Differentiating (26) and using the facts that by definition of
$S_{I}(t|t^{\prime})$,
$P_{I}(t|t^{\prime})=-\frac{d}{dt}(S_{I}(t|t^{\prime}))$ and trivially that
$S_{I}(t|t)=1$, we arrive at

$0=S_{I}(t|t)A(t)+\int_{-\infty}^{t}\frac{d}{dt}S_{I}(t|t^{\prime})A(t^{\prime})dt^{\prime}=A(t)-\int_{-\infty}^{t}P_{I}(t|t^{\prime})A(t^{\prime})dt^{\prime}$

The equation (25) is highly nonlinear since the kernel $P_{I}(t|t^{\prime})$
depends nonlinearly on the total input current which in turn depends
nonlinearly on the population activity $A(t)$.

####  2.2.4 Binary models

The absolute minimal neuron model worth mentioning here for completeness is
the McCulloch-Pitts (MCP) neuron [62] also discussed in the ANN introduction.
In fact, it is so biologically unrealistic that it belongs to the separate
category of artificial neurons. This neuron is basically a threshold function
(a threshold $\theta$ gate in terms of digital logic) like

$v=\left\\{\begin{array}[]{ll}1\text{ if }\sum_{i}u_{i}>\theta\\\ 0\text{
otherwise}\end{array}\right.$ (27)

and it is one of the few artificial neurons that are often used in ANNs. It is
also called a perceptron (more specifically when adding weights $w_{i}$ to the
terms in the sum) which is also used as a name of all networks derived from
this neuron (single or multilayer perceptron as MLP). Such networks are the
first generation of neuron network models according to [58].

###  2.3 Extension to networks

Any neuron model studied previously can be coupled with other neurons of the
same or even different type to form a network with shared dynamics like
oscillations, synchronization, and possibly chaotic behavior [19, 2, 33].
However, the number of neurons that we could couple in a computationally
feasible manner depends significantly on the complexity of the neuron model
with multi-compartment models in one end of the spectrum and minimal models in
the other. Studying networks of possibly simplified neurons as principal units
instead of detailed single neuron models is motivated by the fact that there
are phenomena emerging only from couplings of neuron dynamics that cannot be
reproduced otherwise. General cognitive capabilities like perception,
learning, or memory recall are all within this category.

####  2.3.1 Activity - study of the potential dynamics

We will construct the various types of networks (feedforward, recurrent, etc.)
from the firing-rate neurons (17, 19), eventually taking the middle ground in
the choice of minimal neuron models. The basic idea behind the construction of
these networks is shown on figure 5.

Figure 5: Construction of network models from single neuron models using
vector notation.

For the feedforward network we only have to generalize (19) to a vector of
unknowns therefore obtaining a matrix of weights $\mathbf{W}$, a constant
vector of thresholds $\boldsymbol{\alpha}$, and two layers of units -
$\mathbf{u}$ inputs and $\mathbf{v}$ outputs. For layered recurrent networks
with feedforward and lateral connection only (no feedback connections), we
also add interactions within the layer $\mathbf{v}$ of outputs where
$\mathbf{M}$ contains the strengths or weights of the lateral connections. For
separate handling of excitatory and inhibitory interactions which we should
require if we assume Dale’s law (units can either be only excitatory or only
inhibitory with respect to others), we only need to split the units within the
output layer ($\mathbf{v}$) into excitatory $\mathbf{v_{E}}$ and inhibitory
$\mathbf{v_{I}}$ and the lateral connections ($\mathbf{M}$) into
$\mathbf{M_{EE}}$ and $\mathbf{M_{IE}}$ with weights greater than or equal to
0 and $\mathbf{M_{EI}}$ and $\mathbf{M_{II}}$ with weights less than or equal
to 0.

In the case of fully recurrent networks, the inputs $\mathbf{u}$ would be
unknowns depending on the outputs $\mathbf{v}$ through feedback connections.
However, we don’t need to add more equations and a matrix for the weights of
the feedback connections. What we can do instead is use to our advantage the
fact that adding feedback connections would then include all possible types of
connections and there will be no need of differentiating among them anymore.
Thus, we simply collapse all classes into one and consider some of the $v_{i}$
neurons as input units and the rest as hidden units allowing for any possible
connection between each two $v_{i}$ and $v_{j}$. To generalize the resulting
models, we should add a term $\mathbf{U}$ to represent external input
(sometimes called bias or control/somatic input/current) to each neuron.
Finally, from (17) we can derive a second type of recurrent network called
additive network which reflects the dynamics of total input current or
membrane potential rather than the output of a neuron.

Performing all of the constructions above, we obtain the following firing rate
network models:

  1. 1.

Feedforward networks

$\tau_{r}\frac{d\mathbf{v}}{dt}=-\mathbf{v}+(\mathbf{W}\cdot\mathbf{u}-\boldsymbol{\alpha})^{+}$
(28)

  2. 2.

Recurrent networks

    1. (a)

Layered networks with feedforward and lateral connections

$\tau_{r}\frac{d\mathbf{v}}{dt}=-\mathbf{v}+(\mathbf{W}\cdot\mathbf{u}+\mathbf{M}\cdot\mathbf{v}-\boldsymbol{\alpha})^{+}$
(29)

    2. (b)

Excitatory-inhibitory networks

$\displaystyle\begin{dcases}&\tau_{E}\frac{d\mathbf{v_{E}}}{dt}=-\mathbf{v_{E}}+(\mathbf{W_{E}}\cdot\mathbf{u}+\mathbf{M_{EE}}\cdot\mathbf{v_{E}}+\mathbf{M_{EI}}\cdot\mathbf{v_{I}}-\boldsymbol{\alpha})^{+}\\\
&\tau_{I}\frac{d\mathbf{v_{I}}}{dt}=-\mathbf{v_{I}}+(\mathbf{W_{I}}\cdot\mathbf{u}+\mathbf{M_{IE}}\cdot\mathbf{v_{E}}+\mathbf{M_{II}}\cdot\mathbf{v_{I}}-\boldsymbol{\alpha})^{+}\end{dcases}$
(30)

    3. (c)

Fully recurrent networks

      1. i.

Output-based networks

$\tau_{r}\frac{d\mathbf{v}}{dt}=-\mathbf{v}+(\mathbf{U}+\mathbf{M}\cdot\mathbf{v}-\boldsymbol{\alpha})^{+}$
(31)

      2. ii.

Additive networks

$\tau_{s}\frac{d\mathbf{I}}{dt}=-\mathbf{I}+\mathbf{U}+\mathbf{M}\cdot(\mathbf{I}-\boldsymbol{\alpha})^{+}$
(32)

Due to excessive liberty in the interpretation of $\mathbf{U}$, the fully
recurrent networks (31, 32) can also be viewed as generalizations of all the
rest. In particular, interpreting the external input $U_{i}$ as the total
feedforward input for the neuron $v_{i}$ from another layer $\mathbf{u}$, we
can define $\mathbf{U}=\mathbf{W}\cdot\mathbf{u}$ and obtain (29) from (31).
In addition, setting the lateral connection weights to zero
($\mathbf{M}=\mathbf{0}$) can give us a feedforward network (28) and fixing
the sign of the weights $\mathbf{M}$, setting
$\mathbf{U}=[\mathbf{U_{E}},\mathbf{U_{I}}]^{T}$, and assuming the same time
scale $\tau_{E}=\tau_{I}=\tau_{r}$ can achieve the separation of interactions
in (30), all starting from (31). In general, (31) and all models derivable
from it as special cases are models based on the neuron output dynamics. In
contrast, (32) considers the total input current $\mathbf{I}$ as the dynamical
variable which can be interpreted also as thresholded potential with the
neuron output firing rate given as
$\mathbf{v}=(\mathbf{I}-\boldsymbol{\alpha})^{+}$. While (32) can be obtained
from (31) if the decay coefficients (in this case 1) are equal and (31) from
(32) if the weight matrix $\mathbf{M}$ is invertible, the two are not
equivalent in general [37]. Until the end of this section, we will perform
some mathematical analysis on these network models that will reveal
interesting cognitive level phenomena they reproduce and will visualize
different mathematical techniques used for their study.

We can find an explicit solution for a linear version of (29) with
$\boldsymbol{\alpha}=\mathbf{0}$ and through some technical assumptions extend
our conclusions about their stability and cognitive capabilities to the
nonlinear case [17]. For this purpose, it is necessary to assume the the
matrix $\mathbf{M}$ is symmetric. In particular, consider the network

$\tau_{r}\frac{d\mathbf{v}}{dt}=-\mathbf{v}+\mathbf{W}\cdot\mathbf{u}+\mathbf{M}\cdot\mathbf{v}$
(33)

The symmetry of $\mathbf{M}$ implies that its eigenvalues $\lambda_{i}$ are
real and therefore can be ordered. The eigenvectors $\mathbf{e_{i}}$
(satisfying $\mathbf{M}\cdot\mathbf{e_{i}}=\lambda_{i}\mathbf{e_{i}}$)
associated with any two distinct eigenvalues are orthogonal (or even
orthonormal if normalized, i.e.
$\mathbf{e_{i}}\cdot\mathbf{e_{j}}=\delta_{ij}$ for
$\lambda_{i}\neq\lambda_{j}$) and the eigenvectors associated with an
eigenvalue with greater multiplicity are linearly independent (so we can still
choose a basis of orthonormal vectors). This implies that we can express and
replace $\mathbf{v}$ in (33) as

$\mathbf{v}(t)=\sum_{i=1}^{N}c_{i}(t)\mathbf{e_{i}}\Rightarrow\tau_{r}\sum_{i=1}^{N}\frac{dc_{i}}{dt}\mathbf{e_{i}}=-\sum_{i=1}^{N}(1-\lambda_{i})c_{i}(t)\mathbf{e_{i}}+\mathbf{U}$

where $N=\dim(\mathbf{v})$ and split into decoupled ODEs after a dot product
with $\mathbf{e_{j}}$ and using the orthogonality property

$\tau_{r}\frac{dc_{j}}{dt}=-(1-\lambda_{j})c_{j}(t)+\mathbf{e_{j}}\cdot\mathbf{U}\Rightarrow
c_{j}(t)=\frac{\mathbf{e_{j}}\cdot\mathbf{U}}{1-\lambda_{j}}(1-e^{(-\frac{t(1-\lambda_{j})}{\tau_{r}})})+c_{j}(0)e^{(-\frac{t(1-\lambda_{j})}{\tau_{r}})}$

This analytic solution could be used to identify $\lambda_{j}<1$ as the region
of stability where the system approaches the steady state

$\mathbf{v}(t)=\sum_{i=1}^{N}c_{i}(t)\mathbf{e_{i}}\rightarrow\mathbf{v_{\infty}}=\sum_{i=1}^{N}\frac{\mathbf{e_{i}}\cdot\mathbf{U}}{1-\lambda_{i}}\mathbf{e_{i}}$

The steady state solution hence involves the projection of the input vector
$\mathbf{U}$ onto the eigenvectors of the recurrent matrix $\mathbf{M}$ and is
amplified by a factor of $1/(1-\lambda_{i})>1$ along the $i$-th dimension if
$0<\lambda_{i}<1$. This leads to some cognitive capabilities like selective
amplification when a subset of eigenvectors $\mathbf{e_{k}}$ have eigenvalues
$\lambda_{k}$ less than but close to $1$ in comparison to the rest. In this
case, they dominate the sum and the steady state solution becomes

$\mathbf{v_{\infty}}\approx\sum_{k=1}^{K}\frac{\mathbf{e_{k}}\cdot\mathbf{U}}{1-\lambda_{k}}\mathbf{e_{k}}+\sum_{k=K+1}^{N}\frac{\mathbf{e_{k}}\cdot\mathbf{U}}{1-\lambda_{k}}\mathbf{e_{k}}\approx\sum_{k=1}^{K}\frac{\mathbf{e_{k}}\cdot\mathbf{U}}{1-\lambda_{k}}\mathbf{e_{k}}$

and the network amplified the projection of the input into the subspace
spanned by all $\mathbf{e_{k}}$. Another cognitive capability of such networks
is input integration which can be seen if we consider the eigenvalues
$\lambda_{k}=1$ in which the equations for $c_{k}$ become

$\tau_{r}\frac{dc_{k}}{dt}=-(1-1)c_{k}(t)+\mathbf{e_{k}}\cdot\mathbf{U}=\mathbf{e_{k}}\cdot\mathbf{U}\Rightarrow
c_{k}(t)=c_{k}(0)+\frac{1}{\tau_{r}}\int_{0}^{t}\mathbf{e_{k}}\cdot\mathbf{U}(t^{\prime})dt^{\prime}$

Now $c_{k}$ can remember previous input in the sense that if an input
$\mathbf{U}(t)$ is nonzero initially but zero for an extended period
afterwards, all the remaining $c_{i}$ where $\lambda_{i}<1$ will settle to
zero and we will have

$\mathbf{v}(t)\approx\sum_{k=1}^{K}\frac{1}{\tau_{r}}\int_{0}^{t}\mathbf{e_{k}}\cdot\mathbf{U}(t^{\prime})dt^{\prime}+0\approx\sum_{k=1}^{K}\frac{1}{\tau_{r}}\int_{0}^{t}\mathbf{e_{k}}\cdot\mathbf{U}(t^{\prime})dt^{\prime}$

where we have assumed $c_{j}(0)=0$, $\forall j=1,\cdots N$. This means that
the projection integrals will be preserved in the absence of input and could
be considered as sustained activity or remembrance of the integral of a prior
input.

Our specific interest in network models is in such cognitive level phenomena -
the ones that they could reproduce and the way they would explain them. The
activity dynamics of these networks is studied extensively by [17] who
compiled results from numerous papers before. Feedforward networks with the
above definition can calculate coordinate transformations needed in visually
guided reaching tasks [77]. Recurrent networks have the ability to perform
selective amplification, input integration and sustained activity for instance
in maintaining memory trace of the horizontal eye position [78]. Extending to
the nonlinear recurrent network above (mostly to avoid negative firing rates
which are introduced by (33)) retains this selective amplification of input
stimulus, as well as sustained activity, e.g. memory of a perceived input at
the time of uniform current input. In addition, it is shown to perform
selective dampening and therefore could describe both simple and complex cells
responses in visual processing which are respectively selective and
nonselective towards particular position of the stimulus [79]. The selective
amplification and dampening also produce winner-takes-all perception of
multiple stimuli where one of two overlapping stimuli is ignored.

Stability analysis can be performed on (30) using a phase plane portrait and
linearization around an equilibrium point. Study of the parameters can also
reveal the possibility of a Hopf bifurcation [17]. A deeper analysis
elaborates on limitations of the additive inhibitory interactions and
implements a multiplicative effect inhibition [4]. However, we will skip these
here for the sake of more complex tasks and refer the interested reader to
[17]. A model like (30) is an important alternative when more realism is
needed for the inhibitory interactions since it allows for different
(empirically slower) time scale and (empirically larger) magnitude with
respect to the excitatory interactions.

We are mostly interested in (32), the study of which brings to conclusions
also for the potential-based versions of the previous networks (28, 29, 30).
In particular, by careful choice of the external input $\mathbf{U}$ and the
measured output firing rates $\mathbf{v}$ one can partition the neurons into
input units (units with nonzero external input), output units (units with
measured firing rates), and hidden units (all the rest). These additive
networks are also the most widely studied in literature in comparison to (31)
and any specific cases thereof, most typically with consideration of a
constant or clamping input $\mathbf{U}(t)=\mathbf{U}$ (with some recent
studies on oscillatory input as well [91]). In this case we can reuse our old
notation as

$\tau_{r}\frac{d\mathbf{u}}{dt}=-\mathbf{u}+\mathbf{W}\cdot\mathbf{a(u)}+\mathbf{U}$
(34)

where $\mathbf{u}$ is the neuron potential, $\mathbf{W}$ is the recurrent
weight matrix, $\mathbf{a(u)}$ is a general activation function for the firing
rate such that $\mathbf{a}(\mathbf{u})=(a(u_{1}),\cdots,a(u_{n}))^{T}$ and
$\mathbf{U}$ is an external input to the neuron. The existence and uniqueness
of solutions of (32) can easily be established since the activation function
$a(x)=(x-\alpha)^{+}$ and therefore the entire vector field is globally
Lipschitz-continuous with Lipschitz constant of 1. Further results are well
established for discontinuous activation functions of (34). The stability of
(32, 34) plays a pivotal role in models of memory recall, pattern recognition
and sequence prediction. A stable additive network is called an attractor
network (also autoassociative memory, Hopfield network after [40, 41], or
Cohen-Grossberg-Hopfield network after [13]) due to the existence of an
attractor for its dynamics. It is the ultimate goal for this section and
offers a mathematical explanation for major cognitive phenomena like the ones
just mentioned.

Under certain conditions, the most general recurrent network (34) is stable
and can be shown to converge globally through the use of the energy (Lyapunov)
functional

$E(\mathbf{u})=-\frac{1}{2}\mathbf{a(u)}^{T}\mathbf{W}\mathbf{a(u)}-\mathbf{U}^{T}\mathbf{a(u)}+\sum_{i=1}^{N_{u}}\tilde{E}(u_{i})$
(35)

The functional $E(\mathbf{u})$ was introduced by [13] with last term

$\tilde{E}(u_{i})=\int_{0}^{u_{i}}sa^{\prime}(s)ds$

and by [41] with last term

$\tilde{E}(u_{i})=\int_{0}^{a(u_{i})}a^{-1}(s)ds$

where the second last term can be obtained from the first via change of
variable $y=a(s)$ assuming that $a(0)=0$. This energy has negative
semidefinite time derivative under the assumption that the matrix $\mathbf{W}$
is symmetric and that the activation function is monotone nondecreasing
$a^{\prime}(u)\geq 0$. This can be seen as

$\displaystyle\frac{dE(\mathbf{u})}{dt}=\sum_{i=1}^{N}\frac{dE}{du_{i}}\frac{du_{i}}{dt}=\sum_{i=1}^{N}[(-\sum_{j}^{N}w_{ij}a(u_{j})-U_{i}+u_{i})\frac{da(u_{i})}{du_{i}}]\frac{du_{i}}{dt}=$
$\displaystyle=-\sum_{i=1}^{N}a^{\prime}(u_{i})(\frac{du_{i}}{dt})^{2}\leq
0\text{ if }a^{\prime}(u_{i})\geq 0$

in the first case and as

$\displaystyle\frac{dE(\mathbf{u})}{dt}=\sum_{i=1}^{N}\frac{dE}{da(u_{i})}\frac{da(u_{i})}{du_{i}}\frac{du_{i}}{dt}=$
$\displaystyle=\sum_{i=1}^{N}(-\sum_{j}^{N}w_{ij}a(u_{j})-U_{i}+a^{-1}(a(u_{i})))\frac{da(u_{i})}{du_{i}}\frac{du_{i}}{dt}=$
$\displaystyle=-\sum_{i=1}^{N}a^{\prime}(u_{i})(\frac{du_{i}}{dt})^{2}\leq
0\text{ if }a^{\prime}(u_{i})\geq 0$

in the second case. In order to show that the system is stable using LaSalle’s
invariance principle, we also need $E(\mathbf{u})$ to be continuous and
bounded from below, so that we have convergence to an invariant set. This is
the case for any saturating continuously differentiable activation function.
In the case where the activation function is bounded from below but is at
least locally Lipschitz continuous like $a(x)=(x-\alpha)^{+}$, we need to take
the derivative of the energy along the trajectories as

$\displaystyle D^{+}E(\mathbf{u_{0}})=\lim_{h\rightarrow
0^{+}}\inf\frac{E(\mathbf{u_{0}}+h\mathbf{\dot{u}(u_{0})})-E(\mathbf{u_{0}})}{h}$

and replace the Riemann integral with a Radon integral in the definition of
(35). Finally, if the activation function and therefore the energy functional
is only continuous, its derivative can be taken as

$\displaystyle D^{+}E(\mathbf{u_{0}})=\lim_{h\rightarrow
0^{+}}\inf\frac{E(\mathbf{u}(h,\mathbf{u_{0}}))-E(\mathbf{u_{0}})}{h}$

where $\mathbf{u_{0}}=\mathbf{u}(0)$ and $\mathbf{u}(t,\mathbf{u_{0}})$ is the
trajectory with initial condition $\mathbf{u_{0}}$ [32]. Clearly, the previous
linear case with $a(x)=x$ is smooth but not bounded from below which is why it
can play the role of an integrator.

We are paying large attention to the stability analysis of these network
models because it represents a cognitive phenomenon on its own. When the
activity dynamics of a network converge to a few fixed/equilibrium/steady
states without external influence, we can compare this to the cognitive
process of memory recall. Memory recall can then be considered as a part of
the process of general pattern recognition, or even prediction if the
recognition is of spatio-temporal patterns. More generally, attractor networks
or networks whose dynamics settle towards an attractor, can be seen as a way
to store and retrieve information based on the information’s full or partial
content as initial condition (content-based rather than address-based
information retrieval). Any recalled purely spatial memories are equilibrium
points within the phase space containing any possible network activity as a
state. The retrieved spatio-temporal memories (including sequences of states)
are general attractors within this phase space, like limit cycles for
oscillations (chewing, walking, other rhythmic outputs), lines or other
subspaces, or even fractals in the case of chaotic network behavior.

Let us return to the specific network models described here by considering
spatial (static) memories or equilibrium points. If the initial activity is
within the basin of attraction of a particular memory, the memory will be
gradually recalled as the steady state for any further evolution or a process
of pattern matching. A typical autoassociative network like (32) has the
additional properties that it satisfies conditions to make its convergence to
an equilibrium point possible. Its evolution is studied from a given initial
condition $\mathbf{v}(0)$. If this initial condition is within the basin of
attraction of a memory state $\mathbf{v}_{m}$ for $m\in[1;N]$ with $N$ as the
total number of memories stored in the network, it will converge and the
memory recall will be successful. Recalling the wrong memory is another case
of unsuccessful recall in addition to divergence of the state. Any further
analysis of such associative memory is an optimization problem over the
weights where one maximizes the capacity as the total number of recallable
memories $N$ and the measure of the basin of attraction for each memory.
Another problem that occurs in networks making use of the Lyapunov functional
(35) is the occurrence of spurious equilibrium points which are not memories.

A synthesis of a network like (31) to store $N$ memories is derived by [17] in
the case of binary outputs. A memory $\mathbf{v}_{m}$ is an equilibrium point
of (31) and therefore must satisfy

$\mathbf{v}_{m}=(\mathbf{M}\cdot\mathbf{v}_{m}-\boldsymbol{\alpha})^{+}=\mathbf{a}(\mathbf{M}\cdot\mathbf{v}_{m})$

where we have assumed no external input $\mathbf{U}=\mathbf{0}$. We are
interested in constructing the weights of $\mathbf{M}$ and maximizing the
number $N$ of vectors $\mathbf{v}_{m}$ that can simultaneously satisfy this
equation for an appropriate choice of $\mathbf{M}$. Considering binary values,
we can require that each $\mathbf{v}_{m}$ has $\gamma N$ active units, i.e.
units with value 1, and $(1-\gamma)N$ inactive units, i.e. units with value 0.
We can call $\gamma$ the sparseness of each memory since higher $\gamma$
allows for storing more memories with less information for each. To guarantee
the existence of equilibrium points, we need the Lyapunov function (35) which
in turn must be bounded from below (which can be achieved if the activation
function $\mathbf{a}$ saturates) and requires $\mathbf{M}$ to be symmetric.
Taking the weight matrix to be

$\mathbf{M}=\mathbf{K}-\frac{\mathbf{11}}{\gamma N}\text{ where
}\mathbf{K}\cdot\mathbf{v}_{m}=\lambda\mathbf{v}_{m}$ (36)

i.e. $\mathbf{v}_{m}$ are the degenerate eigenvectors of a matrix
$\mathbf{K}$, we get

$\displaystyle\mathbf{v}_{m}=\mathbf{a}(\mathbf{M}\cdot\mathbf{v}_{m})=\mathbf{a}((\mathbf{K}-\frac{\mathbf{11^{T}}}{\gamma
N})\cdot\mathbf{v}_{m})=\mathbf{a}(\lambda\mathbf{v}_{m}-\frac{\mathbf{1}(\mathbf{1}^{T}\cdot\mathbf{v}_{m})}{\gamma
N})=$ $\displaystyle=\mathbf{a}(\lambda\mathbf{v}_{m}-\frac{\mathbf{1}(\gamma
N)}{\gamma N})=\mathbf{a}(\lambda\mathbf{v}_{m}-\mathbf{1})$

Solving this for each component gives the conditions

$0=a(-1)\qquad 1=a(\lambda-1)$

on the activation function which are satisfied for $a(x)=(x-\alpha)^{+}$ for
$\lambda=2$. It remains to find $\mathbf{K}$ that satisfies at least
approximately $\mathbf{K}\approx\lambda\mathbf{v}_{m}$ but also such that $M$
is a symmetric matrix. This can be done under certain assumptions on
randomness of the patterns as is done in [17].

More thorough survey on synthesis/design methods for attractor networks with
predetermined equilibrium points, basins of attraction, convergence speed,
etc. is presented in [64]. Different approaches vary significantly in terms of
computational requirements, capacity, stability, and restrictions. For
instance, the early synthesis by [40] with

$w_{ij}=\sum_{m}(2u_{m}^{i}-1)(2u_{m}^{j}-1)$ (37)

yields a capacity of $0.15N$. Since this construction of the weight matrix is
done for (34), a memory is represented by $\mathbf{u}_{m}$ instead of
$\mathbf{v}_{m}$. We can see that if $u_{m}^{i}=u_{m}^{j}=0$ or
$u_{m}^{i}=u_{m}^{j}=1$ for some $m$ the weight $w_{ij}$ will increase by $1$
and that if they are different it will decrease by $1$. Further synthesis
methods include the outer product method, the projection learning rule, and
the eigenstructure method [64]. All these methods go in the direction of
creating and modifying attractors in the form of instant or gradual learning
which is a subject of the next section.

####  2.3.2 Connectivity - study of the synapse plasticity

The synaptic strength or synaptic plasticity as well as the overall network
connectivity are determined by activity-independent and activity-dependent
mechanisms. It is a widely accepted belief that activity-independent
mechanisms are responsible for the original targeting of axons, layers of
enervation, etc. while the specialization and development of neuronal input
selectivity and cortical mapping are also determined by activity-dependent
mechanisms. This is suggested by experimental results from the hippocampus,
neocortex, and cerebellum areas in the brain that link synaptic strength with
neural activity [17]. Most of the models encompassing some kind of adaptation
(LTP and or LTD) are using some form of Hebbian learning. Examples of non-
Hebbian synaptic plasticity are also possible and could involve learning based
only on the presynaptic or postsynaptic firing. An example for nonsynaptic
plasticity could be intrinsic plasticity of a neuron based solely on it
activity.

A few important considerations valid for any adaptable models when it comes to
plasticity are:

  * •

We should impose certain bounds on the synaptic strengths - a bound or a
synaptic saturation should prevent uncontrolled growth of the synaptic
strengths as well as the possibility of them crossing zero, i.e. becoming
inhibitory from excitatory or vice versa.

  * •

We should find a way to preserve diversity among the strengths through some
synaptic competition since if all of them converge to zero or a maximum value
the information contained in the connectivity is lost.

  * •

The plasticity is realized by updating a synaptic strength according to a
learning rule which can be a dynamical equation or a discrete update step.

Until the end of this section, we will review the most widely used learning
rules providing dynamics for the weights, starting with the simplest but most
impractical and ending with the most practical but highly specialized among
them.

Expanding on the rate-based networks (28) build earlier, we again consider the
weights $w_{ij}=w_{ij}(t)$ as synaptic strengths however this time not as
constant coefficients. Since $w_{ij}$ now depend on the activity, i.e. the
rates $v_{i}$ and $u_{j}$, we need another equation for its dynamics. However,
since generally the learning/connectivity time scale $\tau_{w}$ is a lot
slower than the perception/activity time scale $\tau_{r}$, we can use QSSA and
assume $v_{i}=\mathbf{w_{i}}\cdot\mathbf{u}$ instantly reaches its steady
state in the first equation so that the simplest network model capable of
adaptation is

$\tau_{w}\frac{d\mathbf{w_{i}}}{dt}=v_{i}\mathbf{u}=(\mathbf{w_{i}}\cdot\mathbf{u})\mathbf{u}$
(38)

where the equations are indexed by $i$, i.e. we have preferred index notation
rather than tensor notation. This network uses a basic Hebbian learning rule

$\tau_{w}\frac{d\mathbf{w_{i}}}{dt}=F(v_{i},\mathbf{u})=v_{i}\mathbf{u}$ (39)

since it can be immediately verified that simultaneous activity on the inputs
$\mathbf{u}$ and the output $v_{i}$ has positive effect on the derivative of
$\mathbf{w_{i}}$. We can see now that we can produce a different network model
(same QSSA but different plasticity dynamics) for each learning rule we
select.

Two prominent learning rules of the Hebbian type are the correlation and
covariance rules which we can obtain if we use respectively

$F(v_{i},\mathbf{u})=\langle
v_{i}\mathbf{u}\rangle\Rightarrow\tau_{w}\frac{d\mathbf{w_{i}}}{dt}=\mathbf{Q}\cdot\mathbf{w_{i}}$
(40) $F(v_{i},\mathbf{u})=(v_{i}-\theta_{v})\mathbf{u}\text{ or
}F(v_{i},\mathbf{u})=v_{i}(\mathbf{u}-\mathbf{\theta_{u}})\Rightarrow\tau_{w}\frac{d\mathbf{w_{i}}}{dt}=\mathbf{C}\cdot\mathbf{w_{i}}$
(41)

with the input correlation matrix
$\mathbf{Q}=\langle\mathbf{u}\mathbf{u}\rangle$ and the input covariance
matrix
$\mathbf{C}=\langle(\mathbf{u}-\langle\mathbf{u}\rangle)(\mathbf{u}-\langle\mathbf{u}\rangle)\rangle=\langle\mathbf{u}\mathbf{u}\rangle-\langle\mathbf{u}\rangle^{2}=\langle(\mathbf{u}-\langle\mathbf{u}\rangle)\mathbf{u}\rangle$.
Some interesting observations to make here are that both covariance rules (41)
contain thresholds (output threshold $\theta_{v}$ or input thresholds
$\theta_{u}$) that separate the LTD from LTP activity (negative to positive
derivative, i.e. a nullcline for $\mathbf{w_{i}}$). The first to the two
equations modifies the synapse strengths only if they have nonzero presynaptic
activities while the second only if the postsynaptic activity is not zero.

All learning rules until here are unstable and therefore not very practical.
An improved and stable rule is the Bienenstock-Cooper-Monroe (BCM) rule

$\displaystyle\begin{dcases}&\tau_{w}\frac{d\mathbf{w_{i}}}{dt}=F(v_{i},\mathbf{u})=v_{i}\mathbf{u}(v_{i}-\theta_{v})\\\
&\tau_{\theta}\frac{d\theta_{v}}{dt}=v_{i}^{2}-\theta_{v}\end{dcases}$ (42)

that also supports LTD in addition to LTP through the threshold $\theta_{v}$.
Allowing $\theta_{v}$ to grow more rapidly than $v_{i}$ with
$\tau_{\theta}<\tau_{w}$ is the way to achieve stability here. One could also
use rules that make use of normalization of the weights (thus also taking care
of competition among the synapses and limited resource within a neuron) but
this makes use of global information outside of the scope of the individual
synapse and is no longer biologically meaningful. A stable, competitive, and
local rule is the Oja rule

$\tau_{w}\frac{d\mathbf{w_{i}}}{dt}=F(v_{i},\mathbf{u})=v_{i}\mathbf{u}-\alpha
v_{i}^{2}\mathbf{w_{i}}$ (43)

where $\alpha>0$ is a constant. It is local because $w_{ij}$ only depends on
$v_{i}$, $u_{j}$, and $w_{ij}$ while in a general normalization we would also
need $w_{ik}$ and $u_{k}$ for $k\neq j$ to complete the computation. It is
stable and competitive because if we multiply the equation (43) by
$\mathbf{w_{i}}$ and use $v_{i}=\mathbf{w_{i}}\cdot\mathbf{u}$, the resulting
equation

$\tau_{w}\frac{d|\mathbf{w_{i}}|^{2}}{dt}=2v_{i}\mathbf{w_{i}}\cdot\mathbf{u}-2\alpha
v_{i}^{2}|\mathbf{w_{i}}|^{2}=2v_{i}^{2}(1-\alpha|\mathbf{w_{i}}|^{2})$

has $|\mathbf{w_{i}}|^{2}=\frac{1}{\alpha}$ as a steady state which is finite
(stability) constant (competition - to increase one weight and preserve the
constant we must decrease another weight).

A final and very important learning rule that we will look at is the timing
based rule

$\tau_{w}\frac{d\mathbf{w_{i}}}{dt}=F(v_{i},\mathbf{u})=\int_{0}^{+\infty}H(\tau)v_{i}(t)\mathbf{u}(t-\tau)+H(-\tau)v_{i}(t-\tau)\mathbf{u}(t)d\tau$
(44)

This rule stems from an attempt to relate rate-based Hebbian learning back to
spike-based empirical results where a significant synaptic potentiation was
observed if the presynaptic spike is followed by a postsynaptic spike within a
sufficiently short window[27]. If the postsynaptic spike is followed by the
presynaptic spike, the synaptic strengths are rather depressed and if the
window is too large no change in the synaptic plasticity takes place. The rule
is an approximation where $\tau$ is the difference between the pre- and
postsynaptic firing rate evaluation times. Here $H(\tau)$ describes the weight
change rate as a function of the difference $\tau$ and a good candidate for it
is anything of the form resembling a bounded version of $\frac{1}{x}$,
$x\in\mathbb{R}$. If this function always has the sign of $\tau$, the first
term in (44) represents LTP and the second LTD. It is stable if we add proper
bounds and competitive if we consider spiking network models. Most
importantly, it depends on the order of the pre- and postsynaptic activity.

Although the above models are only feedforward networks that have been
enhanced with plasticity, they are already capable of a number of cognitive
level learning phenomena like the development of ocular dominance in cells of
the primary visual cortex, orientation selectivity, and trace learning [17].
On the other hand, adding lateral connections like (29) can reproduce the
appearance of ocular dominance strips on the output layer and other additional
phenomena.

There could be different ways of adding lateral connections, namely one could
add fixed lateral connections keeping the feedforward connections plastic, fix
the feedforward connections and add plastic lateral connections or have
plasticity for all connections. In many cases using simply a feedforward
network might lead to redundancy in the output layer where each neuron will
obtain the same weights and input selectivity for instance using (43). One way
to deal with this is by adding plastic lateral connections with anti-Hebbian
learning rule (depressing correlated firing on the second layer to reduce
redundancy in all firing there, also observed in nature) in addition to the
feedforward Oja rule

$\displaystyle\begin{dcases}&\tau_{w}\frac{dw_{ij}}{dt}=v_{i}u_{j}-\alpha
v_{i}^{2}w_{ij}\\\ &\tau_{M}\frac{dm_{ij}}{dt}=-v_{i}v_{j}+\beta
m_{ij}\end{dcases}$ (45)

The second equation is the anti-Hebbian basic rule, i.e. the learning rule
(39) with a reversed sign and with an additional $\beta m_{ij}$ term to make
sure all weights don’t converge to zero (the now reversed instability of
Hebbian learning). This model network is stable yet not trivial across neurons
with a suitable choice of a $\beta>0$ constant.

If the synaptic modification is modelled as a discrete rather than continuous
process which would be preferable if the input consists of temporal sequence
of patterns, we could use a discrete learning/update rule

$\mathbf{w_{i}}\rightarrow\mathbf{w_{i}}+\epsilon F(v_{i},\mathbf{u})$ (46)

To obtain the discrete version of any of the continuous learning rules above,
one only needs to use the correct $F(v_{i},\mathbf{u})$. A this point it’s
very easy to distinguish between unsupervised and supervised learning in terms
of the models’ formalism. In the case of unsupervised learning which is what
we have implicitly considered above, $v_{i}$ will be generated by the network
through $v_{i}=\mathbf{w_{i}}\cdot\mathbf{u}$ which can be replaced in
$F(v_{i},\mathbf{u})$ to obtain an unsupervised learning rule. In the case of
supervised learning, $v_{i}$ is already given as a desired output value
$v_{i}(t)=\bar{v_{i}}$ which is then replaced in
$F(v_{i}(t),\mathbf{u}(t))=F(\bar{v_{i}},\mathbf{u}(t))$.

Two important examples where supervised Hebbian learning rules are not
performing well are binary classification through a perceptron neuron and
functional approximation through a firing rate neuron. If we consider a
perceptron model like (27) as a binary classifier, it is possible for the
perceptron to classify inputs perfectly if and only if there exists a
hyperplane separating the input space into half spaces of inputs with each one
of the binary values. This condition is known as the condition of linear
separability. The case of functional approximation is possible when the firing
rate of a neuron is given by a function of a stimulus parameter
$v(s)=\mathbf{w}\cdot\mathbf{f}(s)=\sum_{i}w_{i}f_{i}(s)$. The approximated
function is a linear combination of the functions of the stimulus for the
firing rates if the input neurons $u_{i}$. In both examples, supervised
Hebbian learning rules will not work well because they do not depend on the
actual performance of the network. For this reason we can introduce error-
correcting learning rules where we would start with initial guess for the
weights, compare the actual neuron output $v(u^{k})$ in response to the input
$u^{k}$ with the desired output $\bar{v}^{k}$ and adapt the weights to reduce
the error between the two. Here we have dropped the $i$ index since we
consider one neuron, and have added the $k$ index since we consider a discrete
sequence of inputs in time. The binary classification can be augmented with
the perceptron learning rule

$\mathbf{w}\rightarrow\mathbf{w}+\frac{\epsilon_{w}}{2}(v^{m}-v(\mathbf{u}^{m}))\mathbf{u}^{m}\text{
and }\theta\rightarrow\theta-\frac{\epsilon_{w}}{2}(v^{m}-v(\mathbf{u}^{m}))$
(47)

and the functional approximation with the delta learning rule

$\mathbf{w}\rightarrow\mathbf{w}-\epsilon_{w}\nabla_{w}E$ (48)

where $\theta$ is the threshold from (27) and $E(s)$ is the square error
function which makes the delta rule a an update step in a gradient descent.

A step beyond the rate-based learning rules discussed above and into
plasticity for spiking models is taken by [27] and is recommended for more
biological realism.

##  3 Proposed model

###  3.1 Formulation

We will now construct and motivate a family of rate-based models of increasing
complexity. These have significant morphological difference with the regular
neural network models in the previous section but can, under suitable
conditions, be reduced to or rather interpreted as such networks therefore
attaining their previous results. The emphasis of course is on new and more
general cognitive phenomena that they should describe.

####  3.1.1 Construction of a family of models

As we have mentioned before, it is essential to pick the right features from
the more complex and biologically realistic models while still abstracting
from everything that could be replaced with computationally effective
mechanism. In the model we are going to propose in this section, we do exactly
this by extracting further possible features from the synaptic integration as
well as the distal dendritic tree.

Figure 6: Comparison between a simple circuit with metaconnections and an
actual neuron’s morphology.

Moving in the direction of multi-compartment models, we are interested in
having a separate dynamic for each branch of the dendritic tree. A separate
activation/threshold function $a(\cdot)$ should capture the properties of
propagating an action potential towards the soma while the local PSP
$e(\cdot)$ should be individually determined from the presynaptic input. To
reduce computational and notational complexity, we will restrict ourselves
only to the case of single branching, i.e. a dendrite can only branch once
where it, together with its subbranches, each have their own dynamic. We shall
exclude having an entire cable PDE as such dynamic and instead construct our
own abstractions to capture some minimal properties we need. We will consider
the neuron somas as simple units while axons and unbranched dendrites will be
allowed to be confused together as connections. Finally, subbranches of
dendrites but also some axon terminals will be called metaconnections since
they can be loosely defined as connections to other connections.

We start constructing a new model from (28) by first defining a differential
equation for each unit, connection, and metaconnection, then adding more
specialized dynamics to each equation and finally defining an input-output
layer feedforward network. Also as before, we will then gradually increase the
complexity of the network topology moving to recurrent and fully recurrent
networks. The basic feedforward model with dynamics similar to (28), $N$
inputs, and $M$ outputs reads

$\displaystyle\begin{dcases}&\tau_{r}\frac{dv_{i}}{dt}=-v_{i}+(\sum_{j=1}^{N}w_{j}^{i}c_{j}^{i}-\alpha)^{+}\\\
&\tau_{r}\frac{dc_{j}^{i}}{dt}=-c_{j}^{i}+(\sum_{\begin{subarray}{c}k=1\\\
k\neq j\end{subarray}}^{N}w_{k}^{ji}c_{k}^{ji}+u_{j}-\alpha)^{+}\\\
&\tau_{r}\frac{dc_{k}^{ji}}{dt}=-c_{k}^{ji}+(u_{k}-\alpha)^{+},\ k\neq
j\end{dcases}$ (49)

where $i\in[1;M]$, $j\in[1;N]$, and $k\in[1;N]$. The notation is similar to
the one used for the basic rate-based networks before, namely $v_{i}$ for the
$i$-th output and $u_{j}$ for $j$-th input unit and the $\mathbf{W}$ tensor
represents constant weights. However, there are the dynamics for the pathways
to the output units denoted by $c_{j}^{i}$ for the regular connections from
input unit $j$ to output unit $i$ and by $c_{k}^{ji}$ for the metaconnections
from input unit $k$ to the connection from $j$ to $i$.

Changing the dynamics from firing rate outputs to internal potentials
$u^{i}=\sum_{j=1}^{N}w_{j}^{i}c_{j}^{i}$, $u_{j}^{i}=\sum_{k=1k\neq
j}^{N}w_{k}^{ji}c_{k}^{ji}+u_{j}$, $u_{k}^{ji}=u_{k}$, and $U_{j}$ as internal
potential of the input unit with its firing output rate determined by
$\tau_{r}\dot{u_{j}}=-u_{j}+(U_{j}-\alpha)^{+}$, we obtain

$\displaystyle\tau_{r}\frac{du^{i}}{dt}=\sum_{j=1}^{N}w_{j}^{i}\tau_{r}\frac{dc_{j}^{i}}{dt}=-\sum_{j=1}^{N}w_{j}^{i}c_{j}^{i}+\sum_{j=1}^{N}w_{j}^{i}(\sum_{\begin{subarray}{c}k=1\\\
k\neq j\end{subarray}}^{N}w_{k}^{ji}c_{k}^{ji}+u_{j}-\alpha)^{+}=$
$\displaystyle=-u^{i}+\sum_{j=1}^{N}w_{j}^{i}(u_{j}^{i}-\alpha)^{+}$
$\displaystyle\tau_{r}\frac{du_{j}^{i}}{dt}=\sum_{\begin{subarray}{c}k=1\\\
k\neq
j\end{subarray}}^{N}w_{k}^{ji}\tau_{r}\frac{dc_{k}^{ji}}{dt}+\tau_{r}\frac{du_{j}}{dt}=\sum_{\begin{subarray}{c}k=1\\\
k\neq
j\end{subarray}}^{N}w_{k}^{ji}[-c_{k}^{ji}+(u_{k}-\alpha)^{+}]+\tau_{r}\frac{du_{j}}{dt}=$
$\displaystyle=-\sum_{\begin{subarray}{c}k=1\\\ k\neq
j\end{subarray}}^{N}w_{k}^{ji}c_{k}^{ji}-u_{j}+u_{j}+\sum_{\begin{subarray}{c}k=1\\\
k\neq
j\end{subarray}}^{N}w_{k}^{ji}(u_{k}^{ji}-\alpha)^{+}-u_{j}+(U_{j}-\alpha)^{+}=$
$\displaystyle=-u_{j}^{i}+\sum_{\begin{subarray}{c}k=1\\\ k\neq
j\end{subarray}}^{N}w_{k}^{ji}(u_{k}^{ji}-\alpha)^{+}+(U_{j}-\alpha)^{+}$
$\displaystyle\tau_{r}\frac{du_{k}^{ji}}{dt}=\tau_{r}\frac{du_{k}}{dt}=-u_{k}+(U_{k}-\alpha)^{+}=-u_{k}^{ji}+(U_{k}-\alpha)^{+}$

While $v_{i}$ stands for the output of the $i$-th output neuron in the basic
rate-based models by [17] and therefore also in (49), i.e. the quantity
obtained after applying the activation function, $u^{i}$ is rather the total
accumulated potential of the $i$-th output neuron. Instead of taking the
activation function of the linear summation of inputs, we now take linear
summation of activated inputs similarly to the transition from (31) to (32).
We will now abuse our notation and use $u_{k}$ to imply internal potential of
an input unit. The resulting model reads

$\displaystyle\begin{dcases}&\tau_{r}\frac{du^{i}}{dt}=-u^{i}+\sum_{j=1}^{N}w_{j}^{i}(u_{j}^{i}-\alpha)^{+}\\\
&\tau_{r}\frac{du_{j}^{i}}{dt}=-u_{j}^{i}+\sum_{\begin{subarray}{c}k=1\\\
k\neq
j\end{subarray}}^{N}w_{k}^{ji}(u_{k}^{ji}-\alpha)^{+}+(u_{j}-\alpha)^{+}\\\
&\tau_{r}\frac{du_{k}^{ji}}{dt}=-u_{k}^{ji}+(u_{k}-\alpha)^{+},\ k\neq
j\end{dcases}$ (50)

Notice that the internal potentials depend on the weighted sums of the instant
output firing rates in the above model. We could also replace
$U_{k}:=(u_{k}-\alpha)^{+}$ as external inputs to the system by a second abuse
of notation this time reusing $U_{k}$. Of course, this replacement is only
valid for feedforward and lateral networks since if there are connections to
the input units, their potentials $u_{i}$ become unknowns for the system.

The first major extension besides the more detailed network structure is the
addition of a ”context” term $e$ which replaces the external input term
$a(u_{k})=(u_{k}-\alpha)^{+}$ of (50) with a function defined by

$e(c,u,C_{i})=\begin{dcases}\frac{ca(u)}{h(C_{i})}\qquad\quad\text{if }\quad
h(C_{i})\neq 0\\\ \frac{a(u)}{MN}\qquad\quad\text{ otherwise}\end{dcases}$
(51)

where $C_{i}=\\{c_{i}^{0,1},\dots,c_{i}^{N,M}\\}$ and

$h(C_{i})=\sum\limits_{l=1}^{M}\sum\limits_{\begin{subarray}{c}m=0\\\ m\neq
i\end{subarray}}^{N}c_{i}^{ml}$ (52)

Similarly to the activation function $a(\cdot)$, $e(\cdot)$ can be called
distribution function since it is used to determine the potential distribution
among outputs to different units. For the description of this distribution
function we have used the convention $c_{i}^{l}=c_{i}^{0,l}$. In essence, the
distribution function $e(\cdot)$ determines the potential at a connection by
splitting the total output of its presynaptic unit according to already
accumulated potential through other metaconnections. If no such potential is
present, it will equidistribute the potential which is analogical to
broadcasting of the presynaptic unit. Using the activation function
$a(u)=(u-\alpha)^{+}$, we get the form

$\displaystyle\begin{dcases}&\tau_{r}\frac{du^{i}}{dt}=-u^{i}+\sum_{j=1}^{N}w_{j}^{i}a(u_{j}^{i})\\\
&\tau_{r}\frac{du_{j}^{i}}{dt}=-u_{j}^{i}+\sum_{\begin{subarray}{c}k=1\\\
k\neq
j\end{subarray}}^{N}w_{k}^{ji}a(u_{k}^{ji})+e(u_{j}^{i},u_{j},u_{j}^{0,1},\dots,u_{j}^{N,M})\\\
&\tau_{r}\frac{du_{k}^{ji}}{dt}=-u_{k}^{ji}+e(u_{k}^{ji},u_{k},u_{k}^{0,1},\dots,u_{k}^{N,M}),\
k\neq j\end{dcases}$ (53)

The second major extension that we will consider is the addition of potential
storage if the threshold condition is not met. We can immediately reformulate
(53) as

$\displaystyle\begin{dcases}&\tau_{r}\frac{du^{i}}{dt}=-a(u^{i})+\sum_{j=1}^{N}w_{j}^{i}a(u_{j}^{i})\\\
&\tau_{r}\frac{du_{j}^{i}}{dt}=-a(u_{j}^{i})+\sum_{\begin{subarray}{c}k=1\\\
k\neq
j\end{subarray}}^{N}w_{k}^{ji}a(u_{k}^{ji})+e(u_{j}^{i},u_{j},u_{j}^{0,1},\dots,u_{j}^{N,M})\\\
&\tau_{r}\frac{du_{k}^{ji}}{dt}=-a(u_{k}^{ji})+e(u_{k}^{ji},u_{k},u_{k}^{0,1},\dots,u_{k}^{N,M}),\
k\neq j\end{dcases}$ (54)

The decay of $u^{i}$ is now being subjected to the threshold condition. If
there is no activation $a(u)=(u-\alpha)^{+}=0$, the potential is stored while
it will simply disappear in the previous cases. This will lead to different
short term behavior which is also affected by activity from previous times.
For a system with potential storage, it makes more sense to use the
discontinuous activation/threshold function

$a(u)=\left\\{\begin{array}[]{ll}u\text{ if }u>\alpha\\\ 0\text{
otherwise}\end{array}\right.$ (55)

which reflects any previously stored potential in its output. Thus, when
referring to the form (54), we will assume the discontinuous activation is
used.

Finally, we could block emitting based on highway priority criterion where a
(meta)connection would emit only if its input unit emits, i.e. use an extended
activation function for all the connections and metaconnections

$\bar{a}(c,u)=\left\\{\begin{array}[]{ll}c\text{ if }c>\alpha\text{ and
}u>\alpha\\\ 0\text{ otherwise}\end{array}\right.$ (56)

The final form of the feedforward network using such dynamics becomes

$\displaystyle\begin{dcases}&\tau_{r}\frac{du^{i}}{dt}=-a(u^{i})+\sum_{j=1}^{N}w_{j}^{i}\bar{a}(u_{j}^{i},u_{j})\\\
&\tau_{r}\frac{du_{j}^{i}}{dt}=-\bar{a}(u_{j}^{i},u_{j})+\sum_{\begin{subarray}{c}k=1\\\
k\neq
j\end{subarray}}^{N}w_{k}^{ji}\bar{a}(u_{k}^{ji},u_{k})+e(u_{j}^{i},u_{j},u_{j}^{0,1},\dots,u_{j}^{N,M})\\\
&\tau_{r}\frac{du_{k}^{ji}}{dt}=-\bar{a}(u_{k}^{ji},u_{k})+e(u_{k}^{ji},u_{k},u_{k}^{0,1},\dots,u_{k}^{N,M}),\
k\neq j\end{dcases}$ (57)

The activation function $\bar{a}(\cdot)$ will propagate the potential further
if the more major dendrite is sufficiently excited and will play the standard
role of a threshold function in case this is true.

Before fully motivating the choices of these extensions, we will complete the
model definitions with the other network variants similarly to the way we did
in the extension of neuron to networks models. In order to make it possible to
consider recurrent networks, we need to elaborate more on both the coupling
among equations as well as the indexing we employ to describe the different
scopes of connectivity. From our construction above, metaconnections can
depend only on units, connections only on metaconnections and units, and units
only on connections. Units then can receive only feedforward, lateral and
feedback connections. Feedforward metaconnections can be defined as
connections to feedforward connections from the same layer or to feedback
connections from the next layer and the same applies for the feedback and
lateral (to lateral inward and outward) metaconnections. Hence, there could be
a much larger variety of recurrent network types depending on acceptable
connection types. In order to derive some recurrent network models similar to
the regular networks before, we must make some simplifying assumptions leading
to reduced coupling:

  * •

We will consider only metaconnections to feedforward connections for
simplicity of indexing and phenomenological study.

  * •

Even though we will state a complete definition in terms of connection types,
we will later on exclude lateral and/or feedback connections again with the
goal of reducing the complexity of our study.

Notice that while in (connection only) regular networks removing feedback
connections would completely decouple any layer vector equations, here this is
no longer the case and the evolution of the network would be different if it
comprises two or more layers.

In light of this, both the horizontal and vertical order of the index will be
used to convey information about the indexed connection as shown in figure 7.
The same indexing approach is immediately extendable to metaconnection
indexing due to our first simplifying assumption.

Figure 7: Indexing cases for a connection in a 2x2 network.

Using this indexing convention and the reduced coupling, we can define a
minimalistic layered network with feedforward and lateral connectivity as

$\displaystyle\begin{dcases}&\tau_{r}\frac{du^{i}}{dt}=-a(u^{i})+\sum_{j=1}^{N}w_{j\
}^{\ i}\bar{a}(u_{j\ }^{\
i},u_{j})+\sum_{l=1}^{M}w^{li}\bar{a}(u^{li},u^{l})\\\ &\tau_{r}\frac{du_{j\
}^{\ i}}{dt}=-\bar{a}(u_{j\ }^{\ i},u_{j})+\sum_{\begin{subarray}{c}k=1\\\
k\neq j\end{subarray}}^{N}w_{k\ }^{\ ji}\bar{a}(u_{k\ }^{\ ji},u_{k})+e(u_{j\
}^{\ i},u_{j},u_{j\ }^{\ 0,1},\dots,u_{j\ }^{\ NM})\\\
&\tau_{r}\frac{du^{li}}{dt}=-\bar{a}(u^{li},u^{l})+\bar{e}(u^{li},u^{l},u^{l1},\dots,u^{lM})\\\
&\tau_{r}\frac{du_{k\ }^{\ ji}}{dt}=-\bar{a}(u_{k\ }^{\ ji},u_{k})+e(u_{k\
}^{\ ji},u_{k},u_{k\ }^{\ 0,1},\dots,u_{k\ }^{\ NM}),\ k\neq j\\\
\end{dcases}$ (58)

where $i,l\in[1;M]$ and $j,k\in[1;N]$. An example of an excitatory-inhibitory
network model which will also be used later on to obtain a feature detection
result is one where we have two groups of neurons where each neuron excites
its neighbors within its group and inhibits its neighbors within the other
group

$\displaystyle\begin{dcases}&\tau_{r}\frac{du^{i}}{dt}=-a(u^{i})+\sum_{j=1}^{N}\bar{a}(u_{j\
}^{\ i},u_{j})\\\ &\tau_{r}\frac{du_{j\ }^{\ i}}{dt}=-\bar{a}(u_{j\ }^{\
i},u_{j})+\sum_{\begin{subarray}{c}k=1\\\ k\neq
j\end{subarray}}^{N}\bar{a}(u_{k\ }^{\
ji},u_{k})-\sum_{\begin{subarray}{c}k=1\\\ k\neq
j\end{subarray}}^{N}\bar{a}(v_{k\ }^{\ ji},v_{k})+e(u_{j\ }^{\ i},\dots)\\\
&\tau_{r}\frac{du_{k\ }^{\ ji}}{dt}=-\bar{a}(u_{k\ }^{\ ji},u_{k})+e(u_{k\
}^{\ ji},u_{k},u_{k\ }^{\ 0,1},\dots,u_{k\ }^{\ NM}),\ k\neq j\end{dcases}$
(59)

with three more equations for the $v$ group respectively and all weight
constants set to $1$. Once again, we have to note that in this model we have
ignored Dale’s law for computational simplicity but we could consider it by
separating the two groups into more groups of neurons some of which will have
only inhibitory and others only excitatory outputs. We will look into more
details for this model in section 3.3.1.

A reduced fully recurrent network using the above indexing

$\displaystyle\begin{dcases}&\tau_{r}\frac{du^{i}}{dt}=-a(u^{i})+\sum_{j=1}^{N}w_{j\
}^{\ i}\bar{a}(u_{j\ }^{\
i},u_{j})+\sum_{l=1}^{M}w^{li}\bar{a}(u^{li},u^{l})\\\ &\tau_{r}\frac{du_{j\
}^{\ i}}{dt}=-\bar{a}(u_{j\ }^{\ i},u_{j})+\sum_{\begin{subarray}{c}k=1\\\
k\neq j\end{subarray}}^{N}w_{k\ }^{\ ji}\bar{a}(u_{k\ }^{\
ji},u_{k})+\sum_{l=1}^{M}w^{l\ }_{\ ji}\bar{a}(u^{l\ }_{\ ji},u^{l})+e(u_{j\
}^{\ i},u_{j},\dots)\\\
&\tau_{r}\frac{du^{li}}{dt}=-\bar{a}(u^{li},u^{l})+\bar{e}(u^{li},u^{l},u^{l1},\dots,u^{lM},u^{l\
}_{\ 0,1},\dots,u^{l\ }_{\ NM})\\\ &\tau_{r}\frac{du^{i\ }_{\
j}}{dt}=-\bar{a}(u^{i\ }_{\ j},u^{i})+\bar{e}(u^{i\ }_{\
j},u^{i},u^{i1},\dots,u^{iM},u^{i\ }_{\ 0,1},\dots,u^{i\ }_{\ NM})\\\
&\tau_{r}\frac{du_{k\ }^{\ ji}}{dt}=-\bar{a}(u_{k\ }^{\ ji},u_{k})+e(u_{k\
}^{\ ji},u_{k},u_{k\ }^{\ 0,1},\dots,u_{k\ }^{\ NM}),\ k\neq j\\\
&\tau_{r}\frac{du^{l\ }_{\ ji}}{dt}=-\bar{a}(u^{l\ }_{\
ji},u^{l})+\bar{e}(u^{l\ }_{\ ji},u^{l},u^{l1},\dots,u^{lM},u^{l\ }_{\
0,1},\dots,u^{l\ }_{\ NM})\\\
&\tau_{r}\frac{du_{j}}{dt}=-a(u_{j})+\sum_{i=1}^{M}w^{i\ }_{\ j}\bar{a}(u^{i\
}_{\ j},u^{i})+U_{j}\\\ \end{dcases}$ (60)

can be useful for studying network motifs still distinguishing among the
layers. The fact that we have excluded connections of the kind $u_{j\ }^{\
ji}$ and haven’t excluded the ones of the kind $u^{j\ }_{\ ji}$ is an
engineering-motivated choice which is not of great importance and can be
ignored. The indexing of feedforward and feedback metaconnections $u_{k\ }^{\
ji}$ and $u^{l\ }_{\ ji}$ is a collapsed version of the indexing we would
otherwise need for a metaconnection because we only allow $ji$ to represent a
feedforward connection $u_{j\ }^{\ i}$. The output layer unit $u^{i}$ can be
replaced with $v_{i}$ as before in which case we can use multilayer (more than
two layer) indexing similar to the one above but with $v^{ji}$ instead of
$u^{ji}$ to contain information about the starting reference unit. With such
extended indexing, $u^{i\ }_{\ j}$ would refer to the feedback connection from
$u_{i}$ in the $\mathbf{u}$ layer while $v^{i\ }_{\ j}$ to the feedback
connection from $v_{i}$ in the $\mathbf{v}$ layer.

The functions $a(\cdot)$ and $e(\cdot)$ are similar for each type of
connection with the only difference that we also add extended version of the
distribution function

$\displaystyle\bar{e}(c,u,C_{i})$
$\displaystyle=\begin{dcases}\frac{ca(u)}{\bar{h}(C_{i})}\text{ if
}\bar{h}(C_{i})\neq 0\\\\[5.0pt] \frac{a(u)}{MN}\text{
otherwise}\end{dcases}\quad\bar{h}(C_{i})$
$\displaystyle=\sum\limits_{m=1}^{M}[u^{im}+\sum\limits_{\begin{subarray}{c}n=0\\\
n\neq i\end{subarray}}^{N}(u_{i\ }^{\ nm}+u^{i\ }_{\ nm})]$ (61)

with $C_{i}=\\{u^{i1},\dots,u^{iM},u_{i\ }^{\ 0,1},\dots,u_{i\ }^{\ NM},u^{i\
}_{\ 0,1},\dots,u^{i\ }_{\ NM}\\}$ where for the network (58) we more
specifically have

$C_{l}=\\{u^{l1},\dots,u^{lM}\\}\qquad h(C_{l})=\sum_{m=1}^{M}u^{lm}$

and for the network (60) we more specifically have

$C_{l}=\\{u^{l1},\dots,u^{lM},u^{l\ }_{\ 0,1},\dots,u^{l\ }_{\ NM}\\}\qquad
h(C_{l})=\sum_{m=1}^{M}[u^{lm}+\sum\limits_{\begin{subarray}{c}n=0\\\ n\neq
j\end{subarray}}^{N}u^{l\ }_{\ nm}]$

Finally, a fully recurrent collapsed version can be constructed from (57) if
we arrange all the unknowns in a tensor using the conventions
$u_{k}^{ji}=u_{kji}$, $u_{j}^{i}=u_{j0i}$, $u^{i}=u_{00i}$, $u_{i}=u_{i00}$
and symmetrize the equations

$\displaystyle\begin{dcases}&\tau_{r}\frac{du_{00i}}{dt}=-a(u_{00i})+\sum_{j=1}^{N}w_{j0i}\bar{a}(u_{j0i},u_{j00}),i\neq
0\\\
&\tau_{r}\frac{du_{j0i}}{dt}=-\bar{a}(u_{j0i},u_{j00})+\sum_{\begin{subarray}{c}k=1\\\
k\neq
j\end{subarray}}^{N}w_{kji}\bar{a}(u_{kji},u_{k00})+e(u_{j0i},u_{j00},\dots),j,i\neq
0\\\
&\tau_{r}\frac{du_{kji}}{dt}=-\bar{a}(u_{kji},u_{k00})+e(u_{kji},u_{k00},u_{k01},\dots,u_{kNM}),k,j,i\neq
0,k\neq j\\\\[10.0pt] &u_{000}=1,\ u_{kj0}=0,\forall j\neq 0,\
u_{0ji}=0,\forall j\neq 0,\text{ and }u_{kki}=0,\forall k,i\neq 0\end{dcases}$

now for $i\in[0;M]$, $j\in[0;N]$, and $k\in[0;N]$. Taking
$a(u_{00i})=\bar{a}(u_{00i},1)$, $w_{jji}=0$, and renumbering all the units so
that we don’t distinguish among input, hidden, and output units (therefore
consider connections, metaconnections, and dynamics also to the input units
$u_{i}$), we can see that

$\displaystyle\begin{split}&\tau_{r}\frac{du_{00i}}{dt}=-\bar{a}(u_{00i},u_{000})+1\sum_{j=1}^{N}w_{j0i}\bar{a}(u_{j0i},u_{00j})+0\bar{e}(u_{00i},u_{000},u_{001},\cdots,u_{0NN})\\\
&\tau_{r}\frac{du_{j0i}}{dt}=-\bar{a}(u_{j0i},u_{00j})+1\sum_{k=1}^{N}w_{kji}\bar{a}(u_{kji},u_{00k})+1\bar{e}(u_{j0i},u_{00j},u_{j01},\dots,u_{jNN})\\\
&\tau_{r}\frac{du_{kji}}{dt}=-\bar{a}(u_{kji},u_{00k})+0\sum_{l=1}^{N}w_{lki}\bar{a}(u_{lki},u_{00l})+1\bar{e}(u_{kji},u_{00k},u_{k01},\dots,u_{kNN})\\\
&u_{000}=1,\ u_{kj0}=0,\forall j\neq 0,\ u_{0ji}=0,\forall j\neq 0,\text{ and
}u_{kki}=0,\forall k,i\neq 0\end{split}$

which can be collapsed to

$\begin{split}\tau_{r}\frac{du_{kji}}{dt}=-\bar{a}(u_{kji},u_{00k})+c_{kji}\sum_{l=1}^{N}w_{lki}\bar{a}(u_{lki},u_{l00})+\\\
+b_{kji}\bar{e}(u_{kji},u_{00k},u_{k01},\dots,u_{kNN})+U_{kji}\end{split}$
(62)

with the constant tensors

$\displaystyle b_{kji}=0\text{ if }k=0\text{ else }1$ $\displaystyle
c_{kji}=0\text{ if }j\neq 0\text{ else }1$

the external input tensor $U_{kji}$ and the same predetermined
$u_{000}=1,u_{kj0}=0,\text{ and }u_{0ji}=0,\forall j\neq 0$. The fully
recurrent model (62) is a generalization of all previous models which can be
derived from it under specific assumptions like $w_{jji}=0$.

Notice that all the introduced models don’t consider delay in their dynamics
which is a must for actual physical implementation and also widely studied in
models of neural networks. The reason for this is that we are more interested
in the purely computational side of these models rather than their biological
and physical realism similarly to the way artificial neural networks have
deviated away from their biological counterparts. The models are also static
in terms of network connectivity or synapse plasticity and no learning takes
place in the evolution of the network activity. This is because we need to