# Experimental generation and characterization of partially spatially coherent
qubits

Preeti Sharma$\dagger$, Sakshi Rao$\dagger$ and Bhaskar Kanseri* Experimental
Quantum Interferometry and Polarization (EQUIP), Department of Physics, Indian
Institute of Technology Delhi, Hauz Khas, New Delhi 110016, India
<EMAIL_ADDRESS>
###### Abstract

Partially spatially coherent qubits are more immune to turbulent atmospheric
conditions than coherent qubits, which makes them excellent candidates for
free-space quantum communication. In this article, we report the generation of
partially spatially coherent qubits in a spontaneous parametric down-
conversion (SPDC) process using a Gaussian Schell model (GSM) pump beam. For
this non-linear process, we demonstrate experimentally for the first time, the
transfer of spatial coherence features of the pump (classical) to the
biphotons (quantum) field. Also, the spatial profiles of partially coherent
qubits generated in type-I and type-II non-collinear SPDC process are
experimentally observed and multi-mode nature of partially coherent photons
(qubit) is ascertained. These investigations pave the way toward the efficient
generation of partially spatially coherent qubits with a tunable degree of
spatial coherence, which lead to wide range of applications in frontier areas
such as quantum cryptography, teleportation, imaging, and lithography.

  * Dec 2022

$\dagger$$\dagger$footnotetext: These authors contributed equally to this work

##  1 Introduction

For decades, researchers have been intrigued by the prospect of free-space
quantum communication, which offers secure information exchange using
principles of quantum physics [1, 2, 3]. However there are many challenges
involved in free space communication, for instance, random variation of the
refractive index of the atmosphere, which usually causes intensity
fluctuations and beam wandering [4]. These factors degrade the quality of
received information and limit the transmission distance to a great extent
through strong atmospheric turbulence. An alternative method is to use
partially coherent beams in such conditions, as they are robust against the
deleterious effects of atmospheric turbulence [5, 6]. Additionally, the
speckle pattern in these beams is diminished as a result of the reduced
spatial coherence [7, 8]. Quantum light, such as single photons and entangled
photons, can possess partial spatial coherence properties and may offer
robustness against such turbulent conditions [9, 10]. Partially coherent beams
have gained attention in the quantum realm due to their numerous applications
in quantum communication, quantum metrology, and ghost imaging.

There are many ways for the generation of the quantum state of light based on
the second-order and third-order non-linearity of the medium [11]. The
spontaneous parametric down-conversion (SPDC) process has been a well
established method to generate entangled photons in the spatial, spectral,
polarization, and hybrid domains. It also provides a means to generate true
single photons by heralding one of the twin photons exhibiting quantum
properties. It is a second-order nonlinear process where a photon of higher
frequency is spontaneously down-converted into two photons (twin
photons/biphoton) of lower frequency, named signal and idler. This process is
governed by the energy and momentum conservation, or phase matching conditions
[12]. These twin photons have been utilized in the investigation of
fundamental features such as quantum interference, quantum cryptography, non-
locality, and quantum teleportation. The study of correlations in spatial,
spectral, and polarization domains with coherent pump is well explored in the
literature. Conventionally, the process is categorized as type-0, type-I and
type-II phase matching, based on the polarization of the interacting photons
(pump, signal and idler). The polarization of signal and idler photons are
identical in type-I and type-0 phase matching, with the difference that it is
orthogonal to the pump in type-I case. On the other hand, in type-II phase
matching, the signal and idler photons have orthogonal polarization. The
transfer of the angular spectra of the pump to the biphotons has been
experimentally observed with a collimated and focused coherent pump in type-I
and type-II SPDC [13, 14]. The transverse distribution of biphotons exhibits
strong spatial correlations with the pump as a coherent beam. Lately, there
have been some studies dealing with biphotons generated with a partially
coherent and incoherent pump which revealed the dependence of pump coherence
on profiles (theoretically) [15], quadrature squeezing [16], quantum
interference [17], spatial entanglement [18, 19] and polarization entanglement
[20, 21]. The study of interference in double-slit experiments using quantum
beams has always been a topic of interest owing to quantum lithography [22],
ghost imaging [23, 24], and fundamental studies [25, 26]. The interference
pattern in transverse positions of biphotons has been reported to tailor
polarization entanglement [27]. The coincidence detection in twin photons with
double-slits placed in one [28] or both paths [29] have been reported with the
finite size of coherent pump [30]. The transfer of coherence features from the
pump to biphotons is of utmost relevance. It is demanding to recognize whether
the quantum field generated in SPDC process and to be used in free-space
quantum communication follows the similar statistical features as the well-
documented classical pump beam. To the best of our knowledge, there is no
experimental study so far revealing the transfer of the spatial coherence
properties of the pump (classical field) to biphotons (quantum field) in the
SPDC process. It is thus quite important and timely to address this problem
both experimentally and theoretically, resulting in the generation of tunable
partially spatially coherent qubits.

This work aims to investigate the coherence properties of photons generated
using the SPDC process in a type-II non-collinear configuration. We
experimentally explore the propagation of coherence features from the
classical beam (pump) to the quantum beam (biphotons) in second-order non-
linear interaction. The coherence properties of the pump are initially
determined followed by the investigation of the coherence of photons generated
using SPDC. The results show that the spatial coherence is transferred from
the pump to biphotons and follows the same characteristics as the classical
partially spatially coherent beam. These experimental findings are supported
by the theoretical formulation provided in this paper. The study reveals that
the partially coherent qubits with a tunable degree of spatial coherence to be
used in free space optical communication can be easily generated using the
proposed experimental configuration. The effect of changing the spatial
coherence of the pump with a fixed beam size on the conditional detection
probability and spatial profiles of biphotons generated in type-I and type-II
non-collinear SPDC processes is also studied. The spatial profiles of the
output beam show asymmetry with degradation in the spatial coherence of the
pump. This would also have an impact on the use of these qubits for quantum
lithography and imaging purposes.

This paper is organized as follows: In Section 2, we demonstrate the
generation of a partially coherent GSM beam with a tunable degree of spatial
coherence experimentally. In Section 3, we theoretically and experimentally
investigate the effect of variation in the spatial coherence of the pump on
the Young’s double-slit interference in biphotons generated using the type-II
SPDC process. In Section 4, the spatial profiles of photons in type-I and II
non-collinear geometry are experimentally determined. The conditional
probability distributions showing coincidence counts in type-II non-collinear
SPDC for the entangled positions are determined with variations in pump
coherence. In Section 5, the study is concluded.

##  2 Generation and characterization of partially spatially coherent pump

Figure 1: Experimental setup for the generation of Gaussian Schell model (GSM)
beam including schemes (a) to measure transverse correlation length and (b) to
measure the beam size. Notations: $L$: lens; RGGD: rotating ground glass
diffuser; GAF: Gaussian amplitude filter; DS: double-slit, the dashed arrow
line represents the translation of $L_{1}$ . Figure 2: (a) Plot of visibility
of interference pattern of GSM pump beam with respect to the slit separation
for variation in spot size at the diffuser. (b) Variation of GSM beam size and
transverse correlation length with the translation of lens $L_{1}$. Data
points correspond to experimental data and solid lines correspond to the
fitted theoretical curve.

To generate partially spatially coherent qubits, quantum state generation
along with a mechanism to control the spatial coherence is required. Our
choice for quantum state generation is the SPDC process, in which the
coherence properties of the classical pump beam can be controlled. To simplify
the analysis, a model for partially coherent beam namely the Gaussian Schell
Model (GSM) beam is utilized which offers an analytical formalism to asses the
coherence control [31, 7]. The experimental setup shown in figure 1 is used
for this purpose. A linearly polarized laser beam (405 nm wavelength) is
focused on a rotating ground glass diffuser (RGGD) using a lens $L_{1}$ with a
focal length of 50 mm to generate an incoherent beam. The spot size on the
diffuser is varied with the translation of lens $L_{1}$. Following the van
Cittert-Zernike theorem, the incoherent beam acquires spatial coherence as it
propagates [7]. It is important to note that the properties of generated
biphotons depend on the size of the pump at the non-linear crystal [32, 33,
15]. To investigate the sole effect of coherence, the size of the pump beam
must remain invariant with the spatial coherence control mechanism. To achieve
this, the RGGD is kept at the back focal plane of a lens $L_{2}$ with a focal
length of 150 mm, and only the position of lens $L_{1}$ is varied within few
millimeters range. The secondary incoherent source thus formed at RGGD changes
size with this lens translation and subsequently the emerging beam is
partially spatially coherent with size determined by the lens $L_{2}$. The
beam is then transformed into a Gaussian Schell model beam by passing it
through a Gaussian amplitude filter (GAF) (Figure 1).

The experimental setup to measure the transverse correlation length of the
output beam with the translation of lens $L_{1}$ using a double-slit
interference experiment is shown in Figure 1(a). Using this setup, we measured
the visibility of fringes for four double-slits having same slit width ($a$ =
0.15 mm) and different slit separations ($d$ = 0.25 mm, 0.5 mm, 0.75 mm and 1
mm) (Figure 2(a)). The fringes are captured using a charged couple device
(CCD) camera with exposure time of 1 ms. The visibility of the fringes is
determined from the respective intensity patterns by making a horizontal scan.
It is apparent from figure 2(a) that the visibility of fringes decreases with
the increment in spot size at RGGD. The visibility also decreases with the
increase in slit separation. For identical intensity at the input of both
slits (of DS), the visibility of the interference pattern is a measure of the
absolute value of the degree of coherence $|\gamma({\bf r}_{1},{\bf
r}_{2},\tau)|$ due to any two points ${\bf r}_{1}$ and ${\bf r}_{2}$ and is
expressed as [34, 35, 36]

$\text{Visibility}=|\gamma({\bf r}_{1},{\bf
r}_{2},\tau)|=\Big{|}\frac{2J_{1}(\nu)}{\nu}\Big{|},$ (1)

where $J_{1}(\nu)$ is first order Bessel function,
$\nu=\frac{k_{p}d_{12}a_{s}}{f}$ with $k_{p}=\frac{2\pi}{\lambda_{p}}$ as
wavevector of the beam, $d_{12}$ is distance between two points, $a_{s}$ is
beam radius (in meters) at the surface of diffuser and $f$ is focal length of
lens $L_{2}$. The equation reveals that by measuring the fringe visibility the
degree of spatial coherence can be deduced. Experimental data in figure 2(a)
fits reasonably well with the theoretical expression given in equation (1).
The transverse correlation length $l_{c}$ of the beam can be calculated using
the expression:

$l_{c}=\frac{3.832f}{k_{p}a_{s}}=\frac{3.70503\times
10^{-8}}{a_{s}}\text{meters}.$ (2)

The size of the beam is then reduced by a factor of eight using lenses $L_{3}$
and $L_{4}$ to achieve the resultant GSM pump of desired beam size ($w_{0}$)
compatible with the crystal size. Throughout this paper, the degree of spatial
coherence is quantified using a parameter $A$, which is expressed as
$A=\frac{\delta}{2w_{0}}$ with
$\frac{1}{\delta^{2}}=\frac{1}{l_{c}^{2}}+\frac{1}{4w_{0}^{2}}$ [20]. The
value of $A$ varies from 0 (incoherent beam) to 1 (coherent beam). The range
of interest for this study is $0<A<1$ which is the partially coherent regime.
The pump beam is subsequently incident on the non-linear crystal to generate
the down-converted photons. The variation in transverse correlation length is
determined using equation (2), and the resulting change in pump beam size on
the crystal’s surface with the translation of the lens ($L_{1}$) is analyzed.
In the configuration shown in figure 1(b), the beam profiler is used in place
of the crystal to measure the beam size of the pump at the crystal. Figure
2(b) illustrates the influence of lens translation on beam size and transverse
correlation length, which is in line with the objectives of our experiment.
From this figure, one can clearly deduce that the beam size remains invariant
with lens translation and a pump beam of fixed beam size and tunable
transverse spatial coherence is achieved at the crystal plane.

##  3  Double-slit interference in SPDC photons generated using partially
spatially coherent pump

In the previous section, the spatial coherence features of the pump beam in
SPDC have been controlled arbitrarily. A natural question here is, whether
these coherence properties translate from the classical pump to the quantum
biphotons generated in the SPDC process, as it involves non-linear
interactions. In the following, we investigate the coherence properties of
biphotons using first order interference, which provides an answer to this
question.

###  3.1 Theory

The theoretical treatment of interference using double-slit in photons
generated by type-II non-collinear degenerate SPDC is carried out with the
arrangement shown in figure 3. The output state of biphotons after passing
through the double-slit, to be detected at positions ${\bf r}_{s}$ and ${\bf
r}_{i}$, can be expressed as [37]

$\psi({\bf r}_{s},{\bf r}_{i})=\int d{\bf r}t({\bf r})h_{1}({\bf r}_{s},{\bf
r})h_{3}({\bf r}_{i},{\bf r}),$ (3)

Figure 3: Schematic to measure interference visibility using double-slit in
signal beam generated from GSM pump in Type-II SPDC.

with $h_{1}({\bf r}_{s},{\bf r})$ is impulse response function and $h_{3}({\bf
r}_{i},{\bf r})$ contains contribution from phase matching function, pump
envelope function and impulse response function. The term $t({\bf r})$ is the
transfer function of the double-slit with slit width as $a$ and slit
separation as $d$ placed in the path of the generated beam. These terms are
expressed as [38]

$\begin{split}&h_{1}({\bf r}_{s},{\bf
r})=-\frac{\iota}{\lambda_{s}z_{1}}\exp\Big{(}-\frac{\iota k_{s}}{2z_{1}}({\bf
r}_{s}-{\bf r})^{2}\Big{)},\\\ &h_{3}({\bf r}_{i},{\bf r})=\int\int d{\bf
q}_{s}d{\bf q}_{i}E_{p}({\bf q}_{s}+{\bf q}_{i})\Phi({\bf q}_{s},{\bf
q}_{i})H_{2}({\bf r},{\bf q}_{s}),\\\ &H_{2}({\bf r},{\bf
q}_{s})=-\frac{1}{2\pi}\exp\Big{(}\frac{\iota
k_{s}}{z}(x^{2}+y^{2})\Big{)}\exp\Big{(}-\iota
2\pi(q_{sx}x+q_{sy}y)\Big{)}\exp\Big{(}\iota\pi\lambda_{s}z(q_{sx}^{2}+q_{sy}^{2})\Big{)},\end{split}$
(4)

where $\lambda_{s}$ is wavelength of signal with its wavevector defined by
$k_{s}$ and $\iota$ (iota) represents the imaginary part. The signal and idler
transverse wavevectors are represented by ${\bf
q}_{s}=q_{sx}\hat{x}+q_{sy}\hat{y}$ and ${\bf
q}_{i}=q_{ix}\hat{x}+q_{iy}\hat{y}$, respectively, while $z_{1}$ represents
the distance of the double-slit from the detector and $z$ represents the
distance of the double-slit from the crystal. The electric field of the pump
and the phase matching function are denoted by $E_{p}({\bf q}_{s}+{\bf
q}_{i})$ and $\Phi({\bf q}_{s},{\bf q}_{i})$ respectively. The interference
visibility is determined using the single-photon probability density. It is
the probability of detecting signal/idler without detecting the conjugate
photon of the generated pair. The single-photon probability density can be
expressed as [39]

$\begin{split}p_{1}({\bf r}_{s})=&\int\int d{\bf r}d{\bf
r}^{{}^{\prime}}h_{1}({\bf r}_{s},{\bf r})h_{1}^{*}({\bf r}_{s},{\bf
r}^{{}^{\prime}})t({\bf r})t^{*}({\bf r}^{{}^{\prime}})h_{3}({\bf r}_{i},{\bf
r})h_{3}^{*}({\bf r}_{i},{\bf r}^{{}^{\prime}}).\end{split}$ (5)

We have considered the GSM beam of size $w_{0}$ and transverse correlation
length $l_{c}$ with the expression for pump envelope function expressed in
momentum basis as [20]

$\langle E_{p}({\bf q})E^{*}_{p}({\bf
q}^{{}^{\prime}})\rangle=A_{c}\exp\Big{(}-b_{1}|{\bf q}|^{2}-b_{1}|{\bf
q}^{{}^{\prime}}|^{2}+2b_{2}{\bf q}.{\bf q}^{{}^{\prime}}\Big{)},$ (6)

where $b_{2}=\frac{w_{0}^{2}}{2b_{0}}$,
$b_{1}=\frac{(l_{c}+2w_{0})^{2}}{4b_{0}}$,
$b_{0}=1+\Big{(}\frac{l_{c}}{2w_{0}}\Big{)}^{2}$,
$A_{c}=\Big{(}\frac{A_{p}w_{0}}{2\pi}\Big{)}^{2}$ and $A_{p}$ is
proportionality constant. The phase matching function $\Phi({\bf q}_{s},{\bf
q}_{i})$ is $\operatorname{sinc}(\frac{\Delta{\bf q}L}{2})$ and the Gaussian
approximation of this function leads to [40]

$\Phi({\bf q}_{s},{\bf q}_{i})=\operatorname{sinc}\Big{(}\frac{\Delta{\bf
q}L}{2}\Big{)}\approx\exp\Big{(}\frac{-\alpha L|{\bf q}_{s}-{\bf
q}_{i}|^{2}}{4k_{p}}\Big{)},$ (7)

with $\Delta{\bf q}=\frac{|{\bf q}_{s}-{\bf q}_{i}|^{2}}{2k_{p}}$, $L$ be the
crystal length and $\alpha=0.455$. Substituting all these expressions in
equation (5) yields

$\begin{split}p_{1}({\bf r}_{s})=&\frac{\iota
A_{c}}{2\pi\lambda_{s}z_{1}}\int\int d{\bf r}d{\bf r}^{{}^{\prime}}t({\bf
r})t^{*}({\bf r}^{{}^{\prime}})\exp\Big{[}-\iota k_{s}\Big{(}\frac{({\bf
r}_{s}-{\bf r})^{2}}{2z_{1}}-\frac{\iota k_{s}}{z}|{\bf
r}|^{2}\Big{)}\Big{]}\int\int d{\bf q}_{s}d{\bf q}_{i}\exp\Big{(}2\\\ &\times
b_{2}{\bf q}.{\bf q}^{{}^{\prime}}-b_{1}|{\bf q}|^{2}-b_{1}|{\bf
q}^{{}^{\prime}}|^{2}-\frac{\alpha L|{\bf q}_{s}-{\bf
q}_{i}|^{2}}{4k_{p}}\Big{)}\exp\Big{(}\iota\pi(\lambda_{s}z|{\bf q}|^{2}-2{\bf
q}_{s}.{\bf r})\Big{)}\end{split}$ (8)

For simplicity, the expression is reduced to 1-D (x-axis). The transfer
function of the double-slit is given by [17]

$t(x)=\begin{cases}1&\text{ $\frac{d-a}{2}\leq x\leq\frac{d+a}{2}$ and
$-\big{(}\frac{d+a}{2}\big{)}\leq x\leq-\big{(}\frac{d-a}{2}\big{)}$}\\\
0&\text{otherwise}.\end{cases}$ (9)

Substituting the value of the transfer function for a double-slit, the final
expression for single photon probability density is obtained. The final
expression is numerically integrated to get the intensity profile at the
detector position. The intensity depends on the slit width, slit separation,
the degree of spatial coherence of the pump beam, and the distance of the
double-slit from the crystal and the detector, respectively.

###  3.2 Experiment

Figure 4: Experimental setup for the double-slit interference in photons
generated using non-collinear degenerate type-II SPDC with GSM pump beam.
Notations: $L_{5}$: lens; DM: dichroic mirror; IF: interference filter; DS:
double-slit. Figure 5: Intensity profile of interference fringes obtained for
the signal beam captured using EMCCD camera with variation in spatial
coherence of pump ($A$). Markers with error bars correspond to experimental
data and solid lines are fitted theoretical curves.

The experimental setup to investigate the double-slit interference of photons
produced in the non-collinear degenerate type-II SPDC process is depicted in
figure 4. The generated partially spatially coherent pump beam (see Section 2)
is incident on type-II beta Barium Borate (BBO) crystal which generates down-
converted photons represented by two cones of orthogonal polarizations. The
residual pump beam is blocked using a dichroic mirror (DM). The interference
filter (IF) with a central wavelength at 810 nm and bandwidth of 10 nm is used
for the spectral selection of biphotons. The selection of the vertically
polarized signal beam for the double-slit interference is done by filtering
the idler beam from the output of type-II SPDC. This is done with the aid of a
polarizer placed at $90^{\circ}$ with respect to the pass axis which separates
the ordinary polarized beam (signal) from the extraordinary beam (idler). The
ordinary beam (signal) is then directed through a double-slit of 0.15 mm slit
width ($a$) placed before a 200 mm focal length lens ($L_{5}$). Using an
electron-multiplying charge-coupled device (EMCCD) set in a photon counting
setting, 20000 images are accumulated with an exposure time of 20 ms to record
the interference pattern of the output beam. A zero point or reference point
is designated to a pixel position of the signal in the $x$ direction. Image
processing is used to extract the intensity of each individual pixel in an
image.

Figure 6: Plot of visibility of interference patterns of photons vs slit
separation (of double-slit) for different values of spatial coherence of pump.
Data points represent experimental observation and solid lines correspond to
the fitted theoretical curve.

The experimentally observed intensity at interference plane is normalized and
fitted with the theory as shown in figure 5. This figure shows the
interference pattern for various values of pump spatial coherence with double-
slit having 0.25 mm slit separation ($d$). The achieved interference patterns
degrade with the decrease in the spatial coherence of the pump (A). As
indicated by the solid curves, the theoretically expected patterns using
equation (8) are found in good agreement with the observed data points
(markers with error bars). The resulting visibilities of the interference
patterns for different values of spatial coherence of the pump are obtained
for changing the slit separation (0.25 mm, 0.5 mm, and 0.75 mm) with the same
slit width (0.15 mm) and are plotted in figure 6. The experimentally observed
values in the figure are fitted with a Gaussian function. It is evident that
when the slit separation increases, the visibility of the interference pattern
diminishes, which is consistent with the theoretical expression within the
experimental errors. Furthermore, the visibility decreases with a decrease in
the spatial coherence of the pump. The results show that the correlations in
photons generated using the GSM pump beam in the SPDC process are directly
proportional to the spatial coherence of the pump beam. This is the key result
of the paper as it shows that the spatial coherence is transferred from the
pump (classical) to biphotons (quantum field). Thus, the partially spatially
coherent qubits with a tunable degree of spatial coherence are generated which
could be used in free space optical communication due to robustness against
atmospheric turbulence.

##  4  Transverse distribution of twin photons generated in non-collinear
configuration

In addition to the effect of the spatial coherence of the GSM pump on double-
slit interference, the distribution of generated photons in the transverse
plane is observed for type-I and type-II non-collinear configurations. The
spatial profiles are measured using the scheme shown in figure 4, with the
exception of the polarizer and double-slit.

Figure 7: (a) Simulated and (b) experimentally observed spatial distribution
of output beam generated in type-II SPDC for different values of spatial
coherence of pump ($A$).

###  4.1 Type-II SPDC

The biphotons are generated in the type-II SPDC process with a GSM pump using
beta Barium Borate (BBO) crystal in a non-collinear configuration. The BBO
crystal of 2 mm crystal length interacts with a horizontally polarized pump
with a wavelength of 405 nm, which produces frequency degenerate biphotons in
type-II non-collinear geometry. The spectral selection of the biphotons using
an interference filter (IF) with a bandwidth of 10 nm and a central wavelength
of 810 nm is the same as is discussed in the former section. Utilizing a lens
$L_{5}$ with a 200 mm focal length, the output is collimated before being
detected by an EMCCD camera (Figure 4).

Figure 8: Conditional detection probability of idler in the horizontal
direction with variation in spatial coherence of pump. The vertical pixel
position for both signal and idler beams, and the horizontal pixel position of
the signal are fixed. The markers represent experimental values and are fitted
with theoretical curves (solid lines).

Ordinary signal (V polarized) and extraordinary idler (H polarized) photons
are created by the extraordinary pump (H polarized). The extraordinary pump
and idler encounter spatial walk-off in this configuration. By using a photon
counting setting with an exposure time of 10 ms and integrating 20,000 images,
the experimental profiles are acquired (figure 7(b)). The corresponding
simulated profiles are shown in figure 7(a). In the profiles, the two rings
indicate ordinary and extraordinary polarizations. The location of the
polarization-entangled photons is at the points where these two rings overlap.
With the degradation of the spatial coherence of the pump, the rings broaden
asymmetrically due to the spatial walk-off of the pump beam and idler beam and
the overlap region becomes fuzzier. This is consistent with the theoretical
formulation within experimental limits [15].

The conditional probability of idler [41] with signal fixed at a particular
position is also calculated and is shown in figure 8. The coincidence counts
between $i$ and $j$ pixel position in captured images is calculated using
$C=\langle n_{i}n_{j}\rangle-\langle n_{i}\rangle\langle n_{j}\rangle$ [42]
where $n_{i}$ and $n_{j}$ are average counts at position $i$ and $j$
respectively. The signal is fixed at pixel position $x_{s},y_{s}$ and $y$
position of idler ($y_{i}$) is also fixed . These positions are the two points
at which polarization entanglement can be observed. The scan of coincidence
counts of idler along $x$ axis with signal conditioned at one of the entangled
positions is fitted with Gaussian distribution. The increase in full-width
half maxima (FWHM) of the fitted curve indicates the increase in the area of
correlation at the observational point. It is clear from figure 8 that with a
decrease in spatial coherence of the pump, the strength (conditional
probability) of correlations decreases and the area of spatial correlation
broadens (increase in FWHM). This broadened distribution infers correlation
with more number of spatial modes (multi-mode nature) as the partial coherence
of the signal beam increases. Owing to this multi-mode nature, alike classical
partially coherent beams, these partially coherent qubits are expected to show
robustness against atmospheric turbulence and higher detection probability on
free-space propagation [6]. Similarly, with a rise in coherence of pump, the
conditional probability increases and the distribution becomes narrower,
resulting in stronger correlations with fewer spatial modes (single mode
nature) of coherent beams. This is another significant outcome of this study.

###  4.2 Type-I SPDC

In this case, Bismuth Borate (BiBO) crystal is employed as non-linear medium
to generate the biphotons in the type-I degenerate SPDC process in non-
collinear geometry. The spatial profile of these biphotons is captured using
an EMCCD camera using the experimental setup depicted in figure 4, excluding
the polarizer and double-slit, and with bi-axial BiBO (type-I) crystal
replacing BBO (type-II) crystal. The extraordinary pump (H polarized)
generates ordinary polarized (V polarized) photons. The simulated and
experimental results of spatial profile in non-collinear type-I SPDC with
variation in the degree of spatial coherence of pump are shown in figure 9(a)
and (b) respectively. The ring displays the biphoton distribution along the
pump’s transverse axis. The profiles are recorded for various values of the
pump’s spatial coherence with lens translation $L_{1}$. A careful look into
the profiles reveals the broadening of SPDC rings in a particular direction,
along with a decrease in the transverse correlation length of the pump beam.
This verifies the theoretical framework proposed in our earlier study that
spatial correlations in twin photons deteriorate with the loss in spatial
coherence of the pump [15], and shows asymmetric broadening in the opposite
direction of spatial walk-off of the pump. These results show that one can
generate a robust beam with variable spatial coherence, which could be
applicable for communication purposes. The study also reveals that the partial
spatial coherence should be optimized carefully so that it does not affect
other features such as polarization entanglement to a great extent.

Figure 9: Spatial profiles of biphotons generated in type-I non-collinear SPDC
for different values of spatial coherence of pump ($A$). (a) Simulated results
(b) Experimental observations.

##  5 Conclusion

The article illustrates an experimental method to generate the spatial
coherence controlled qubits using type-I and type-II non-collinear SPDC
processes. The classical double-slit interference method is used to
characterize the spatial coherence of photons generated in type-II non-
collinear SPDC pumped with a GSM beam having a tunable degree of spatial
coherence. The experimental results are verified with the theoretical
formulation for different values of slit separation and variation in spatial
coherence of the pump for a fixed beam size. This study shows the direct
dependence of the spatial coherence of the pump on the visibility in the
interference pattern of biphotons. In addition, the spatial profiles of down-
converted photons in type-I and type-II SPDC in non-collinear geometry are
experimentally observed with a partially spatially coherent pump beam. The
asymmetric broadening of SPDC rings in the direction opposite to the spatial
walk-off of the extraordinary beam clearly depicts the significant impact of
the spatial coherence of the pump on generated photons. The study of the
conditional probability of spatial correlations at the location of
polarization entanglement shows the enhancement in spatial modes with
degradation in pump coherence. This implies that the strong correlation at two
conjugate points expected with a coherent pump decreases with a decrease in
the spatial coherence of the pump. The polarization entanglement would also be
impacted by the decreased spatial correlations and multiple spatial modes in
type-II SPDC. These findings are crucial in selection of appropriate
parameters of a partially coherent pump beam that will produce entangled
photon sources while maintaining a high degree of spatial correlations. The
robust nature of these partially spatially coherent qubits against atmospheric
turbulence makes them useful for long-distance free-space communication,
quantum teleportation, and reduction in speckles for quantum imaging.

The authors acknowledge the funding received from the Department of Science
and Technology (DST), India under the QUEST scheme
(DST/ICPS/QUEST/Theme-I/2019) to carry out this work. Author PS and SR
acknowledge financial support from the Council of Scientific and Industrial
Research (CSIR) India and DST INSPIRE, respectively for Senior Research
Fellowship (SRF). The authors are grateful to Hemant Kumar Singh for his
suggestions in the experimental section involving pump generation.

Conflicts of interest The authors declare that they have no conflict of
interest.

Data Availability Statement The data that support the findings of the study
are available upon request from the authors.

## References

## References

  * [1] Peng C Z, Yang T, Bao X H, Zhang J, Jin X M, Feng F Y, Yang B, Yang J, Yin J, Zhang Q et al. 2005 Phys. Rev. Lett. 94 150501
  * [2] Aspelmeyer M, Bohm H R, Gyatso T, Jennewein T, Kaltenbaek R, Lindenthal M, Molina-Terriza G, Poppe A, Resch K, Taraba M et al. 2003 Science 301 621–623
  * [3] Jin X M, Ren J G, Yang B, Yi Z H, Zhou F, Xu X F, Wang S K, Yang D, Hu Y F, Jiang S et al. 2010 Nat. Photonics 4 376–381
  * [4] Churnside J H and Lataitis R J 1990 Appl. Opt. 29 926–930
  * [5] Dogariu A and Amarande S 2003 Opt. Lett. 28 10–12
  * [6] Gbur G 2014 J. Opt. Soc. Am. A 31 2038–2045
  * [7] Mandel L and Wolf E 1995 Optical coherence and quantum optics (Cambridge university press)
  * [8] Lajunen H and Saastamoinen T 2011 Opt. Lett. 36 4104–4106
  * [9] Shekel R, Lib O, Sardas A and Bromberg Y 2021 OSA Continuum 4 2339
  * [10] Jha A K and Boyd R W 2010 Phys. Rev. A 81 013828
  * [11] Boyd R W 2020 Nonlinear optics (Academic press)
  * [12] Gerry C, Knight P and Knight P L 2005 Introductory quantum optics (Cambridge university press)
  * [13] S K Lee P, P van Exter M and P Woerdman J 2005 Phys. Rev. A 72 033803
  * [14] Ramírez-Alarcón R, Cruz-Ramírez H and U’Ren A B 2013 Laser Phys. 23 055204
  * [15] Kanseri B and Sharma P 2020 J. Opt. Soc. Am. B 37 505–512
  * [16] Rao S, Sharma P and Kanseri B 2022 J. Opt. Soc. Am. B 39 2280–2288
  * [17] Joshi S and Kanseri B 2020 Optik 217 164941
  * [18] Defienne H and Gigan S 2019 Phys. Rev. A 99 053831
  * [19] Zhang W, Fickler R, Giese E, Chen L and Boyd R W 2019 Opt. Express 27 20745–20753
  * [20] Sharma P, Pathak N K and Kanseri B 2021 Results Phys. 27 104506 ISSN 2211-3797
  * [21] Zhang W, Xu D and Chen L 2022 arXiv preprint arXiv:2211.00841
  * [22] Shimizu R, Edamatsu K and Itoh T 2003 Phys. Rev. A 67(4) 041805
  * [23] Strekalov D, Sergienko A, Klyshko D and Shih Y 1995 Phys. Rev. Lett. 74 3600
  * [24] Li J, Gao W, Qian J, Guo Q, Xi J and Ritz C H 2019 Sensors 19 192
  * [25] Gatti A, Brambilla E and Lugiato L A 2003 Phys. Rev. Lett. 90(13) 133603
  * [26] Zhang D J, Wu S, Li H G, Wang H B, Xiong J and Wang K 2017 Sci. Rep. 7 1–7
  * [27] Fran ça Santos M, Milman P, Khoury A Z and Souto Ribeiro P H 2001 Phys. Rev. A 64(2) 023804
  * [28] Strekalov D V, Sergienko A V, Klyshko D N and Shih Y H 1995 Phys. Rev. Lett. 74(18) 3600–3603
  * [29] D’Angelo M, Chekhova M V and Shih Y 2001 Phys. Rev. Lett. 87(1) 013602
  * [30] Souto Ribeiro P H 1997 Phys. Rev. A 56(5) 4111–4117
  * [31] Gori F, Santarsiero M, Borghi R and Ramírez-Sánchez V 2008 J. Opt. Soc. Am. A: Optics and Image Science, and Vision 25 1016 – 1021
  * [32] Lee J C and Kim Y H 2016 Opt. Comm. 366 442–450
  * [33] Jeronimo-Moreno Y and Jáuregui R 2014 J. Opt. 16 065201
  * [34] Kanseri B and Singh H K 2020 Optik 206 163747 ISSN 0030-4026
  * [35] Thompson B J and Sudol R 1984 J. Opt. Soc. Am. A 1 598–604
  * [36] Pearson B J, Ferris N, Strauss R, Li H and Jackson D P 2018 OSA Contin. 1 755–763
  * [37] Abouraddy A F, Saleh B E A, Sergienko A V and Teich M C 2002 J. Opt. Soc. Am. B 19 1174–1184
  * [38] Saleh B E A, Abouraddy A F, Sergienko A V and Teich M C 2000 Phys. Rev. A 62(4) 043816
  * [39] Ribeiro P H S, Pádua S and Monken C H 1999 Phys. Rev. A 60(6) 5074–5078
  * [40] Chan K W, Torres J P and Eberly J H 2007 Phys. Rev. A 75(5) 050101
  * [41] Defienne H, Reichert M and Fleischer J W 2018 Phys. Rev. Lett. 120(20) 203604
  * [42] Avella A, Ruo-Berchera I, Degiovanni I P, Brida G and Genovese M 2016 Opt. Lett. 41 1841–1844


# SWITCHING ACTIVITY IN AN ENSEMBLE OF EXCITABLE NEURONS

A. G. Korotkov  S. Y. Zagrebin  G. V. Osipov

###### Abstract

The paper proposes two dynamical systems based on the generalized Lotka-
Volterra model of three excitable elements interacting through excitatory
couplings. It is shown that for some values of the coupling parameters in the
phase space of systems, there are heteroclinic cycles containing three or six
saddle equilibrium states and heteroclinic curves connecting them. Under
certain external stimuli that transfer the system from a stable zero
equilibrium state to a small neighborhood of the heteroclinic cycle, the phase
trajectory will alternately visit the neighborhood of saddle equilibrium
states (possibly more than once), after which it will return to its initial
state. The described behavior is proposed to be used to simulate switching
activity in neural ensembles. Different transients are determined by different
external stimuli. The passage of the phase point of the system near the saddle
equilibrium states included in the heteroclinic circuit is proposed to be
interpreted as activation of the corresponding element.

##  1 Introduction

Often, dynamic processes in neural ensembles are sequential switches of
activity between individual neurons and/or groups of neurons. This type of
activity is associated with many physiological functions of the nervous
system. The generation and propagation of excitation sequences between neurons
or groups of them play a crucial role in the brain functioning. For example,
sequential activity is observed in neural networks of sensory and motor
systems of animals [1]. In addition, the regions of the brain responsible for
the reproduction of singing in the brain of birds is the one that generates
sequences of bursting activity [2]. A number of examples where excitation
sequences play a key role can be continued [3, 4].

Currently, an urgent task is to study sequential neural activity through the
theory of dynamical systems. In [1] it was hypothesized that sequential
activity may be a consequence of the presence of a stable heteroclinic contour
in the phase space of a dynamic system modeling the dynamics of a neural
network [5, 6, 7, 8]. According to it, the basis of the transient dynamics is
the winnerless competition principle [9, 10, 11]. In short, the essence of
this principle is the existence of a stable heteroclinic contour between
trajectories of saddle type (these can be saddle equilibrium states, saddle
limit cycles, etc.) [12, 13, 14, 15]. When the representative point is in a
small neighborhood of the saddle trajectory, it is assumed that a certain
element of the network (neuron or group of neurons) corresponding to this
saddle trajectory is activated. Thus, a stable heteroclinic contour in phase
space is a mathematical image of sequential switching of activity in neural
motifs. It was shown in [12] that such a contour exists in ensembles with
inhibitory couplings.

Traditionally, in dynamical systems, the main attention is paid to asymptotic
modes (at $t\rightarrow\infty$), and transition processes are discarded.
However, in neural motifs, the transition process is often more important than
the attractor [16, 17, 18]. This fact brings some specifics to the study of
models of neural ensembles: for example, it may matter whether it is a stable
equilibrium state, a node or a focus (this difference is usually neglected
when studying only the asymptotic properties of a dynamical system). The
presence/absence of saddle limit cycles, etc., may also be important.

The main attention in the work will be paid to the study of transition
processes associated with the presence of heteroclinic contours: closed
contours consisting of one-dimensional heteroclinic trajectories connecting
saddle special trajectories. In the models under consideration, these are
saddle equilibrium states.

Mathematical models of neural ensembles can be divided into two types:
microscopic level models (cell level) and large-scale models (neural ensembles
level). In the first type models each neuron of an ensemble is modeled by a
system of differential or difference equations, and connections are
established between neurons (which can also be modeled by a system of
differential or difference equations).

The second type includes Wilson-Cowan models and generalized Lotka-Volterra
models. These models describe the interaction between populations of neurons.
In this paper, a generalized Lotka-Volterra model is used to model a neural
motif.

##  2 The first type system

In this section we will study the dynamics of the activity of neural ensembles
consisting of elements modeled by the following equation

$\dot{\rho}=\rho(\rho-1).$ (1)

Here and further, the variable $\rho$ describes the activity of a neuron. The
neuron is inactive when $\rho=0$ and active when $\rho>0$. The maximum
activity takes place at $\rho=1$. The paper assumes to consider the interval
of change of $\rho$ from 0 to 1. The dynamics of the system described by the
Eq. (1) is trivial. The equilibrium state $\rho_{1}=0$ is stable. The
equilibrium state $\rho_{2}=1$ is unstable. At $\rho=0$, the individual neuron
is in an unexcited (inactive) state. In the study of the dynamics of neural
ensembles, it will be assumed that all the neurons of the ensemble are not
excited at the initial moment of time. The excitation of neurons will occur
due to the supply of external stimuli.

We will analyze the dynamics of the ensemble of their three elements. The
dynamical system in this case is a special case of the generalized Lotka-
Volterra model and has the form

$\begin{cases}\dot{\rho}_{1}=\rho_{1}(\rho_{1}-1+\alpha\rho_{2}+\beta\rho_{3})+a\,\delta(\rho_{1})\\\
\dot{\rho}_{2}=\rho_{2}(\rho_{2}-1+\alpha\rho_{3}+\beta\rho_{1})+b\,\delta(\rho_{2})\\\
\dot{\rho}_{3}=\rho_{3}(\rho_{3}-1+\alpha\rho_{1}+\beta\rho_{2})+c\,\delta(\rho_{3})\end{cases},$
(2)

where $\delta(x)$ is Dirac delta function. $a,b,c,$ parameters are amplitudes
of the external stimulus on each element. $\alpha$ and $\beta$ parameters
define the couplings between the elements. The paper will consider only the
case of $\alpha>0$, $\beta>0$, which corresponds to the case of excitatory
couplings.

###  2.1 A study of a two-dimensional system on an invariant plane
$\rho_{3}=0$

Since one of the conditions $\rho_{i}=0$ is met, the $i$-th equation of Eqs.
(2) turns into an identity, then the three planes $\rho_{i}=0$ are invariant.
To understand the dynamics in a three-dimensional space, we first consider the
behavior of two-dimensional systems on the invariant planes $\rho_{i}=0$.

On the invariant plane $\rho_{3}=0$ Eqs. (2) takes the form

$\begin{cases}\dot{\rho}_{1}=\rho_{1}(\rho_{1}-1+\alpha\rho_{2})\\\
\dot{\rho}_{2}=\rho_{2}(\rho_{2}-1+\beta\rho_{1})\end{cases}.$ (3)

This system describes the interaction of two interconnected elements.

The equilibrium states of Eqs. (3) and their eigenvalues are given in Table 1.

Equilibrium state | Eigenvalues
---|---
$O_{1}(0,0)$ | $-1$ | $-1$
$O_{2}(1,0)$ | $1$ | $\beta-1$
$O_{3}(0,1)$ | $1$ | $\alpha-1$
$O_{4}(\frac{\alpha-1}{\alpha\beta-1},\frac{\beta-1}{\alpha\beta-1})$ | $1$ | $\frac{(1-\alpha)(\beta-1)}{\alpha\beta-1}$
Table 1: Steady states of Eqs. (3) and it’s eigenvalues.

Let us define the types of equilibrium states of Eqs. (3) depending on the
values of the $\alpha$ and $\beta$ parameters.

  1. 1.

$O_{1}$ is always a stable node;

  2. 2.

In the region A ($0<\alpha,\beta<1$), the equilibria $O_{2}$ and $O_{3}$ are
saddles, $O_{4}$ is an unstable node.

  3. 3.

In the region B ($0<\beta<1<\alpha$), $O_{2}$ is a saddle, $O_{3}$ is an
unstable node, $O_{4}$ has one negative coordinate (at $\alpha\beta\neq 1$).

  4. 4.

In the region C ($0<\alpha<1<\beta$), $O_{2}$ is an unstable node, $O_{3}$ is
a saddle, $O_{4}$ has one negative coordinate (at $\alpha\beta\neq 1$).

  5. 5.

In the region D ($\alpha,\beta>1$), the equilibria $O_{2}$ and $O_{3}$ are
unstable nodes, $O_{4}$ is a saddle.

Since the negative values of $\rho_{i}$ are not considered, there are three
equilibrium states in the regions B and C, and they lie on the $\rho_{1}$ and
$\rho_{2}$ axes.

Figure 1 shows the partitioning of the parameter plane into the regions A, B,
C, D and the corresponding phase portraits.

Figure 1: Different types of phase portraits of Eqs. (3).

In all subsequent figures the following notation is used: blue dots correspond
to saddle equilibrium states, green - to stable, red - to unstable.

Using the Bendixon-Dulac theorem, we show that there are no limit cycles in
the phase space of the model (3) in the region $\rho_{1,2}>0$. Indeed,

$\frac{\partial}{\partial\rho_{1}}\frac{\dot{\rho}_{1}}{\rho_{1}\rho_{2}}+\frac{\partial}{\partial\rho_{2}}\frac{\dot{\rho}_{2}}{\rho_{1}\rho_{2}}=\frac{\partial}{\partial\rho_{1}}\frac{\rho_{1}-1+\alpha\rho_{2}}{\rho_{2}}+\frac{\partial}{\partial\rho_{2}}\frac{\rho_{2}-1+\beta\rho_{1}}{\rho_{1}}=\frac{1}{\rho_{2}}+\frac{1}{\rho_{1}}.$

Thus,
$\frac{\partial}{\partial\rho_{1}}\frac{\dot{\rho}_{1}}{\rho_{1}\rho_{2}}+\frac{\partial}{\partial\rho_{2}}\frac{\dot{\rho}_{2}}{\rho_{1}\rho_{2}}>0$
in the region $\rho_{1,2}>0$, this means that there are no limit cycles in
this region.

Let us consider the parameter region $C$ (similarly, the region $B$ is
considered). On the plane ($\rho_{1},\rho_{2}$) there is a trajectory $g_{1}$
going from the unstable node $O_{2}$ to the saddle $O_{3}$. For the three-
dimensional system considered below, this trajectory $g_{1}$ is heteroclinic,
since in a three-dimensional system the equilibrium state $O_{3}$ is a saddle.
Similarly, it is possible to show the presence of heteroclinic trajectories on
the invariant planes $\rho_{1}=0$ and $\rho_{2}=0$. The heteroclinic cycle in
a three-dimensional system consists precisely of such heteroclinic
trajectories.

The system described by Eqs. (3) has one drawback. For the initial conditions
lying below the curve $g_{1}$, all trajectories come to a stable node $O_{1}$.
And for the initial conditions lying above the curve $g_{1}$, all trajectories
go to infinity. This disadvantage disappears in the second model, discussed
below.

####  2.1.1 Bifurcations in a two-dimensional system

The analysis of the eigenvalues of the equilibria given in Table 1 allows us
to determine the bifurcations occurring in Eqs. (3). At $\beta=1$, the
equilibrium states $O_{2}$ and $O_{4}$ undergo transcritical bifurcation. At
$\alpha=1$, the equilibrium states $O_{3}$ and $O_{4}$ also undergo
transcritical bifurcation. At $\alpha\beta=1$, the equilibrium state $O_{4}$
does not exist; when passing through this curve, the sign of one of the
coordinates changes to the opposite. The bifurcation diagram of Eqs. (3) is
shown in Fig. 1.

###  2.2 Dynamics of a three-neuron ensemble (Eqs. (2))

Eqs. (2) can have 8 equilibrium states. The coordinates of the equilibrium
states and their eigenvalues are given in Table 2.

Equilibrium state | Eigenvalues
---|---
$O_{1}(0,0,0)$ | $-1$ | $-1$ | $-1$
$O_{2}(1,0,0)$ | $1$ | $\alpha-1$ | $\beta-1$
$O_{3}(0,1,0)$ | $1$ | $\alpha-1$ | $\beta-1$
$O_{4}(0,0,1)$ | $1$ | $\alpha-1$ | $\beta-1$
$O_{5}(0,\frac{\alpha-1}{\alpha\beta-1},\frac{\beta-1}{\alpha\beta-1})$ | $1$ | $\frac{\alpha^{2}+\beta^{2}-\alpha\beta-\alpha-\beta+1}{\alpha\beta-1}$ | $\frac{(\alpha-1)(1-\beta)}{\alpha\beta-1}$
$O_{6}(\frac{\beta-1}{\alpha\beta-1},0,\frac{\alpha-1}{\alpha\beta-1})$ | $1$ | $\frac{\alpha^{2}+\beta^{2}-\alpha\beta-\alpha-\beta+1}{\alpha\beta-1}$ | $\frac{(\alpha-1)(1-\beta)}{\alpha\beta-1}$
$O_{7}(\frac{\alpha-1}{\alpha\beta-1},\frac{\beta-1}{\alpha\beta-1},0)$ | $1$ | $\frac{\alpha^{2}+\beta^{2}-\alpha\beta-\alpha-\beta+1}{\alpha\beta-1}$ | $\frac{(\alpha-1)(1-\beta)}{\alpha\beta-1}$
$O_{8}(\frac{1}{\alpha+\beta+1},\frac{1}{\alpha+\beta+1},\frac{1}{\alpha+\beta+1})$ | $1$ | $\frac{2-\alpha-\beta+(\alpha-\beta)\sqrt{3}i}{2(\alpha+\beta+1)}$ | $\frac{2-\alpha-\beta-(\alpha-\beta)\sqrt{3}i}{2(\alpha+\beta+1)}$
Table 2: The equilibrium states of Eqs. (2) and their eigenvalues.

Figure 2 shows the bifurcation diagram of Eqs. (2). The following are
descriptions of equilibrium states in all regions of the diagram 2.

  1. 1.

$O_{1}$ is a stable node.

  2. 2.

In the region $A$ ($0<\alpha,\beta<1$) the equilibria $O_{2}$, $O_{3}$ and
$O_{4}$ are saddles with a two-dimensional stable manifold and a one-
dimensional unstable manifold, the equilibria $O_{5}$, $O_{6}$ and $O_{7}$ are
saddles with a one-dimensional stable manifold and a two-dimensional unstable
manifold, the equilibrium state $O_{8}$ is an unstable focus (at
$2-\alpha-\beta<1$) or a node.

  3. 3.

In the regions $B_{1}$ ($\alpha>1$, $\alpha+\beta<2$) and $C_{1}$ ($\beta>1$,
$\alpha+\beta<2$) the equilibria $O_{2}$, $O_{3}$ and $O_{4}$ are saddles with
a two-dimensional unstable manifold and a one-dimensional stable manifold, the
equilibria $O_{5}$, $O_{6}$ and $O_{7}$ are saddles with a one-dimensional
unstable manifold and a two-dimensional stable manifold, the equilibrium state
$O_{8}$ is an unstable focus.

  4. 4.

In the regions $B_{2}$ ($\alpha>1$, $\alpha+\beta>2$, $\alpha\beta<1$) and
$C_{2}$ ($\beta>1$, $\alpha+\beta>2$, $\alpha\beta<1$) the equilibria $O_{2}$,
$O_{3}$ and $O_{4}$ are saddles with a two-dimensional unstable manifold and a
one-dimensional stable manifold, the equilibria $O_{5}$, $O_{6}$ and $O_{7}$
are saddles with a one-dimensional unstable manifold and a two-dimensional
stable manifold, the equilibrium state $O_{8}$ is a saddle-focus.

  5. 5.

In the regions $B_{3}$ ($\beta<1$, $\alpha\beta>1$) and $C_{3}$ ($\alpha<1$,
$\alpha\beta>1$) the equilibria $O_{2}$, $O_{3}$ and $O_{4}$ are saddles with
a two-dimensional unstable manifold and a one-dimensional stable manifold, the
equilibria $O_{5}$, $O_{6}$ and $O_{7}$ are unstable nodes, the equilibrium
state $O_{8}$ is a saddle-focus.

  6. 6.

In the region $D$ ($\alpha,\beta>1$) the equilibria $O_{2}$, $O_{3}$ and
$O_{4}$ are unstable nodes, the equilibria $O_{5}$, $O_{6}$ and $O_{7}$ are
saddles with a one-dimensional stable manifold and a two-dimensional unstable
manifold, the equilibrium state $O_{8}$ is a saddle-focus.

When $\alpha\beta=1$, the equilibria $O_{5}$, $O_{6}$, $O_{7}$ do not exist.
In addition, at $\alpha+\beta=2$ the real part of the two eigenvalues of the
equilibrium state $O_{8}$ turns into $0$. At $\alpha+\beta<2$, this
equilibrium state is an unstable saddle-focus with three unstable manifolds,
at $\alpha+\beta>2$, it is a saddle-focus with one-dimensional unstable and
two-dimensional stable manifolds. Let’s examine this bifurcation in detail.

When the condition is met $\alpha+\beta=2$ we have

$\dot{\rho}_{1}+\dot{\rho}_{2}+\dot{\rho}_{3}=\rho_{1}(\rho_{1}-1+\alpha\rho_{2}+(2-\alpha)\rho_{3})+\rho_{2}(\rho_{2}-1+\alpha\rho_{3}+(2-\alpha)\rho_{1})+\rho_{3}(\rho_{3}-1+\alpha\rho_{1}+(2-\alpha)\rho_{2})=$
$=\rho_{1}^{2}+\rho_{2}^{2}+\rho_{3}^{2}-\rho_{1}-\rho_{2}-\rho_{3}+2(\rho_{1}\rho_{3}+\rho_{2}\rho_{1}+\rho_{3}\rho_{2})=(\rho_{1}+\rho_{2}+\rho_{3})^{2}-\rho_{1}-\rho_{2}-\rho_{3}.$

If $\rho_{1}+\rho_{2}+\rho_{3}=1$, we get the identity $0=0$. Therefore, the
plane $\rho_{1}+\rho_{2}+\rho_{3}=1$ is an invariant set.

On the invariant plane $\rho_{1}+\rho_{2}+\rho_{3}=1$, provided
$\alpha+\beta=2$, Eqs. (2) takes the form

$\begin{cases}\dot{\rho}_{1}=(\alpha-1)\rho_{1}(\rho_{1}-1+2\rho_{2})\\\
\dot{\rho}_{2}=(1-\alpha)\rho_{2}(\rho_{2}-1+2\rho_{1})\end{cases}.$ (4)

This system can be transformed by substitution of time into

$\begin{cases}\dot{\rho}_{1}=-\rho_{1}(\rho_{1}-1+2\rho_{2})\\\
\dot{\rho}_{2}=\rho_{2}(\rho_{2}-1+2\rho_{1})\end{cases}.$ (5)

The solution of this system is

$2\rho_{1}\rho_{2}^{2}-2\rho_{1}\rho_{2}(1-\rho_{1})-C=0.$

The diagrams of this function inside the region $\rho_{1,2}>0$,
$\rho_{1}+\rho_{2}<1$ are closed curves. Thus, under the condition of
$\alpha+\beta=2$ on the invariant plane $\rho_{1}+\rho_{2}+\rho_{3}=1$ through
each point of the set $\rho_{1,2,3}>0$ a closed curve passes.

Figure 3 shows all the different types of phase portraits of Eqs. (2) for
parameters from the regions shown in Fig. 2.

Figure 2: Bifurcation diagram of Eqs. (2).

(a) The region $A$ (b) The region $C_{1}$

(c) The regions $C_{2}$ and $C_{3}$ (d) The region $D$

Figure 3: Various types of phase portraits of Eqs. (2). The heteroclinic
trajectories forming the heteroclinic cycle $G$ are shown in red in Figs. b)
and c) the invariant manifold $M$ is highlighted in gray. In the region
$C_{1}$, the manifold $M$ is convex, and in the regions $C_{2}$ and $C_{3}$ is
concave.

The dynamics of an ensemble of three neurons, depending on the parameters, can
be as follows:

  1. 1.

For all regions (Fig. 3), depending on the external stimulus, the phase
trajectories either come to a stable equilibrium state $O_{1}$ or go to
infinity.

  2. 2.

For the regions $C_{1}$, $C_{2}$ and $C_{3}$, the boundary between the basins
of attraction of the equilibrium and infinity is the integral manifold $M$
containing the saddle points $O_{2}$, $O_{3}$ and $O_{4}$ and the saddle-focus
$O_{8}$.

In the previous section, it was shown that in the regions $C_{1}$, $C_{2}$ and
$C_{3}$ in the non-negative quadrant of the invariant plane $\rho_{3}=0$ there
are three equilibrium states: a stable node $O_{1}(0,0)$, an unstable node
$O_{3}(0,1)$ and a saddle $O_{2}(1,0)$. When adding the third coordinate
$\rho_{3}$, the equilibrium state $O_{1}(0,0,0)$ is a stable node, the
equilibrium states $O_{3}(0,1,0)$ and $O_{2}(1,0,0)$ are saddles with a two-
dimensional unstable manifold and a one-dimensional stable manifold. Between
the saddle equilibrium states $O_{2}(1,0,0)$ and $O_{3}(0,1,0)$ there is a
heteroclinic trajectory $g_{1}$. Due to the symmetry of Eqs. (2), between
pairs of saddles $O_{3}(0,1,0)$, $O_{4}(0,0,1)$ and $O_{4}(0,0,1)$,
$O_{2}(1,0,0)$ there are also heteroclinic trajectories $g_{2}$ and $g_{3}$,
respectively. Thus, in the phase space of Eqs. (2) there is a heteroclinic
cycle $G$ connecting the saddles $O_{2}$, $O_{3}$ and $O_{4}$. Like the usual
limit cycle, the heteroclinic cycle can be stable, unstable and saddle. This
is determined by the stability (instability) of the manifolds forming this
cycle at the intersection. The heteroclinic cycle $G$ unstable at
$\alpha+\beta<2$ (manifolds $\rho_{1}$, $\rho_{2}$, $\rho_{3}$ and $M$ are
unstable) and saddle at $\alpha+\beta>2$ (the manifold $M$ is stable).

####  2.2.1 Sequential activity

As mentioned above, one of the behaviors of a neural ensemble is the alternate
activation of neurons as a result of a response to an external stimulus,
followed by the return of all neurons to their original non-excited state.
Mathematical images of such behavior can be phase trajectories passing near
the heteroclinic cycle $G$. In the model under consideration, this cycle is a
closed trajectory consisting of heteroclinic curves connecting the saddle
equilibrium states $O_{2}$, $O_{3}$, $O_{4}$, i.e. the saddles lying on the
axes. Only in the regions $B_{i}$ and $C_{i}$ these saddles have one-
dimensional unstable manifolds that connect them. Thus, the sequential
activity of neurons can be realized only in the case of asymmetric
connections. Note that the sequential neural activity was also observed in
ensembles with inhibitory asymmetric connections [19, 20, 21].

Further, Figs. 4-10 show phase portraits and alternating realizations of
several possible behaviors of a neural ensemble depending on the external
stimulus determined by the values $a$, $b$ and $c$. We will assume that the
$i$-th neuron is active if the value of $\rho_{i}$ exceeds the value of 0.5.
Initially, all neurons are inactive. The image point is in the equilibrium
state $O_{1}$. In order to activate the neuron, you need to send such an
external stimulus to transfer the initial phase point from the equilibrium
state $O_{1}$ to a small neighborhood of an unstable or saddle heteroclinic
cycle $G$ or to a small neighborhood of an invariant manifold $M$ in the case
of a saddle heteroclinic cycle $G$. Depending on the external stimulus, the
following dynamics can be obtained:

  * •

excitation of a single neuron (Fig. 4)

  * •

sequential excitation of two neurons (Fig. 5)

  * •

sequential one-time excitation of all three neurons (Figs. 6, 7 and 10)

  * •

sequential multiple excitation of all three neurons (Fig. 8)

  * •

phase trajectory goes to infinity (Fig. 9)

Figure 4: The process of activation and damping. $a$, $b$, $c$ parameters are
selected so that only one element is activated - the second one. (a) Phase
portrait. (b) Time diagram. Here and further, the red curve corresponds to
$\rho_{1}(t)$, the green curve is $\rho_{2}(t)$, the blue curve is
$\rho_{3}(t)$. Parameter values: $\alpha=0.1$, $\beta=1.5$, $a=0.6917$,
$b=0.4$, $c=0.01$.

Figure 5: Activation, sequential activity and damping. $a$, $b$, $c$
parameters are selected so that only two elements are activated: the second
and the third ones. (a) Phase portrait. (b) Time diagram. Parameter values:
$\alpha=0.1$, $\beta=1.5$, $a=0.69173179$, $b=0.4$, $c=0.01$.

Figure 6: Activation, sequential activity and damping. All three neurons are
activated sequentially once. (a) Phase portrait: the phase trajectory takes
place in a vicinity of the heteroclinic cycle due to the selection of an
external stimuli. (b) Time diagram. Parameter values: $\alpha=0.1$,
$\beta=1.5$, $a=0.691731794346462$, $b=0.4$, $c=0.01$.

Figure 7: Activation, sequential activity and damping. All three neurons are
activated sequentially once. (a) Phase portrait: the phase trajectory does not
take place in a vicinity of the heteroclinic cycle due to the selection of an
external stimuli. (b) Time diagram. Parameter values: $\alpha=0.1$,
$\beta=1.5$, $\alpha+\beta<2$, $a=0.62313692$, $b=0.4$, $c=0.1$.

Figure 8: Activation, sequential activity and damping. The second and third
neurons are activated twice. The first neuron is activated once. (a) Phase
portrait: the phase trajectory takes place in a vicinity of the heteroclinic
cycle due to the selection of an external stimuli moving the starting point to
the integral manifold $M$. (b) Time diagram. Parameter values: $\alpha=0.1$,
$\beta=1.5$, $a=0.540537789125142$, $b=0.4$, $c=0.2$.

Figure 9: Sequential activity. The phase trajectory goes to infinity. (a)
Phase portrait. (b) Time diagram. Parameter values: $\alpha=0.1$, $\beta=1.5$,
$a=0.691731794346463$, $b=0.4$, $c=0.01$.

Figure 10: Activation, sequential activity and damping.(a) Phase portrait. (b)
Time diagram. Parameter values: $\alpha=0.4$, $\beta=1.9$, $\alpha+\beta>2$,
$a=0.52382217$, $b=0.4$, $c=0.01$.

##  3 The second type system

To eliminate the disadvantage of the first system, the possibility of
trajectories going to infinity, in this section we will analyze the dynamics
of the ensemble modeled by the following equation

$\dot{\rho}=\rho\left(\gamma-\rho\right)(\rho-1).$ (6)

The $\gamma$ parameter is the threshold value of the ensemble’s activity. It’s
value satisfies the inequalities $0<\gamma<1$.

Eqs. (6) has three equilibrium states: $\rho=0$ and $\rho=1$ are stable and
$\rho=\gamma$ is unstable. If the initial condition is $\rho_{0}<\gamma$, then
the activity of the $\rho$ ensemble decreases and tends to zero. If
$\rho_{0}>\gamma$, then the activity of the $\rho$ ensemble increases and
tends to one. It follows from the above that the set $0\leq\rho\leq 1$ is
invariant.

The model of an ensemble of three connected elements modeled by Eq. (6) has
the form

$\begin{cases}\dot{\rho_{1}}=\rho_{1}(\gamma-\rho_{1})(\rho_{1}-1+\alpha\rho_{2}+\beta\rho_{3})+a\,\delta(\rho_{1})\\\
\dot{\rho_{2}}=\rho_{2}(\gamma-\rho_{2})(\rho_{2}-1+\alpha\rho_{3}+\beta\rho_{1})+b\,\delta(\rho_{2})\\\
\dot{\rho_{3}}=\rho_{3}(\gamma-\rho_{3})(\rho_{3}-1+\alpha\rho_{1}+\beta\rho_{2})+c\,\delta(\rho_{3})\end{cases}.$
(7)

###  3.1 Study of a two-dimensional system on an invariant plane $\rho_{3}=0$

Since when one of the conditions $\rho_{i}=0$ is met, the i-th equation of
Eqs. (7) turns into an identity, then the three planes $\rho_{i}=0$ are
invariant.

On the invariant plane $\rho_{3}=0$, Eqs. (7) takes the form

$\begin{cases}\dot{\rho}_{1}=\rho_{1}(\gamma-\rho_{1})(\rho_{1}-1+\alpha\rho_{2})+a\,\delta(\rho_{1})\\\
\dot{\rho}_{2}=\rho_{2}(\gamma-\rho_{2})(\rho_{2}-1+\beta\rho_{1})+b\,\delta(\rho_{2})\end{cases}.$
(8)

The equilibrium states of Eqs. (8) and their eigenvalues are given in Table 3:

Equilibrium state | Eigenvalues
---|---
$O_{1}(0,0)$ | $-\gamma$ | $-\gamma$
$O_{2}(1,0)$ | $\gamma(\beta-1)$ | $\gamma-1$
$O_{3}(0,1)$ | $\gamma(\alpha-1)$ | $\gamma-1$
$O_{4}(1-\alpha\gamma,\gamma)$ | $\gamma(\beta(\alpha\gamma-1)-\gamma+1)$ | $(1-\alpha\gamma)(\gamma+\alpha\gamma-1)$
$O_{5}(\gamma,1-\beta\gamma)$ | $\gamma(\alpha(\beta\gamma-1)-\gamma+1)$ | $(1-\beta\gamma)(\gamma+\beta\gamma-1)$
$O_{6}(\gamma,0)$ | $\gamma(\beta\gamma-1)$ | $\gamma(1-\gamma)$
$O_{7}(0,\gamma)$ | $\gamma(\alpha\gamma-1)$ | $\gamma(1-\gamma)$
$O_{8}(\gamma,\gamma)$ | $\gamma(1-\gamma-\alpha\gamma)$ | $\gamma(1-\gamma-\beta\gamma)$
$O_{9}(\frac{\alpha-1}{\alpha\beta-1},\frac{\beta-1}{\alpha\beta-1})$ | $\lambda_{1}$ | $\lambda_{2}$
Table 3: The equilibrium states of Eqs. (8) and their eigenvalues.

The eigenvalues $\lambda_{1,2}$ of the equilibrium state $O_{9}$ are too
cumbersome and therefore are not given in the table. They are found from the
characteristic equation

${\textstyle\lambda^{2}+\frac{(\alpha-1)(\alpha+\gamma-\alpha\beta\gamma-1)+(\beta-1)(\beta+\gamma-\alpha\beta\gamma-1)}{(\alpha\beta-1)^{2}}\lambda-\frac{(\alpha-1)(\alpha+\gamma-\alpha\beta\gamma-1)(\beta-1)(\beta+\gamma-\alpha\beta\gamma-1)}{(\alpha\beta-1)^{3}}=0.}$

$O_{9}$ undergoes Andronov-Hopf bifurcation under the following conditions.

$\begin{cases}\alpha\beta\gamma(\alpha+\beta-2)-\alpha^{2}-\beta^{2}+(2-\gamma)(\alpha+\beta)+2\gamma-2=0\\\
\frac{\left(\alpha-1\right)\left(\beta-1\right)\left(\alpha\beta\gamma(\alpha\beta\gamma-\alpha-\beta-2\gamma+2)+(\gamma-1)(\alpha+\beta)+\alpha\beta+\gamma^{2}-2\gamma+1\right)}{(1-\alpha\beta)^{3}}>0\end{cases}.$
(9)

$\beta$ parameter can be expressed from the first equation:

$\beta=\frac{-\alpha^{2}\gamma+2\alpha\gamma+\gamma-2\pm(\alpha-1)\sqrt{\alpha^{2}\gamma^{2}+2\alpha\gamma(2-\gamma)+\gamma^{2}+4\gamma-4}}{2(\alpha\gamma-1)}.$

At the intersection of this curve in the direction of decreasing $\beta$
parameter, a stable limit cycle is born in the phase space, which further
increases in size and gets into the heteroclinic contour formed by the saddles
$O_{3}$, $O_{4}$ and $O_{7}$ (Fig. 11).

(a) $\beta=1.3055$ (b) $\beta=1.324802754125605$

(c) $\beta=1.3325$ (d) $\beta=1.34$

Figure 11: Andronov–Hopf bifurcation and heteroclinic bifurcation. In Figure
(b), the heteroclinic cycle connecting the saddles $O_{3}$, $O_{4}$ and
$O_{7}$ is highlighted in red. In Figure (c), a stable limit cycle is
highlighted in red. In Figure (d), the limit cycle does not exist, while
$O_{9}$ is stable. Parameter values: $\alpha=1.1$, $\gamma=0.6$.

Similarly with the study of the first type system on invariant planes
$\rho_{i}=0$, we determine the types of equilibrium states of Eqs. (8)
depending on the values of the $\alpha$ and $\beta$ parameters. The number of
regions of the parameter space with different number and type of equilibrium
states in this model is greater. It is enough to describe only those regions
of the parameter space that lie either above the diagonal $\alpha=\beta$, or
intersect it, since in the regions symmetrical with respect to the diagonal,
the phase portraits differ only in the numbering of the variables $\rho_{1}$
and $\rho_{2}$.

  1. 1.

$O_{1}$ is always a stable node.

  2. 2.

$O_{2}$ is a stable node in the regions $A$, $B$, $C$, $H$, $I$ ($0<\beta<1$)
and a saddle in the remaining regions.

  3. 3.

$O_{3}$ is a stable node in the regions $A$, $B$, $C$, $D$, $E$, $F$, $G$,
$H$, $I$, $J$, $K$, $L$ ($0<\alpha<1$) and a saddle in the remaining regions.

  4. 4.

$O_{4}$ is a stable node in the region $A$ and has one negative coordinate in
the region $S$, a saddle in the regions $B$, $C$, $D$, $E$, $F$, $G$, $Q$,
$R$, and an unstable node in the remaining regions.

  5. 5.

$O_{5}$ is stable node in the regions $A$, $B$ and has one negative coordinate
in the regions $F$, $G$, $L$, $N$, $P$, $R$, $S$, an unstable node in the
region $H$ and a saddle in the remaining regions.

  6. 6.

$O_{6}$ is an unstable node in the regions $F$, $G$, $L$, $N$, $P$, $R$, $S$
($\beta>\frac{1}{\gamma}$) and a saddle in the remaining regions.

  7. 7.

$O_{7}$ is an unstable node in the region $S$ ($\alpha>\frac{1}{\gamma}$) and
a saddle in the remaining regions.

  8. 8.

$O_{8}$ is an unstable node in the regions $A$ and $B$, a saddle in the
regions $C$, $D$, $E$, $F$ and $G$ and a stable node in the remaining regions.

  9. 9.

$O_{9}$ is stable node in the region $A$, a saddle in the regions $B$, $C$,
$I$, $Q$, $R$ and $S$, an unstable node in the region $H$, a stable focus in
the regions $M$ and $N$, an unstable focus in the regions $O$ and $P$ and has
one negative coordinate in the remaining regions.

Table 4 provides a summary list of the types of all nine equilibria and the
regions in which they exist.

Region | $O_{1}$ | $O_{2}$ | $O_{3}$ | $O_{4}$ | $O_{5}$ | $O_{6}$ | $O_{7}$ | $O_{8}$ | $O_{9}$
---|---|---|---|---|---|---|---|---|---
$A$ | SN | SN | SN | SN | SN | S | S | UN | SN
$B$ | SN | SN | SN | S | SN | S | S | UN | S
$C$ | SN | SN | SN | S | S | S | S | S | S
$D$ | SN | S | SN | S | S | S | S | S | -
$E$ | SN | S | SN | S | S | S | S | S | -
$F$ | SN | S | SN | S | - | UN | S | S | -
$G$ | SN | S | SN | S | - | UN | S | S | -
$H$ | SN | SN | SN | UN | UN | S | S | SN | UN
$I$ | SN | SN | SN | UN | S | S | S | SN | S
$J$ | SN | S | SN | UN | S | S | S | SN | -
$K$ | SN | S | SN | UN | S | S | S | SN | -
$L$ | SN | S | SN | UN | - | UN | S | SN | -
$M$ | SN | S | S | UN | S | S | S | SN | SF
$N$ | SN | S | S | UN | - | UN | S | SN | SF
$O$ | SN | S | S | UN | S | S | S | SN | UF
$P$ | SN | S | S | UN | - | UN | S | SN | UF
$Q$ | SN | S | S | S | S | S | S | SN | S
$R$ | SN | S | S | S | - | UN | S | SN | S
$S$ | SN | S | S | - | - | UN | UN | SN | S
Table 4: Types of equilibrium states in all regions of the parameter space
(the partitioning into regions is shown in Fig. 12). A dash stands for a
presence of a negative coordinate of the equilibrium. SN - stable node, SF -
stable focus, UN - unstable node, UF - unstable focus, S - saddle.

In the system under consideration, two bifurcations of equilibrium states are
possible: the Andronov-Hopf bifurcation described above and the transcritical
bifurcation. Table 5 shows a list of transcritical bifurcations of the
equilibrium states of the system.

| $O_{1}$ | $O_{2}$ | $O_{3}$ | $O_{4}$ | $O_{5}$ | $O_{6}$ | $O_{7}$ | $O_{8}$ | $O_{9}$
---|---|---|---|---|---|---|---|---|---
$O_{1}(0,0)$ | - | - | - | - | - | - | - | - | -
$O_{2}(1,0)$ | - | - | - | - | - | - | - | - | $\beta=1$
$O_{3}(0,1)$ | - | - | - | - | - | - | - | - | $\alpha=1$
$O_{4}(1-\alpha\gamma,\gamma)$ | - | - | - | - | - | - | $\alpha=\frac{1}{\gamma}$ | $\alpha=\frac{1}{\gamma}-1$ | $\beta=\frac{\gamma-1}{\alpha\gamma-1}$
$O_{5}(\gamma,1-\beta\gamma)$ | - | - | - | - | - | $\beta=\frac{1}{\gamma}$ | - | $\beta=\frac{1}{\gamma}-1$ | $\beta=\frac{\alpha+\gamma-1}{\alpha\gamma}$
$O_{6}(\gamma,0)$ | - | - | - | - | $\beta=\frac{1}{\gamma}$ | - | - | - | -
$O_{7}(0,\gamma)$ | - | - | - | $\alpha=\frac{1}{\gamma}$ | - | - | - | - | -
$O_{8}(\gamma,\gamma)$ | - | - | - | $\alpha=\frac{1}{\gamma}-1$ | $\beta=\frac{1}{\gamma}-1$ | - | - | - | -
$O_{9}(\frac{\alpha-1}{\alpha\beta-1},\frac{\beta-1}{\alpha\beta-1})$ | - | $\beta=1$ | $\alpha=1$ | $\beta=\frac{\gamma-1}{\alpha\gamma-1}$ | $\beta=\frac{\alpha+\gamma-1}{\alpha\gamma}$ | - | - | - | -
Table 5: Lines on which transcritical bifurcations of equilibrium states
occur.

In Figure 12 the bifurcation diagram of the system described by Eqs. (8) is
given for $\gamma=0.6$.

Figure 12: The bifurcation diagram of Eqs. (8) at $\gamma=0.6$. The black
lines correspond to transcritical bifurcations, the red line corresponds to
bifurcations of (9).

Phase portraits of Eqs. (7) for all marked areas are shown in Fig. 13.
Similarly with the first system, the study of two-dimensional Eqs. (8) should
provide information about the dynamics of a three-dimensional system, in which
the main attention will be paid again to the existence of a heteroclinic
cycle. The trajectories of two-dimensional systems on invariant planes
$\rho_{i}=0$ will be components of a three-dimensional heteroclinic cycle. As
for the first system on invariant planes $\rho_{i}=0$, closed trajectories
including an unstable node, saddle equilibrium states and a stable node at the
origin are of interest. Moreover, the unstable node and one of two saddles
should lie on different coordinate axes. Among all other regions on the
$\alpha\beta-plane$ $F$, $G$ and $L$ meet all the requirements. In these
regions, there are trajectories connecting unstable nodes, saddle equilibrium
states and a stable node $O_{1}$. For example, in Fig. 13 (e), for the regions
$F$ and $G$, there is a contour containing the equilibrium states $O_{6}$,
$O_{8}$, $O_{7}$ and $O_{1}$. This contour will be an integral part of the
heteroclinic cycle in the three-dimensional phase space of the initial system.

(a) Region A (b) Region B (c) Region C (d) Regions D and E

(e) Regions F and G (f) Region H (g) Region I (h) Regions J and K

(i) Region L (j) Region M (k) Region N (l) Region O

(m) Region P (n) Region Q (o) Region R (p) Region S

Figure 13: Phase portraits in various regions of Fig. 12.

###  3.2 A study of Eqs. (7)

The equilibria of Eqs. (7) and their eigenvalues are given in the Table 6.

№ | Equilibrium state | Eigenvalues
---|---|---
$1$ | $(0,0,0)$ | $-\gamma$ | $-\gamma$ | $-\gamma$
$2$ | $(1,0,0)$ | $\gamma-1$ | $\gamma(\alpha-1)$ | $\gamma(\beta-1)$
$3$ | $(0,1,0)$ | $\gamma-1$ | $\gamma(\alpha-1)$ | $\gamma(\beta-1)$
$4$ | $(0,0,1)$ | $\gamma-1$ | $\gamma(\alpha-1)$ | $\gamma(\beta-1)$
$5$ | $(\gamma,0,0)$ | $\gamma(\alpha\gamma-1)$ | $\gamma(\beta\gamma-1)$ | $\gamma(1-\gamma)$
$6$ | $(0,\gamma,0)$ | $\gamma(\alpha\gamma-1)$ | $\gamma(\beta\gamma-1)$ | $\gamma(1-\gamma)$
$7$ | $(0,0,\gamma)$ | $\gamma(\alpha\gamma-1)$ | $\gamma(\beta\gamma-1)$ | $\gamma(1-\gamma)$
$8$ | $(\gamma,\gamma,0)$ | $\gamma(\alpha\gamma+\beta\gamma-1)$ | $\gamma(1-\gamma-\alpha\gamma)$ | $\gamma(1-\gamma-\beta\gamma)$
$9$ | $(\gamma,0,\gamma)$ | $\gamma(\alpha\gamma+\beta\gamma-1)$ | $\gamma(1-\gamma-\alpha\gamma)$ | $\gamma(1-\gamma-\beta\gamma)$
$10$ | $(0,\gamma,\gamma)$ | $\gamma(\alpha\gamma+\beta\gamma-1)$ | $\gamma(1-\gamma-\alpha\gamma)$ | $\gamma(1-\gamma-\beta\gamma)$
$11$ | $(\gamma,\gamma,\gamma)$ | $\gamma(1-\gamma-\alpha\gamma-\beta\gamma)$ | $\gamma(1-\gamma-\alpha\gamma-\beta\gamma)$ | $\gamma(1-\gamma-\alpha\gamma-\beta\gamma)$
$12$ | $(1-\alpha\gamma,\gamma,0)$ | $\gamma(\alpha+\beta\gamma-\alpha^{2}\gamma-1)$ | $\gamma(\alpha\beta\gamma-\beta-\gamma+1)$ | $-(\gamma-\frac{1}{\alpha})(\gamma-\frac{1}{\alpha+1})$
$13$ | $(\gamma,0,1-\alpha\gamma)$ | $\gamma(\alpha+\beta\gamma-\alpha^{2}\gamma-1)$ | $\gamma(\alpha\beta\gamma-\beta-\gamma+1)$ | $-(\gamma-\frac{1}{\alpha})(\gamma-\frac{1}{\alpha+1})$
$14$ | $(0,1-\alpha\gamma,\gamma)$ | $\gamma(\alpha+\beta\gamma-\alpha^{2}\gamma-1)$ | $\gamma(\alpha\beta\gamma-\beta-\gamma+1)$ | $-(\gamma-\frac{1}{\alpha})(\gamma-\frac{1}{\alpha+1})$
$15$ | $(\gamma,1-\beta\gamma,0)$ | $\gamma(\beta+\alpha\gamma-\beta^{2}\gamma-1)$ | $\gamma(\alpha\beta\gamma-\alpha-\gamma+1)$ | $-(\gamma-\frac{1}{\beta})(\gamma-\frac{1}{\beta+1})$
$16$ | $(1-\beta\gamma,0,\gamma)$ | $\gamma(\beta+\alpha\gamma-\beta^{2}\gamma-1)$ | $\gamma(\alpha\beta\gamma-\alpha-\gamma+1)$ | $-(\gamma-\frac{1}{\beta})(\gamma-\frac{1}{\beta+1})$
$17$ | $(0,\gamma,1-\beta\gamma)$ | $\gamma(\beta+\alpha\gamma-\beta^{2}\gamma-1)$ | $\gamma(\alpha\beta\gamma-\alpha-\gamma+1)$ | $-(\gamma-\frac{1}{\beta})(\gamma-\frac{1}{\beta+1})$
$18$ | $(\frac{\alpha+\beta\gamma-\alpha^{2}\gamma-1}{\alpha\beta-1},\frac{\beta+\alpha\gamma-\beta^{2}\gamma-1}{\alpha\beta-1},\gamma)$ | $\lambda_{1}$ | $\lambda_{2}$ | $\lambda_{3}$
$19$ | $(\frac{\beta+\alpha\gamma-\beta^{2}\gamma-1}{\alpha\beta-1},\gamma,\frac{\alpha+\beta\gamma-\alpha^{2}\gamma-1}{\alpha\beta-1})$ | $\lambda_{1}$ | $\lambda_{2}$ | $\lambda_{3}$
$20$ | $(\gamma,\frac{\alpha+\beta\gamma-\alpha^{2}\gamma-1}{\alpha\beta-1},\frac{\beta+\alpha\gamma-\beta^{2}\gamma-1}{\alpha\beta-1})$ | $\lambda_{1}$ | $\lambda_{2}$ | $\lambda_{3}$
$21$ | $(\frac{\alpha-1}{\alpha\beta-1},\frac{\beta-1}{\alpha\beta-1},0)$ | $\mu_{1}$ | $\mu_{2}$ | $\mu_{3}$
$22$ | $(\frac{\beta-1}{\alpha\beta-1},0,\frac{\alpha-1}{\alpha\beta-1})$ | $\mu_{1}$ | $\mu_{2}$ | $\mu_{3}$
$23$ | $(0,\frac{\alpha-1}{\alpha\beta-1},\frac{\beta-1}{\alpha\beta-1})$ | $\mu_{1}$ | $\mu_{2}$ | $\mu_{3}$
$24$ | $(1-\gamma(\alpha+\beta),\gamma,\gamma)$ | $\nu_{1}$ | $\nu_{2}$ | $\nu_{3}$
$25$ | $(\gamma,1-\gamma(\alpha+\beta),\gamma)$ | $\nu_{1}$ | $\nu_{2}$ | $\nu_{3}$
$26$ | $(\gamma,\gamma,1-\gamma(\alpha+\beta))$ | $\nu_{1}$ | $\nu_{2}$ | $\nu_{3}$
$27$ | $(\frac{1}{\alpha+\beta+1},\frac{1}{\alpha+\beta+1},\frac{1}{\alpha+\beta+1})$ | $\frac{\gamma(1+\alpha+\beta)-1}{\alpha+\beta+1}$ | ${\scriptstyle\frac{(\gamma(1+\alpha+\beta)-1)(\alpha+\beta-2+(\beta-\alpha)\sqrt{3}i)}{2(\alpha+\beta+1)^{2}}}$ | ${\scriptstyle\frac{(\gamma(1+\alpha+\beta)-1)(\alpha+\beta-2-(\beta-\alpha)\sqrt{3}i)}{2(\alpha+\beta+1)^{2}}}$
Table 6: Equilibria of Eqs. (7).

Expressions for the eigenvalues $\lambda_{1,2,3}$, as well as the
characteristic equation from which they are determined, are too cumbersome.
The eigenvalues $\mu_{1,2,3}$ are determined from the characteristic equation

$\scriptstyle\left(\mu^{2}+\frac{(\alpha-1)(\alpha+\gamma-\alpha\beta\gamma-1)+(\beta-1)(\beta+\gamma-\alpha\beta\gamma-1)}{(\alpha\beta-1)^{2}}\mu-\frac{(\alpha-1)(\alpha+\gamma-\alpha\beta\gamma-1)(\beta-1)(\beta+\gamma-\alpha\beta\gamma-1)}{(\alpha\beta-1)^{3}}\right)\left(\mu-\frac{\alpha^{2}+\beta^{2}+\alpha\beta+\alpha+\beta-1}{\alpha\beta-1}\right)=0.$

The eigenvalues $\nu_{1,2,3}$ are determined from the characteristic equation

${\scriptstyle\left(\nu+(\alpha\gamma+\beta\gamma-1)(\alpha\gamma+\beta\gamma+\gamma-1)\right)\left(\nu^{2}-\gamma(\alpha\gamma+\beta\gamma+\gamma-1)(\alpha+\beta-2)\nu+\gamma^{2}(\alpha\gamma+\beta\gamma+\gamma-1)^{2}(\alpha-1)(\beta-1)\right)=0.}$

The subset $0\leq\rho_{1},\rho_{2},\rho_{3}\leq\gamma$ of the phase space of
Eqs. (7) is invariant. In order to simplify the study of the system, we impose
conditions on $\alpha,\beta$ parameters under which the equilibria
$O_{12}$-$O_{26}$ are not in this subset.

$O_{12}$-$O_{14}$ equilibria have a coordinate $1-\alpha\gamma$ that goes
beyond at $\alpha>\frac{1}{\gamma}$ (region S) and at
$\alpha<\frac{1}{\gamma}-1$ (regions $A$-$L$).

$O_{15}$-$O_{17}$ equilibria have a coordinate $1-\beta\gamma$ that goes
beyond at $\beta>\frac{1}{\gamma}$ (regions $F$, $G$, $L$, $N$, $P$, $R$ and
$S$) and at $\beta<\frac{1}{\gamma}-1$ (regions $A$ and $B$).

$O_{18}$-$O_{20}$ equilibria have a coordinate
$\frac{\alpha+\beta\gamma-\alpha^{2}\gamma-1}{\alpha\beta-1}$ that goes beyond
in the regions $E$, $F$, $G$, $J$-$N$, $P$ and $R$.

$O_{21}$-$O_{23}$ equilibria have a coordinates
$\frac{\alpha-1}{\alpha\beta-1}$ and $\frac{\beta-1}{\alpha\beta-1}$ that go
beyond in the regions $A$-$G$ and $I$-$P$.

$O_{24}$-$O_{26}$ equilibria have a coordinate $1-\gamma(\alpha+\beta)$ that
goes beyond at $\beta>\frac{1}{\gamma}-\alpha$ (regions $F$, $G$, $J$-$S$) and
at $\beta<\frac{1}{\gamma}-\alpha-1$ that is never met in the current regions.

Equilibrium state | Eigenvalues
---|---
$O_{1}(0,0,0)$ | $-\gamma$ | $-\gamma$ | $-\gamma$
$O_{5}(\gamma,0,0)$ | $\gamma(\alpha\gamma-1)$ | $\gamma(\beta\gamma-1)$ | $\gamma(1-\gamma)$
$O_{6}(0,\gamma,0)$ | $\gamma(\alpha\gamma-1)$ | $\gamma(\beta\gamma-1)$ | $\gamma(1-\gamma)$
$O_{7}(0,0,\gamma)$ | $\gamma(\alpha\gamma-1)$ | $\gamma(\beta\gamma-1)$ | $\gamma(1-\gamma)$
$O_{8}(\gamma,\gamma,0)$ | $\gamma(\alpha\gamma+\beta\gamma-1)$ | $\gamma(1-\gamma-\alpha\gamma)$ | $\gamma(1-\gamma-\beta\gamma)$
$O_{9}(\gamma,0,\gamma)$ | $\gamma(\alpha\gamma+\beta\gamma-1)$ | $\gamma(1-\gamma-\alpha\gamma)$ | $\gamma(1-\gamma-\beta\gamma)$
$O_{10}(0,\gamma,\gamma)$ | $\gamma(\alpha\gamma+\beta\gamma-1)$ | $\gamma(1-\gamma-\alpha\gamma)$ | $\gamma(1-\gamma-\beta\gamma)$
$O_{11}(\gamma,\gamma,\gamma)$ | $\gamma(1-\gamma-\alpha\gamma-\beta\gamma)$ | $\gamma(1-\gamma-\alpha\gamma-\beta\gamma)$ | $\gamma(1-\gamma-\alpha\gamma-\beta\gamma)$
$O_{27}(\frac{1}{\alpha+\beta+1},\frac{1}{\alpha+\beta+1},\frac{1}{\alpha+\beta+1})$ | $\frac{\gamma(1+\alpha+\beta)-1}{\alpha+\beta+1}$ | ${\scriptstyle\frac{(\gamma(1+\alpha+\beta)-1)(\alpha+\beta-2+(\beta-\alpha)\sqrt{3}i)}{2(\alpha+\beta+1)^{2}}}$ | ${\scriptstyle\frac{(\gamma(1+\alpha+\beta)-1)(\alpha+\beta-2-(\beta-\alpha)\sqrt{3}i)}{2(\alpha+\beta+1)^{2}}}$
Table 7: Equilibria of Eqs. (7) located in the invariant set
$0\leq\rho_{1},\rho_{2},\rho_{3}\leq 1$.

Table 7 contains a list of equilibria and the corresponding eigenvalues
located in the invariant set $0\leq\rho_{1},\rho_{2},\rho_{3}\leq 1$. It is
these equilibrium states that will be considered in the study of switching
activity associated with the presence of heteroclinic cycles $G$. Unstable
heteroclinic cycles will contain equilibria located on the planes
$\rho_{i}=0$. The attractors are the equilibrium states $O_{1}$ – all neurons
are deactivated and $O_{27}$ – all neurons are active. Depending on the
amplitude of the external stimuli, the neural ensemble, which was initially in
an unexcited state, will resume its initial state or all the neurons of the
ensemble will become active. Transients to one state or another depend on the
amplitude of the external stimuli. The boundary of the regions of attraction
of stable equilibrium states $O_{1}$ and $O_{11}$ is an invariant manifold $M$
containing equilibrium states $O_{5}$, $O_{8}$, $O_{6}$, $O_{10}$, $O_{7}$,
$O_{9}$ and $O_{27}$. This manifold intersects $\rho_{i}$ planes along the
trajectories belonging to the heteroclinic cycle $G$. Switching sequential
neural activity is most pronounced in cases when an external stimuli throws
the starting point into a small neighborhood of the manifold $M$. It should be
noted that heteroclinic cycles include six (and not three as for the first
system) saddle states of equilibrium.

####  3.2.1 Sequential activity

Three-dimensional system described by Eqs. (7), like Eqs. (2), can exhibit
damping sequential activity. Figures. 14 and 15 show examples of such
activity. Depending on the external stimuli, the following dynamics can be
obtained:

  * •

sequential one-time excitation of all three neurons (Fig. 14 with the
subsequent transition of all neurons into an inactive state defined by the
equilibrium state $O_{1}$ ((a) and (b)) or into the active state determined by
the equilibrium state $O_{11}$ ((c) and (d)).

  * •

sequential multiple excitation of all three neurons (Fig. 15 with subsequent
damping ((a) and (b)) or with the transition to the active state of all
neurons ((c) and (d)).

In addition, as for the first system, activation of one or two neurons is
possible.

Figure 14: Sequential neural activity. The heteroclinic cycle is marked in red
(in Figures (a) and (b)) here and further. Depending on the amplitude of the
external stimuli, after each of the elements is activated sequentially once,
either deactivation occurs (Figure (a) and (b)), or activation of all elements
(Figure (c) and (d)). Parameter values: $\alpha=0.2$, $\beta=1.9$,
$\gamma=0.6$. (a), (b) $a=0.05$, $b=0.4$, $c=0.5857$. (c), (d) $a=0.05$,
$b=0.4$, $c=0.5858$.

Figure 15: Sequential neural activity. Depending on the amplitude of the
external stimuli, after each of the elements is activated sequentially twice,
either deactivation occurs (Figure (a) and (b)), or activation of all elements
(Figure (c) and (d)). Parameter values: $\alpha=0.2$, $\beta=1.9$,
$\gamma=0.6$. (a), (b) $a=0.05$, $b=0.4$, $c=0.585745$. (c), (d) $a=0.05$,
$b=0.4$, $c=0.585746$.

##  4 Summary and conclusions

In this paper

  * •

Two dynamical systems based on the generalized Lotka-Volterra were proposed.
Ensembles consist of excitable elements with excitatory couplings.

  * •

The existence of heteroclinic cycles in the phase space of models is shown.
The presence of a heteroclinic cycle in the phase space makes it possible to
simulate the effect of winnerless competition between modes – sequential
switching of ensemble activity between individual neurons or groups of
neurons.

  * •

In the proposed models, heteroclinic cycles are not stable, which leads to a
complete damping of the ensemble’s activity after some time.

  * •

Depending on the external stimuli, one, two or three neurons can be activated.

  * •

The second model, depending on the amplitude of the external stimuli, can
demonstrate a damping sequential neural activity with an exit to an inactive
state, or to a permanent active state.

The proposed models will be used to construct networks of neural elements and
to study their dynamics.

This paper is devoted to Prof. V. N. Belykh, who is outstanding scientist in
the field of nonlinear dynamics.

## Appendix. Relationship of systems (2) and (7) with the population dynamics
equation

When intraspecific competition is taken into account , the dynamics of
population size can be described by the following equation [22]

$\dot{x}=x\left(\frac{b_{0}x}{N+x}-d_{0}-ex\right),$ (1)

where $x$ is the population size, $b_{0}x$ is the specific rate of
reproduction (fertility), $d_{0}$ is the specific rate of death of individuals
(mortality), $e$ is the coefficient characterizing the influence of
competition on fertility and mortality, $N$ is the population density at which
the average time between consecutive contacts of one the individual is equal
to the minimum time between successive fertilizations of the female.

If we decompose the right side of Eq. (1) into a power series, we get

$\dot{x}=x\left(-d_{0}+\left(\frac{b_{0}}{N}-e\right)x-\frac{b_{0}x^{2}}{N^{2}}+\frac{b_{0}x^{3}}{N^{3}}-\frac{b_{0}x^{4}}{N^{4}}+\frac{b_{0}x^{5}}{N^{5}}-...\right),$
(2)

where the coefficient $\frac{b_{0}}{N}-e>0$.

If in Eq. (2) we restrict ourselves to terms of the second order, then we get
the equation

$\dot{x}=x\left(-d_{0}+\left(\frac{b_{0}}{N}-e\right)x\right)=d_{0}x\left(-1+\frac{b_{0}-eN}{d_{0}N}x\right).$
(3)

By the substitution $\frac{b_{0}-eN}{d_{0}N}x=\rho$, we get

$\frac{d\rho}{dt}=d_{0}\rho(-1+\rho).$ (4)

If we now substitute the variable $t$ with $\frac{\tau}{d_{0}}$, then we get
Eq. (1).

If the equation (2) takes into account terms up to and including the third
order, then we get the equation

$\dot{x}=x\left(-d_{0}+\left(\frac{b_{0}}{N}-e\right)x-\frac{b_{0}x^{2}}{N^{2}}\right).$
(5)

On condition $\frac{b_{0}}{N}-e>\frac{2}{N}\sqrt{b_{0}d_{0}}$ polynomial
$f(x)=-d_{0}+\left(\frac{b_{0}}{N}-e\right)x-\frac{b_{0}x^{2}}{N^{2}}$ has two
different real roots:

$x_{1,2}=\frac{\frac{b_{0}}{N}-e\pm\sqrt{D}}{\frac{2b_{0}}{N^{2}}},$

where $D=\left(\frac{b_{0}}{N}-e\right)^{2}-\frac{4b_{0}d_{0}}{N^{2}}$ is the
discriminant of the polynomial $f(x)$. That is $x_{1}>x_{2}>0$. Then Eq. (5)
can be written as

$\dot{x}=\frac{b_{0}}{N^{2}}x(x_{2}-x)(x-x_{1})=\frac{b_{0}}{N^{2}}x_{1}^{2}x\left(\frac{x_{2}}{x_{1}}-\frac{x}{x_{1}}\right)\left(\frac{x}{x_{1}}-1\right).$
(6)

By the substitution $\frac{x}{x_{1}}=\rho$, we get

$x_{1}\frac{d\rho}{dt}=\frac{b_{0}}{N^{2}}x_{1}^{3}\rho\left(\frac{x_{2}}{x_{1}}-\rho\right)(\rho-1).$
(7)

By the substitution $\frac{x_{1}^{2}b_{0}}{N^{2}}t=\tau$,
$\frac{x_{2}}{x_{1}}=\gamma$ we get Eq. (6).

## Funding

The results in Section 2 were supported by the Ministry of Science and Higher
Education of Russian Federation (Project No. 0729-2020-0036), the results in
Section 3 were supported by RSF grant # 23-12-00180.

## References

  * [1] Rabinovich, M. I., Varona, P., Selverston, A. I., and Abarbanel, H. D., Dynamical principles in neuroscience, Reviews of modern physics, 2006, vol. 78, no. 4, p. 1213.
  * [2] Hahnloser, R. H., Kozhevnikov, A. A., and Fee, M. S., An ultra-sparse code underliesthe generation of neural sequences in a songbird, Nature, 2002, vol. 419, no. 6902, pp. 65–70.
  * [3] Galán, R. F., Sachse, S., Galizia, C. G., and Herz, A. V., Odor-driven attractor dynamics in the antennal lobe allow for simple and rapid olfactory pattern classification, Neural computation, 2004, vol. 16, no. 5, pp. 999–1012.
  * [4] Levi, R., Varona, P., Arshavsky, Y. I., Rabinovich, M. I., and Selverston, A. I., Dual sensory-motor function for a molluskan statocyst network, Journal of neurophysiology, 2004, vol. 91, no. 1, pp. 336–345.
  * [5] Guckenheimer, J., and Holmes, P., Structurally stable heteroclinic cycles, in textitMathematical Proceedings of the Cambridge Philosophical Society, Cambridge University Press, 1988, vol. 103, no. 1, pp. 189–192.
  * [6] Stone, E. and Holmes, P., Random perturbations of heteroclinic attractors, SIAM Journal on Applied Mathematics, 1990, vol. 50, no. 3, pp. 726–743.
  * [7] Postlethwaite, C. M., and Dawes, J. H., Resonance bifurcations from robust homoclinic cycles, Nonlinearity, 2010, vol. 23, no. 3, p, 621.
  * [8] Driesse, R. and Homburg, A. J., Resonance bifurcation from homoclinic cycles, Journal of Differential Equations, 2009, vol. 246, no. 7, pp. 2681–2705.
  * [9] Seliger, P., Tsimring, L. S., and Rabinovich, M. I., Dynamics-based sequential memory: winnerless competition of patterns, Physical Review E, 2003, vol. 67, no. 1, p. 011905.
  * [10] Afraimovich, V. S., Tristan, I., Huerta, R., and Rabinovich, M. I., Winnerless competition principle and prediction of the transient dynamics in a Lotka–Volterra model, Chaos: An Interdisciplinary Journal of Nonlinear Science, 2008, vol. 18, no. 4.
  * [11] Rabinovich, M. I., Zaks, M. A., and Varona, P., Sequential dynamics of complex networks in mind: Consciousness and creativity, Physics Reports, 2020, vol. 883, pp. 1–32.
  * [12] Afraimovich, V. S., Rabinovich, M. I., and Varona, P., Heteroclinic contours in neural ensembles and the winnerless competition principle, International Journal of Bifurcation and Chaos, 2004, vol. 14, no. 4, pp. 1195–1208.
  * [13] Afraimovich, V. S., Zhigulin, V. P., and Rabinovich, M. I., and Varona, P., On the origin of reproducible sequential activity in neural circuits, Chaos: An Interdisciplinary Journal of Nonlinear Science, 2004, vol. 14, no. 4, pp. 1123–1129.
  * [14] Rabinovich, M. I., Varona, P., Tristan, I., and Afraimovich, V. S., Chunking dynamics: heteroclinics in mind, Frontiers in computational neuroscience, 2014, vol. 8, p. 22.
  * [15] Latorre, R., Varona, P., and Rabinovich, M. I., Rhythmic control of oscillatory sequential dynamics in heteroclinic motifs, Neurocomputing, 2019, vol. 331, pp. 108–120.
  * [16] Durstewitz, D. and Deco, G., Computational significance of transient dynamics in cortical networks, European Journal of Neuroscience, 2008, vol. 27, no. 1, pp. 217–227.
  * [17] Muezzinoglu, M. K., Tristan, I., Huerta, R., Afraimovich, V. S., and Rabinovich, M. I., Transients versus attractors in complex networks, International Journal of Bifurcation and Chaos, 2010, vol. 20, no. 06, pp. 1653–1675.
  * [18] Brinkman, B. A., Yan, H., Maffei, A., Park, I. M., Fontanini, A., Wang, J., and La Camera, G., Metastable dynamics of neural circuits and networks, Applied Physics Reviews, 2022, vol. 9, no. 1.
  * [19] Komarov, M. A., Osipov, G. V., and Suykens, J. A. K., Sequentially activated groups in neural networks, Europhysics Letters, 2009, vol. 86, no. 6, p. 60006.
  * [20] Komarov, M. A., Osipov, G. V., and Zhou, C. S., Heteroclinic contours in oscillatory ensembles, Physical Review E, 2013, vol. 87, no. 2, p. 022909.
  * [21] Rabinovich, M. I., Afraimovich, V. S., and Varona, P., Heteroclinic binding, Dynamical Systems, 2010, vol. 25, no. 3, pp. 433–442.
  * [22] Bazykin, A. D., Nonlinear dynamics of interacting populations, World Scientific, 1998.