text:
\begin{center}
{\LARGE Gradient-Like Flows and Self-Indexing}

\vspace{.08in}

{\LARGE in Stratified Morse Theory}

\vspace{.23in}

\large{Mikhail Grinberg\footnotemark}

\vspace{.15in}

\large{August 15, 2000}
\end{center}

\footnotetext{Research partially supported by NSF grant
$\#$ DMS-9971030.}

\vspace{.2in}

\begin{abstract}
We develop the idea of self-indexing and the technology of
gradient-like vector fields in the setting of Morse theory on a
complex algebraic stratification.  Our main result is the local
existence, near a Morse critical point, of gradient-like vector
fields satisfying certain ``stratified dimension bounds up to
fuzz'' for the ascending and descending sets.  As a global
consequence of this, we derive the existence of self-indexing
Morse functions.
\end{abstract}

\vspace{.3in}

\tableofcontents

\vfill

\section{Introduction}

The goal of this paper is to develop the idea of self-indexing and
the technology of gradient-like vector fields in the setting of Morse
theory on a complex algebraic stratification.  Our main result
(Theorem \ref{main}) is the local existence, near a Morse critical
point, of gradient-like vector fields satisfying certain ``stratified
dimension bounds up to fuzz'' for the ascending and descending sets.
As a global consequence of this, we derive the existence of
self-indexing Morse functions (Theorem \ref{siinsm}).

This paper traces its roots to the informal lecture notes
``Intersection Homology and Perverse Sheaves'' by R. MacPherson
\cite{MP}.  These notes outline a vision for developing the theory of
middle perversity perverse sheaves on a stratified complex variety,
starting with a definition of a perverse sheaf in the style of the
Eilenberg-Steenrod axioms.  This approach to perverse sheaves relies
on self-indexing Morse functions as a key technical tool.  We now
proceed to introduce this circle of ideas at a leisurely pace.

\subsection{Classical Morse Theory}

Let $X^d$ be a compact smooth manifold, and let $f: X \to \R$ be a
Morse function.  Write $\Sigma_f$ for the set of critical points of
$X$.  The function $f$ is called self-indexing if $f(p) =
\index_f (p)$ for every $p \in \Sigma_f$.  The significance of this
definition is the following.  A self-indexing $f$ gives rise
canonically to a cochain complex $\cC_f = (C^*, d)$ of $\C$-vector
spaces such that:

(a)  $C^i$ has a natural basis up to sign parameterized by
$\{ \, p \in \Sigma_f \; | \; \index_f (p) = i \, \}$;

(b)  $H^*(\cC_f) = H^*(X; \C)$.

Often, one uses an arbitrary Morse function plus another piece
of structure (e.g., a generic metric) to define a complex  for
computing $H^*(X; \C)$.  However, one needs self-indexing to produce
a canonical complex starting with a Morse function alone.  The
complex $\cC_f$ is called the Morse-Smale complex.  It plays a
central role in Milnor's beautiful exposition \cite{Mi} of Smale's
proof of the $h$-cobordism theorem in dimensions $> 4$.  The first
step of that proof is the following theorem.

\begin{thm}\label{siincm} {\em [16, Theorem B]}
Every compact, smooth  manifold $X$ admits a self-indexing Morse
function $f: X \to \R$.
\end{thm}

We recall an outline of the proof of this in \cite{Mi}.  Start with
any Morse function $g : X \to \R$.  A $\nabla g$-like vector field
$V$ on $X$ is a $C^\infty$ vector field such that:

(a)  $V_p = 0$ for all $p \in \Sigma_g$;

(b)  $V_x \, g > 0$ for all $x \notin \Sigma_g$.

\noindent
(This definition is slightly different from the one in \cite{Mi}.)
Given such a $V$, let $\psi_V : X \times \R \to X$ denote its flow.
Then for $p \in \Sigma_g$, let
$$M^\pm_V (p) = \{ \, x \in X \; | \; \lim_{t \to \mp \infty}
\psi_V (x,t) = p \, \}.$$
The sets $M^\pm_V (p)$ are called the ascending and descending sets
of $p$ (relative to $V$).  Here now are the three main steps of the
proof of Theorem \ref{siincm}.

{\em Step 1: Morse lemma.}  By Morse lemma, there exists a
$\nabla g$-like vector field $V$ such that, for every $p \in \Sigma_g$,
the set $M^-_V (p)$ is a manifold of dimension $\index_g(p)$, and
$M^+_V (p)$ is a manifold of dimension $d - \index_g(p)$.

{\em Step 2: General position.}  By perturbing the vector field $V$,
we can ensure that $M^-_V (p) \cap M^+_V (q) = \emptyset$ whenever
$\index_g (p) \leq \index_g(q)$.

{\em Step 3: Modifying $g$.}  There exists a Morse function
$f : X \to \R$ such that:

(a)  near each $p \in \Sigma_g$, we have $f(x) = g(x) + \index_g(p)
- g(p)$;

(b)  $V$ is $\nabla f$-like.

\noindent
Conditions (a) and (b) guarantee that $f$ is self-indexing.  The role
of the vector field $V$ is to ensure that $\Sigma_f = \Sigma_g$.

The main result of this paper (Theorem \ref{main}) is a replacement
for Step 1 in the complex stratified setting.  Steps 2 and 3 carry
over in a more or less straightforward fashion.

\subsection{Stratified Morse Theory}

Let now $X^d$ be a smooth complex algebraic variety, and let $\cS$ be
an algebraic Whitney stratification of $X$.  Stratified Morse theory,
pioneered by Goresky and MacPherson (see \cite{GM1}, \cite{GM2}), aims
to study the topology of the pair $(X, \cS)$ by means of a real
$C^\infty$ function $f: X \to \R$.  The following definitions go back
to the original paper \cite{GM1}.

\begin{defn}
(i)   Let $p \in X$ be a point contained in a stratum $S$.  We say
that $p$ is critical for $f$ ($p \in \Sigma_f$) if it is critical for
the restriction $f|^{}_S$.

(ii)  For $S \in \cS$, let $\Lambda_S$ be the conormal bundle
$T^*_S X \subset T^*X$, and let $\Lambda = \Lambda_\cS = \bigcup_S
\Lambda_S$.  This $\Lambda$ is called the conormal variety to $\cS$;
it is a closed Lagrangian subvariety of $T^*X$ (the fact that
$\Lambda$ is closed is equivalent to one of the Whitney conditions).
Note that $p \in \Sigma_f$ if and only if $d_p f \in \Lambda$.

(iii) Let $p \in \Sigma_f$.  We say that $p$ is Morse for $f$ if it
is Morse for the restriction $f|^{}_S$ (where $S$ is the stratum
containing $p$) and $d_p f \in \Lambda^0$, the smooth part of $\Lambda$.
The set $\Lambda^0$ is called the set of generic conormal vectors to
$\cS$.
\end{defn}

It is important to note that there is nothing specifically complex
about the above definitions: they would apply equally well to any
Whitney stratification of a real $C^\infty$ manifold.  By contrast,
the following definition is only justified in the complex algebraic
(or analytic) setting.

\begin{defn}
Let $p \in \Sigma_f$ be Morse and let $S$ be the stratum of $p$.
We define
$$\index_f (p) = \index_{f|^{}_S} (p) - \dim^{}_\C S =
\index_{f|^{}_S} (p) + \codim^{}_\C S - d.$$
\end{defn}

The normalization constant $d = \dim_\C X$ is subtracted to make the
possible range of values of the index symmetric about the origin:
$\index_f (p) \in \{ -d, \dots, d \}$.  Ignoring this normalization,
the index is attempting to count the (real) ``descending directions''
of $f$ at $p$, by counting the tangent directions in the obvious way
and then postulating that exactly half of the normal directions are
descending.  The point of this paper is that this method of counting
is justified.  The justification comes from a local result concerning
Morse theory near a point stratum.  Before stating the result itself,
we present a somewhat sharper conjecture.

\begin{conj}\label{mainconj}
Let $X = \C^d$, $p \in X$ be the origin, and $\cS$ be an algebraic
Whitney stratification of $X$ such that $\{p\}$ is a stratum.  Let
$\Lambda^0_p = \Lambda_{ \{ p \} } \cap \Lambda^0$ be the set of
generic covectors at $p$, and let $f \in \Lambda^0_p$.  Regard $f$
as a linear function $f : X \to \R$.  Then there exists a closed
ball $B$ around $p$ and an $\cS$-preserving $\nabla f$-like vector
field $V$ on $B$ such that the ascending and descending sets
$M^\pm_V(p)$ satisfy:
$$\dim^{}_\R M^\pm_V(p) \cap S \leq \dim^{}_\C S \;\;\;
\mbox{for every} \;\;\; S \in \cS.$$
\end{conj}

The reader is referred to Section 2.2 for precise definitions of an
$\cS$-preserving $\nabla f$-like vector field and of the sets
$M^\pm_V(p)$.  However, the truth and the level of difficulty of this
conjecture might be sensitive to the exact class of vector fields
chosen.  Therefore, it is best to view the conjecture as being
somewhat imprecise, with the phrase ``$\cS$-preserving $\nabla f$-like
vector field'' open to interpretation.  Theorem \ref{main} should be
seen as Conjecture \ref{mainconj} ``up to fuzz''.

\begin{thm}\label{main}
Let $X$, $p$, $\cS$ be as in Conjecture \ref{mainconj}, and let
$\Delta \subset \Lambda^0_p$ be an open set.  Then there exit an
$f \in \Delta$ (which we view as a linear function $f : X \to \R$),
a closed ball $B$ around $p$, and a closed, real semi-algebraic
$K \subset B$ such that:

(i)   $\dim^{}_\R K \cap S \leq \dim^{}_\C S$ for every $S \in \cS$;

(ii)  $f^{-1}(0) \cap K = \{ p \}$; and

(iii) for every open $\cU \supset K$, there exists an $\cS$-preserving
$\nabla f$-like vector field $V$ on $B$ with $M^\pm_V(p) \subset \cU$.
\end{thm}

Theorem \ref{main} will be proved in Sections 4 and 5.
It is important to note that for both Conjecture \ref{mainconj}
and Theorem \ref{main} it is crucial that the pair $(X, \cS)$ is
assumed to be complex.  For contrast, consider $X = \R^3$, let $\cS$
be the obvious stratification with six strata whose $2$-skeleton is
the cone $x^2+y^2=z^2$, and take $f(x,y,z) = z$.  Then every
$\cS$-preserving $\nabla f$-like vector field $V$ on $X$ will have
$\dim M^+_V(p) = \dim M^-_V(p) = \dim X = 3$.  (An ice cream cone can
hold a 3d amount of ice cream!)  This example shows that the concept
of index does not generalize well to Morse theory on a real
stratification.  It is also easy to modify the above example to one
with all strata of even dimension (just put the same cone in $\R^4$).
The following consequence of Theorem \ref{main} will be proved in
Section 6.

\begin{thm}\label{siinsm}
Every proper, Whitney stratified complex algebraic variety $(X, \cS)$
admits a self-indexing Morse function $f : X \to \R$.
\end{thm}

\subsection{Perverse Sheaves}

We now explain the motivation behind Theorems \ref{main} and
\ref{siinsm} coming from the theory of middle perversity perverse
sheaves.  Let $(X^d, \cS)$ be as in Section 1.2.  Associated to the
pair $(X, \cS)$ is the category $\cP(X, \cS)$ of middle perversity
perverse sheaves on $X$ constructible with respect to $\cS$.  This
category was first introduced by Beilinson, Bernstein, and Deligne in
the seminal paper \cite{BBD} in 1982.  They defined $\cP(X, \cS)$ as a
full subcategory of $D^b_\cS (X)$, the bounded $\cS$-constructible
derived category of sheaves on $X$.  This way $\cP(X, \cS)$ inherits
an additive structure from $D^b_\cS (X)$.  In fact, it turns out to
be an abelian category, unlike the derived category $D^b_\cS (X)$.
Perverse sheaves are ``simpler'' than the derived category in several
other respects.  For example, they form a stack (i.e., objects and
morphisms are locally defined), which is not the case for
$D^b_\cS (X)$.

For these reasons, it was long felt desirable to have a definition
of $\cP(X, \cS)$ which does not rely on the more complicated
derived category, and which elucidates the very simple formal
properties of perverse sheaves.  Such a definition was proposed by
MacPherson in his AMS lectures in San Francisco in January 1991.
Informal notes of those lectures \cite{MP} produced by MacPherson (and
dated 1990) were distributed at the meeting, but were never published.
A slightly modified version of the same definition has now appeared
in print in the lecture notes \cite{GrM} from the summer institute
held in Park City, Utah in July 1997.

Roughly, MacPherson's definition (the 1997 version) proceeds as
follows.  A perverse sheaf $P \in \cP(X, \cS)$ is an assignment
$(Y,Z) \mapsto H^*(Y, Z; P)$, plus coboundary and pull-back maps,
subject to a set of axioms modeled after the Eilenberg-Steenrod
axioms.  Here $(Y,Z)$ is a {\em standard pair} in $X$, that is a
pair of compact subsets of a certain special kind.  Namely, $Y$ is
a smooth, real $2d$ dimensional submanifold of $X$ with corners of
codimension two, and $Z$ is a union of some of the walls of $Y$,
so that each corner is formed by one wall from $Z$ and one wall
not from $Z$, and all the boundary strata of $Y$ are transverse
to $\cS$.  It is a very pleasant feature of this definition that we
only need to consider standard pairs.  The structure group
$H^*(Y, Z; P)$ is called the cohomology of $Y$ relative to $Z$ with
coefficients in $P$; it is a finite dimensional, graded vector space
over $\C$.

Turning now to the axioms, all but one (the dimension axiom)
represent a more or less straightforward adaptation of the
classical Eilenberg-Steenrod to the setting of pairs inside a
fixed space.  The remaining axiom runs as follows.

\vspace{.1in}

\noindent
{\bf Dimension Axiom:}  Let $f : X \to \R$ be a proper Morse
function, and $c \in \R$ be a value such that $f^{-1}(c)$
contains exactly one critical point: $f^{-1}(c) \cap \Sigma_f =
\{ p \}$.  Then
$$H^i (\{ x \in X \; | \; c-\epsilon \leq f(x) \leq c+\epsilon\}, \,
\{ x \in X \; | \; f(x) = c-\epsilon\}; P) = 0$$
for $i \neq \index_f(p)$ and $0 < \epsilon \ll 1$.

\vspace{1in}

Here is how this definition fits with the idea of self-indexing.
Given a standard pair $(Y,Z)$, one can use a variant of Theorem
\ref{siinsm} to find a self-indexing function $$f: Y \setminus
\{ \mbox{corners} \} \to [-d-1, d+1]$$ which is {\em adapted to}
the pair $(Y,Z)$.  More precisely, this means that
$f^{-1}(-d-1) = Z \setminus \{ \mbox{corners} \}$,
$f^{-1}(d+1) = \partial Y \setminus Z$, and the level sets of $f$
``look like pages of a book'' near the corners.  Now, just as in
classical Morse theory, such a function  $f$ gives rise to a functor
$$\cC_f : \cP(X, \cS) \to \{ \mbox{cochain complexes} \}$$
with the property that $H^*(\cC_f(P)) = H^*(Y, Z; P)$ for every
$P \in \cP(X, \cS)$.  If we write $\cC_f = (C^*, d)$, then the
functor $$C^i : \cP(X, \cS) \to \{ \mbox{vector spaces} \}$$ breaks
up into a direct sum over $\{ \, p \in \Sigma_f \; | \; \index_f (p)
= i \, \}$ of the so-called {\em Morse group} functors
$$M_{d_p f} : \cP(X, \cS) \to \{ \mbox{vector spaces} \}.$$
Each of the $M_{d_p f}$ is localized near $p$, i.e., factors through
the restriction to any open neighborhood of $p$.  (Perverse sheaves
on open subsets of $X$ and restriction functors between them are
defined in the obvious way.)  Moreover, as implied by the notation,
the functor $M_{d_p f}$ depends on $f$ only through the differential
$d_p f \in \Lambda^0$.

The existence of the lift $\cC_f$ from cohomology to cochain
complexes with the above properties, by itself, has important
consequences for the structure of the category $\cP(X, \cS)$.
For example, it immediately implies that morphisms of perverse
sheaves are locally defined.  More is true however: the functor
$\cC_f$ is, in some sense, independent of $f$; the only essential
dependence is on the pair $(Y,Z)$.  To be precise, given any two
self-indexing functions $f_1, f_2$ adapted to $(Y,Z)$, the functors
$\cC_{f_1}$ and $\cC_{f_2}$ are related by a quasi-isomorphism which
is itself canonical up to chain homotopy.  A reader familiar with
Floer homology has, no doubt, anticipated this.  However, unlike in
Floer homology, this independence of the function requires a proof
which is rather more involved than the proof of $d^2 = 0$.
(The later is a straightforward consequence of the axioms.)

Once the independence of the function is established in a suitably
flexible form, all of the (soft) formal properties of perverse
sheaves can be derived from considering the complexes $\cC_f$ for
different functions and standard pairs.  For example, the abelian
property of $\cP(X, \cS)$ follows easily from the abelian property
of the category of cochain complexes.  Finally, one can develop the
theory of perverse sheaves, starting with MacPherson's definition,
to the point of proving the following.

\begin{thm}\label{indep}
MacPherson's definition of $\cP(X, \cS)$ agrees with the original
definition due to Beilinson-Bernstein-Deligne.
\end{thm}

The present author has recently had an opportunity to work through
a proof of Theorem \ref{indep} in a graduate course at MIT, and a
future paper describing this proof is being planned.  The goal of
the present paper is to provide the main geometric tool:
self-indexing Morse functions.  We conclude this introduction with a
conjecture due to MacPherson.

\begin{conj}\label{fary}
Omitting the dimension axiom from MacPherson's definition of
$\cP(X, \cS)$ gives a category which is equivalent to $D^b_\cS (X)$.
\end{conj}

\noindent
{\bf Acknowledgments.}  I would like to thank the IHES in
Bures-sur-Yvette for their hospitality during August 99 when most of
this work was done.  My intellectual debt to Bob MacPherson is clear
from the preceding introduction.  I am also grateful to the students
in my course at MIT for patiently (but not uncritically) allowing me
to think through these ideas out loud.

config: arxiv


text:
\begin{document}

\begin{frontmatter}

\begin{fmbox}
\dochead{Research}

\title{Simulating brain rhythms using an ODE with stochastically
varying coefficients}

\author[
   addressref={aff1},
   email={benjamin.ambrosio@univ-lehavre.fr}
]{\inits{B}\fnm{Benjamin} \snm{Ambrosio}}
\author[
   addressref={aff2},
   email={lsy@cims.nyu.edu}
]{\inits{L-S}\fnm{Lai-Sang} \snm{Young}}

\address[id=aff1]{
  \orgname{Normandie Univ, UNIHAVRE,LMAH,FR CNRS 3335,
ISCN},

  \postcode{76600},
  \city{Le Havre},
  \cny{France}
}
\address[id=aff2]{
  \orgname{Courant Institute of Mathematical Sciences, New York University},

  \postcode{10003},
  \city{New York City},
  \cny{USA}
}

\address[id=aff3]{
  \orgname{Center for Neural Science, New York University},

  \postcode{10003}
  \city{New York City},
  \cny{USA}
}

\begin{artnotes}

\end{artnotes}

\end{fmbox}

\begin{abstractbox}

\begin{abstract}
\parttitle{First part title}
Text for this section.

\parttitle{Second part title}
Text for this section.
\end{abstract}

\begin{keyword}
  neural excitability, signal propagation, neural rhythms
\kwd{sample}
\kwd{article}
\kwd{author}
\end{keyword}

\end{abstractbox}

\end{frontmatter}

\section{Introduction}  Rhythms, or oscillatory patterns of neural activity, occur ubiquitously
 in many parts of the central nervous system. One typically classifies
them by their frequency bands.
For example, $\alpha$-band activity  ($8$-$12$ Hz) is detected in the occipital lobes and is
associated with relaxed wakefulness, $\beta$-band rhythms ($12$-$30$ Hz) are related to muscles and movements,
and  $\gamma$-rhythms ($30$-$90$ Hz) are implicated in  information
transfer and are associated with cognitive processes.
The origins and functional roles of these rhythms are active topics of current research, though much remains to be understood.

This paper is about simulating brain rhythms; the challenge is to generate mathematically signals that resemble the rhythms produced naturally by the brain. We will focus on gamma rhythms, and
for definiteness, we will base our study on gamma-band activity in the visual cortex,
which has been the subject of detailed experimental studies e.g. \cite{Gra-1989,Hen-2005,Xin-2012}
and for which mechanisms are better understood (see below). See also the review article \cite{Car-2016},
and the modeling paper \cite{Cha-2018}.

Experimental data show that there are two aspects to the character of gamma rhythms: one is their oscillatory nature,
and the other is their irregularity. In spite of their being called ``a rhythm", gamma rhythms are
far from being periodic in the sense of mathematics. There is a recognizable beat, to be sure,
but spectral power density studies
show that gamma rhythms are broad-band, with wandering
frequencies and phases. Activity patterns are episodic; the beats are uneven in magnitude,
degrading from time to time before the resumption of oscillatory behavior.
One of the challenges to the dynamical systems community is to understand what kind of underlying dynamics can produce
signals with these characteristics.

Before explaining what is new, let us first review what has been done.
There have been a number of theoretical studies of gamma rhythms. An early model
is PING \cite{Whi-2000,Bor-2003}; similar models include \cite{Erm-1998,
Tie-2001} and others. These models produce highly regular population spikes, capturing
the oscillatory behavior of gamma rhythms well but not their irregular character. Two other much cited works are \cite{Bru-1999,
Bru-2000,Bru-2001}, which modeled gamma dynamics as periodic plus a noise term; that is not quite
right still as it presumes a single intrinsic frequency.
The irregular, episodic nature of gamma rhythms was captured in the multiple-firing events of
\cite{Ran-2013} and further studied in \cite{Cha-2016,Cha-2018} using networks of hundreds to thousands of integrate-and-fire neurons.

Our aim in this paper is not to build a biologically realistic model to reproduce the gamma activity in
real cortex. We seek to model only the phenomenon, or the signal,
not the underlying biology, and to the extent that is possible, to do so analytically using as simple
an equation as we can.
The solution we came up with is defined by
a nonlinear ODE in two variables, with stochastically varying coefficients.
In much the same way that the FitzHugh-Nagumo (FHN)
equations describe oscillatory behavior of a system with a voltage and a recovery variable,
the two variables in our system can be thought of as describing the Excitatory and Inhibitory conductances produced
by a local neuronal population. We found that allowing the coefficients of our ODE to wander
randomly produces results that best match the irregular character of gamma rhythms in real cortex.

The three main sections of this paper are Section 2, which
describes the deterministic model, Section 3, where we introduce the stochasticity, and
Section 4, which examines the spectral properties and  frequencies of the signal produced,
to be compared to corresponding results in gamma-band activity in cortex.

\section{The deterministic model}

In Sect. 2.1 we present the model equations together with basic features of the dynamical
systems they generate. Sects. 2.2 and 2.3 study the dependence on the three parameters that are
allowed to vary.

\subsection{Equations and basic dynamical features}

 To motivate the equations we propose, we review quickly the {\it Recurrent Excitation-Inhibition} (REI) mechanism for gamma-band rhythms proposed in \cite{Cha-2018}. From the crossing of threshold by a few E-cells, recurrent
excitation leads to the elevation of subthreshold activity and the spiking of more E and I-cells,
activating both populations essentially simultaneously. Depending on how the membrane potentials of
neurons postsynaptic to the spiking cells are positioned, I-cells can step in to stop this developing
event quickly, or the recruitment of more cells to join the event continues for a millisecond
or two longer before the suppressive effects of I-spikes are felt. Since I-neurons are fairly densely connected to the rest of the local population, this pushback tends to hyperpolarize a good majority
of the cells. The decay of the suppressive effect, and the depolarization of E-cells,
helped along by the external drive, allows the scenario to be repeated again. This explains
the production of a rhythm. Biophysical time constants cause the cycles to have
frequencies in the gamma band.

In our model $u$ and $v$ represent the absolute values of the
conductances of a typical E and a typical I-neuron in  a local population
We propose that their dynamics are described by the following equations:
\begin{equation}
\label{eq:a-y}
\left\{\begin{array}{rcl}
\epsilon u_t&=&u(-K(u-a_1)(u-a_2)-v)\\
v_t&=&\gamma v(bu-v+c)
\end{array}
\right.
\end{equation}
where  $a_1, a_2, b$ and $c$ are fixed parameters with
\begin{equation}
\label{eq:parameters-ini}
\begin{array}{c}
 a_1=-0.01, \ \ a_2=0.1, \ \ b= 11.9, \ \ c=6.6 \times 10^{-4},
\end{array}
\end{equation}
while $\epsilon, \gamma$ and $K$ are parameters taking values, initially at least, in
$$\epsilon \in [0.01,1], \ \ \gamma \in [1,25] \quad \mbox{and} \quad K\in [30,100]\ .
$$
These equations are inspired by the FHN system, which also has a slow-fast structure.
The cubic nullcline in FHN is replaced by a parabola in the first equation. Factors of $u$ and $v$ are added respectively in front of the first and second equations to constrain the solutions to the upper-right quadrant as well as to better fit conductance oscillations. These modifications are also inspired by the Leslie-Gower model, see \cite{Amb-2018}.

In the rest of this paper, we will adopt the following notation:
\[F(u,v)=u(-K(u-a_1)(u-a_2)-v),\,\ \ G(u,v)=\gamma v(bu-v+c)\]
\[f(u)=-K(u-a_1)(u-a_2),\,\ \ g(u)=(bu+c)\ .\]

The nullclines of Eqns (\ref{eq:a-y}) are then given by
\[u=0, \, v=f(u)\]
for the first equation, and
\[v=0, v=g(u)\]
for the second equation. Note that the  polynomial $f(u)=-K(u-a_1)(u-a_2)$ reaches its maximum for $u=0.5(a_1+a_2)$ with a value of $0.25K(a_2-a_1)^2$.

Our choice of parameters allows only one intersection between the quadratic function $f(u)$
and the linear function $g(u)$ in the positive quadrant. This is ensured by the condition:
\[b\frac{a_1+a_2}{2}+c>0.25(a_2-a_1)^2K\]
which gives
\[K<176.8.\]
There are four stationary points in the positive quadrant, the region of interest. They are
\[(0,0), (0,c), (a_2,0), (u^*,v^*)\]
where $u^*$ is the positive solution of
\[K(u-a_1)(u-a_2)+bu+c=0 \]
and
\[v^*=bu^*+c.\]
See Figure 1, which gives a sense of the global dynamics and basic structures
of the system.

\begin{figure}
\begin{center}
 \label{fig:FP}
\includegraphics[height=9cm, width=12cm]{FP-VF-3.png}

\end{center}
\caption{This figure shows the fixed points, nullclines, and vector field, as well as the
trajectories lying in the sets $u=0$, $v=0$ and the limit-cycle for $K=60$ and $\epsilon=0.1$. For these values of parameter $ (a_2,0)$ is a saddle and $(u^*,v^*)$ is a source.  The two fixed points
$(0,0)$ and $(0,c)$ are too close to be discernible in this figure. In the right square,
we zoom in to visualize them: $(0,0)$  is a source, and $(0,c)$  is a saddle.
}
\end{figure}

The next result confirms that the region depicted in Figure 1 is indeed positively invariant.

\begin{theorem}
\label{th:lc}
The positive $(u,v)$-quadrant is invariant under the dynamics defined by Eq. (\ref{eq:a-y}),
and there exists a bounded absorbing set to which all solutions enter.
\end{theorem}

\begin{proof}
The fact that the positive quadrant is positively invariant follows from the fact that $u=0, v_t=\gamma v(-v+c)$ and $\epsilon u_t=-u(K(u-a_1)(u-a_2)), v=0$, are solutions lying on the $u$ and $v$-axes. For the existence of an absorbing set, we compute
\[\frac{d}{dt}(\epsilon u^2+v^2).\]
We find:
\begin{equation}
\label{eq:absset}
\frac{d}{dt}(\epsilon u^2+v^2)=2\big[-Ku^4+K(a_1+a_2)u^3-Ka_1a_2u^2-u^2v+b\gamma uv^2-\gamma v^3+c\gamma v^2   \big]
\end{equation}
By using Young inequality, for $u,v >0$,
\[uv\leq \frac{u^p}{p}+\frac{v^q}{q}, \, \frac{1}{p}+\frac{1}{q}=1 .\]
We can split the right hand side of \eqref{eq:absset} into a  polynomial in $u$ of degree $4$, and a polynomial in $v$ of degree $3$, with both negative leading coefficients. Since the solution lies in the positive quadrant, by polynomial comparison, we obtain:
\begin{equation}
\label{eq:ieqas}
\frac{d}{dt}(\epsilon u^2+v^2)\leq -A(\epsilon u^2+v^2)+B
\end{equation}
where $A>0$ and $B>0$ can be chosen independently of initial conditions. Integrating \eqref{eq:ieqas} leads to the existence of an absorbing set attracting all trajectories.

\end{proof}

\subsection{Varying the parameter $\epsilon$}

\begin{figure}
 \label{fig:eps=0001-10K=60}
\includegraphics[height=5cm, width=6cm]{eps=0001K=60.png}
\includegraphics[height=5cm, width=6cm]{eps=01K=60.png}\\
\includegraphics[height=5cm, width=6cm]{eps=1K=60.png}
\includegraphics[height=5cm, width=6cm]{eps=10K=60.png}
\label{fig:eps=1K=60}
\caption{This figure illustrates the effect of varying $\epsilon$, from small to large values. The parameter $K$ is set to $K=60$. The top left panel illustrates $\epsilon=0.001$, the top right, $\epsilon=0.1$, bottom left, $\epsilon=1$ and bottom right $\epsilon=10$.}
\end{figure}

Recall that we have three parameters: $\epsilon, \gamma$ and $K$. We first explain the role
of $\epsilon$, fixing for now $\gamma=1$, and studying the dynamics as $\epsilon$ is varied
for each value of $K$. The dynamics of the system from $\epsilon$ very small to very large
for $K=60$ (a fairly typical value of $K$) are summarized in Figure 2.
For $\epsilon \ll 1$, Eq. (\ref{eq:a-y}) describes a slow-fast system
with a limit cycle as can be seen in the top two panels of Fig 2. As $\epsilon$ is increased,
this limit cycle turns into a sink somewhere between $\epsilon=0.1$ and $\epsilon=1$.
For large $\epsilon$, e.g at $\epsilon=10$, one can show that the the critical attractive
manifold is the line $\Delta: v=g(u)$. All the trajectories reach this line fast and then follow it slowly
toward the fixed point $(u^*,v^*)$ as can be seen in the bottom right panel.

Below we will give the analysis for $\epsilon$ very small, as well as the Hopf bifurcation that takes
the limit cycle to the sink.

\subsubsection{The case of $\epsilon \ll 1$}

For $\epsilon$ small enough a slow-fast analysis allows to compute the limit-cycle up to an $O(\epsilon)$ order. The behavior can be described geometrically as follows. We denote by $\mathcal{C}$ the curve  $v=f(u)$. For $\epsilon$ small enough, a trajectory starting from the right side of $\mathcal{C}$ will increase along the curve ($v_t>0$ there) until it reaches the maximum point $A=(\frac{a_1+a_2}{2},f(\frac{a_1+a_2}{2}))$. This is a jump point, see \cite{Kru-2001}, i.e., from there the trajectory leaves $\mathcal{C}$ and goes
at high speed to reach a neighborhood of the point $B=(0,f(\frac{a_1+a_2}{2}))$. After that, since at first $u_t<0$, the trajectory remains stuck near the line $u=0$. It goes down ($v_t<0$) until it crosses the point $(0,f(0))=(0,-Ka_1a_2)$, at which point $u_t$ becomes positive. This is a fold point but not a jump point, see \cite{Kru-2001}.  Dynamics near this point have been analyzed in \cite{Amb-2018}. See also \cite{Wan-2019} and references therein cited. The trajectory continues to follow the axis $u=0$ until it reaches a point $C$ on the axis $u=0$ which is significantly below $(0,f(0))$.  Here there is the possibility of the so-called \textit{canard} phenomenon,  see \cite{Ben-1981,Kru-2001,Szm-2001}.
At $C$, the trajectory leaves the axis $u=0$ and goes very quickly toward the point $D$ on $\mathcal{C}$ with the same ordinate as $C$. This gives a qualitative description of the limit-cycle. For $\epsilon$
sufficiently small, precise statements can be rigorously deduced from Geometrical Singular Perturbation  Theory. Good reviews can be found in \cite{Hek-2010,Jon-1995,Kap-1999,Kru-2001}.

Let $\Gamma'$ be the closed curve defined by:
\[\Gamma'=[A,B]\cup [B,C]\cup [C,D] \cup \zeta \]
where, $\zeta \subset \mathcal{C}$ is the arc from $D$ to $A$.
\begin{theorem}
\label{th:eps-small}
The limit-cycle $\Gamma$ introduced in theorem \ref{th:lc} is $O(\epsilon)$ close of $\Gamma'$.
\end{theorem}
\begin{remark}
We did not considered the uniqueness of the limit-cycle here. For a proof in the case $\epsilon$ small, see \cite{Wan-2019}.
\end{remark}

We point out that the system defined by \eqref{eq:a-y} provides a simple example, in a Neuroscience context, in which canard solutions emerge and can be computed explicitly. Thanks to the polynomial expression of  the vector field, the computations performed in \cite{Amb-2018}  become simpler and explicit around the point $(0,f(0))$. Details are given in the Appendix for the convenience of the reader.

\subsubsection{Hopf bifurcations}

As indicated in Sect. 2.1, the range of $\epsilon$ of interest  is $[0.1, 1]$, and it is in this range
of $\epsilon$ that the limit cycle turns into a sink as shown in Fig 2. We now give more detail
on this bifurcation, specifically the Hopf bifurcation that occurs at $(u^*, v^*)$ where $(u^*, v^*)$
is the unique fixed point in the interior of the positive quadrant; see Sect. 2.1.

\begin{figure}
 \label{fig:FromSinkToLC}
 \includegraphics[height=4cm, width=5cm]{BifDiagKeps.png}\\
\includegraphics[height=4cm, width=4cm]{eps=04K=60.png}
\includegraphics[height=4cm, width=4cm]{eps=036K=60.png}
\includegraphics[height=4cm, width=4cm]{eps=03K=60.png}
\caption{Illustration of the Hopf bifurcation. In the top panel we have plotted the Hopf bifurcation diagram in the $K,\epsilon$ plane. The bottom panels illustrate the bifurcation for $K=60$. As $\epsilon$ is decreased, trajectory ensuing from a fixed initial condition, evolves first towards a spiral sink and then to a limit-cycle. Illustrations, from left to right are for $K=60$, $\epsilon=0.4, 0.36$ and $0.3$.}
\end{figure}

The Jacobian matrix at fixed points is
\[J=\begin{pmatrix}
\frac{1}{\epsilon}(-3Ku^2+2K(a_1+a_2)u-Ka_1a_2-v)& -\frac{1}{\epsilon}u\\
bv & -2v+(bu+c)\
\end{pmatrix}.\]
Substituting in $v^*=bu^*+c$, we obtain at $(u^*,v^*)$, that
\[J_{(u^*,v^*)}=J^*=\begin{pmatrix}
-\frac{1}{\epsilon}Ku^*(2u^*-(a_1+a_2))& -\frac{1}{\epsilon}u^*\\
b(bu^*+c) & -(bu^*+c)
\end{pmatrix},\]
which gives

\[\det(J^*)=\frac{1}{\epsilon}u^*(bu^*+c)(K(2u^*-(a_1+a_2))+b)\]
while
\[tr(J^*) = -\frac{1}{\epsilon}Ku^*(2u^*-(a_1+a_2)) - (bu^*+c)\ .\]
From the above expressions, we deduce the following proposition.

\begin{proposition}
For $K\in [30,100], det(J^*)>0$. It follows that for each $K$ there exists a value of $\epsilon$ at which
a Hopf bifurcation occurs. This value is given by:
\[\epsilon=Ku^*\frac{a_1+a_2-2u^*}{bu^*+c}.\]
\end{proposition}

We close this section by an application of the Poincare-Bendixon theorem to our system.
\begin{theorem}
Each trajectory starting in the region $\{u>0,v>0\}$ either converges to $(u^*,v^*)$ or evolves towards a limit-cycle. For $\epsilon<Ku^*\frac{a_1+a_2-2u^*}{bu^*+c}$, it converges towards a limit-cycle.
\end{theorem}
\begin{proof}
The proof follows from the analysis of the nullclines and the nature of fixed points.

\end{proof}

\subsection{Dependence of dynamics on the parameters $\epsilon, K$ and $\gamma$}

Continuing to keep $\gamma=1$, we first examine the dynamics of
Eq. (\ref{eq:a-y}) as functions of $K$ and $\epsilon$. Simulation results are shown in Fig 4.
Notice first that these results are consistent with those in Fig 3 with regard to increasing
$\epsilon$ for fixed $K$. What is new here is the effect of varying $K$ for each $\epsilon$.
Fig 4 shows clearly that larger $K$ corresponds to larger excursions by $u$ and $v$.
This means
\begin{itemize}
\item[(i)] when solutions are attracted to a limit cycle, the limit cycle has a larger diameter for larger $K$; and
\item[(ii)] whether the $\omega$-limit set is a  limit cycle or a sink, it is located closer to $u=0, v=0$
for smaller values of $K$.
\end{itemize}

\medskip
Finally, we examine the effect of varying $\gamma$. From the equations, it is clear that trajectories
of Eq. (\ref{eq:a-y}) will trace out the same curves as long as $\gamma\epsilon$ remains constant;
and that varying $\gamma$ keeping $\gamma\epsilon$ fixed corresponds to changing the speed with which one moves along these curves. For example, at $K=60$, for values of $\epsilon=0.1$
and $\gamma=1$,  numerical simulation gives a limit cycle with period $\sim 44$  ms (equivalently
frequency around 22 Hz). For $\epsilon=0.01, \gamma=10$, the period becomes $4.4$ ms (frequency around 225 Hz).

\begin{proposition}
\label{prop:freq} For each fixed $K$, the curves traced out by the trajectories of Eq. (\ref{eq:a-y})
depend only on  $\epsilon \gamma$. Fixing $K$ and $\epsilon \gamma$, and varying $\gamma$,
velocities are proportional to $\gamma$; in particular, the frequency of the limit cycle is proportional to $\gamma^{-1}$.
\end{proposition}

\medskip
The meaning and main general effects of the variation of parameters $\epsilon, K$ and $\gamma$ in Eq. (\ref{eq:a-y})
can be summarized as follows:
{\it \begin{itemize}
\item increasing $\epsilon$ changes the dynamical regime from one with a limit cycle
in a slow-fast
system to one with an attractive fixed point;
\item $K$ controls the sizes of the excursion of $(u,v)$ in the system's oscillatory behavior:
in general, the larger $K$, the larger the excursions; while
\item
 for each fixed value of $\epsilon \gamma$, the magnitude of $\gamma$ controls the frequencies
of the limit cycle.
\end{itemize}}

\medskip
As we will show momentarily, these are the parameters we need to vary to produce the irregularities
seen in gamma rhythms.

\begin{figure}
 \label{fig:eps=01-04K=30-100}
\includegraphics[height=4cm, width=5cm]{eps=01K=30-100.png}
\includegraphics[height=4cm, width=5cm]{eps=02K=30-100.png}\\
\includegraphics[height=4cm, width=5cm]{eps=03K=30-100.png}
\includegraphics[height=4cm, width=5cm]{eps=04K=30-100.png}
\caption{This figure gives a panorama of sample trajectories within the paremeters range of interest.
Four panels corresponding to $\epsilon=0.1, 0.2, 0.3$ and $0.4$ are shown. In each panel,
trajectories for different values of $K$ are depicted by different colors:
$K=30$ (red), $50$ (cyan), $70$ (green) and $90$ (purple).}
\end{figure}

\section{The stochastic model}

As discussed in the Introduction, there are two facets to gamma rhythms as observed
in real cortex: one is their oscillatory nature; the other is their irregular, episodic character.
The deterministic system in Section 2 provided the underlying oscillations. Here we create
irregularity by adding randomness to the deterministic model.
We begin by explaining the motivation for what is done.

A standard way to model random dynamical systems is to use a stochastic differential equation (SDE).
These ideas were used in \cite{Bru-1999,Bru-2000}, who studied networks of sparsely coupled
integrate-and-fire neurons, focusing on the limit where system size tends to infinity
and the number of connections is infinitely small compared to system size. Arguing that distinct
neurons are likely to have disjoint sets of presynaptic cells (which is not the case in real cortex),
the authors of
\cite{Bru-1999,Bru-2000} modeled neuronal dynamics by an equation consisting of
a deterministic part describing meanfield activity plus a Gaussian noise that is independent
from neuron to neuron, and gamma rhythms were modeled as regimes following a supercritical
Hopf bifurcation. Another relevant paper is \cite{Bru-2003}. Here the authors assumed that gamma rhythms consisted of random noise
on top of a periodic motion; they were primarily interested in the periodicity.

Since the publication of these theoretical results, new experimental data were obtained, including
the broad-band nature of gamma rhythms in real cortex \cite{Hen-2005}.
With this new understanding, we believe gamma activity might fit better with
dynamics that are not tied to a single period, something assumed in previous works.
As to the nature of the irregularity, it was shown in real data that
gamma power and frequencies {\it wander} \cite{Xin-2012}, i.e. these quantities drift,
but similar gamma patterns are often repeated for tens, sometimes up to one or two hundred,
milliseconds.
The same was observed in simulations
using networks of integrate-and-fire neurons \cite{Cha-2018}, Fig 1, and
the phenomenon can be explained by the
REI mechanism (see Sect. 2.1): In the aftermath of a large firing event that involves a larger than usual
fraction of the E and I-population, inhibitory conductances  tend to
be strong. Disinhibition occurs at roughly the same time for many neurons,
increasing the likelihood that another major event will be precipitated when the march towards
threshold is too synchronized. Similarly, when the rhythm degrades, such as when the nucleation
of a firing event is interrupted too quickly by I-firing, the situation will likely remain ambiguous
for some time before another large event can occur.

To model the wandering character of gamma power and frequency, as well as the occasional
degradation of the rhythm, we propose to use an ODE whose coefficients are not fixed but
are allowed to drift freely and
randomly -- effectively performing random walks -- within certain designated ranges.

Specifically, we allow $K$ to evolve between  $K_{\min}\geq 30$ and $K_{\max}\leq100$. Within this range, each $0.1$ ms, we update $K$ by letting
\[K=K(1+0.1\mathcal{U}(-1,1))\]
where  $\mathcal{U}$ is a random uniform variable in $[-1,1]$. If at some point, the random walk sets the value of $K$ we move backward in the other direction by setting
\[K(1-0.1\mathcal{U}(-1,1)),\]
where the value $\mathcal{U}(-1,1))$ is kept as in the step before.
Analogously, each $0.1$ ms, we update $\epsilon$ by setting
\[\epsilon=\epsilon+0.01\mathcal{U}(-1,1)\]
using an independent uniform realization. If the value $\epsilon=\epsilon+0.01\mathcal{U}(-1,1)$ is outside a range $[\epsilon_{\min},\epsilon_{\max}]$, with $\epsilon_{\min}\geq 0.04$ and  $\epsilon_{\max}\leq 0.1$, we set $\epsilon$ to the value
\[\epsilon=\epsilon-0.01\mathcal{U}(-1,1).\]

Finally, each $0.1$ ms, we evaluate $\gamma$ as
\[\gamma=\gamma+0.1\mathcal{U}(-1,1).\]
and draw repeatedly with the uniform law until $\epsilon \gamma$ belongs in the range $[f_{\min},f_{\max}]$, with $f_{\min}\geq 0.2$ and  $f_{\max}\leq 0.5$.
 To simulate gamma rhythms, the parameters above
are chosen so that most but not all of the time,
the dynamics are oscillatory. Once parameters that
produce suitable qualitative behaviors are located, it is generally simpler to adjust the values
of $u, v$  or the mean frequencies of the oscillations by modifying slightly the two equations of Eq. \eqref{eq:a-y} (e.g. by inserting a scaling coefficient in front).

\begin{figure}
 \label{fig:uvsto}
\includegraphics[height=3.5cm, width=12cm]{utvtsto500-1000.png}\\
\includegraphics[height=3.5cm, width=12cm]{utvtsto1000-1500.png}\\
\includegraphics[height=3.5cm, width=12cm]{utvtsto1500-2000.png}\\
\includegraphics[height=3.5cm, width=12cm]{utvtsto2000-2500.png}\\
\caption{This figure represents the evolution of the stochastic version of the system \ref{eq:a-y}. The time window (in ms) is [500-1000],[1000-1500],[1500-2000], [2000-2500],[2500-3000]. $K \in [30-50],$ $\epsilon \in [0.2, 0.5]$, $\gamma \epsilon \in  [0.2,0.5]$. We plotted
$\bar u =1.96u+0.00672$ in solid purple, to be thought of as representing E-conductance,
$v$ in green, representing I-conductance. Since the ratio of E-current to E-conductance
is roughly three to four times that of the ratio of I-current to I-conductance, we have plotted
also $3.5 \bar u$ in dashed purple. Note the tight relationship between $3.5 \bar u$ and $v$.
}
\end{figure}

 Figure 5 shows a solution of the stochastic version of Eq. \eqref{eq:a-y}
adjusted as indicated to simulate the E and I-conductances
of a typical neuron in a local population under drive thus producing gamma rhythms.
The irregular, episodic nature of the rhythm is
clearly visible. By definition, the E-current entering a neuron is defined to be its E-conductance
times a factor proportional to the distance of membrane potential
to the E-reversal potential, and the same is true for I-currents. As this factor for E-current is
$3$ to $4$ times that for I,
we have also plotted (in dash) a graph that is
$3.5$ times the height of the E-conductance. Modulo a multiplicative constant, then,
the dashed purple and green plots can be thought of as approximations of
 E and I-currents respectively.

We remark on the tightness with which the green plots (I-current) follow the dashed purple plots
(E-current). There is a well known theory of balanced states \cite{Vre-1998}
that asserts that in the limit as system size tends to infinity, E-currents and I-currents are balanced
when averaged over time. But the experimental results of \cite{Oku-2008} and subsequently
the modeling paper \cite{Cha-2018} show that much more than that is true, namely that these
currents are in fact roughly balanced from moment to moment, not just when averaged over time.
This is entirely consistent with the REI mechanism, which tells us that in most
firing events, both E and I-neurons are recruited, with a few E crossing threshold first, hence
E-firing possibly leading by a little bit.
The tight relationship between our dashed purple
and green curves captured well  this phenomenon.

\section{Control of frequencies and spectrum}

\begin{figure}
 \label{fig:FTMoy}
\includegraphics[height=5cm,width=12cm]{TFMoy2.png}\\
\includegraphics[height=5cm,width=12cm]{Spectro.png}
\caption{Panel A. This panel illustrates the mean value of Fourier coefficients computed along the time interval $[0.5s, 2.5s]$ with periods of $200 ms$. In other words, the mean value of Fourier coefficients was computed over time intervals $[t_0+idt,t_0+idt+T]$ with $i\in \{0,20000\}$, $t_0=500ms$, $T=200 ms$ and $dt=0.1$. Original signal comes from simulation of the stochastic system \eqref{eq:a-y}. The main frequency is around 65 Hz, which belongs to the $\gamma$-band frequency.
Panel B.
This panel illustrates the Fourier coefficients of signal on a moving frame. Again, Fourier coefficients are computed over time intervals $[t_0+idt,t_0+idt+T]$ with $i\in \{0,20000\}$, $t_0=500ms$, $T=200 ms$ and $dt=0.1$. The original signal is the same as in panel A. However, here, we do not represent the mean value. Instead, for each $i\in \{0,20000\}$, and for each frequency in hertz in ordinates,  the modulus of coefficients are represented with specific colors.  Again, typical frequencies in the $gamma$-band show up. }
\end{figure}
 For periodic signals, Fourier coefficients provide the right mathematical tool for extracting the main frequencies. To capture the pseudo-periodicity of gamma rhythms, neuroscientists have used
the following computational tool \cite{Hen-2005,Cha-2018}, which we review below.
It is useful but perhaps not standard for mathematicians.
The idea is to fix a time interval of suitable length $T$ and
compute Fourier coefficients on $[0,T]$ as if the signal was periodic with period $T$.
This is repeated with many time shifts, i.e., we sample the signal on other time intervals of
length $T$, replacing $[0,T]$ above by $[t, t+T]$ for
$t=t_0, t_0+dt, t_0+3dt, \cdots$ for small $dt$,
and the computed Fourier coefficients are averaged over all of these samples.
To capture gamma-band frequencies,
the time interval $T$ is usually chosen to be between $200$ and $500$ ms: too short of
an interval will fail to capture the relevant frequencies, and too long of an interval is ineffective
since the signal is not truly periodic.

The usual formula for the Fourier coefficients of a function $f$ is given by:
\begin{equation}
\label{eq:TF-k}
\hat{f}(k)=\frac{1}{T}\int_0^Tf(t)e^{-2i\pi k \frac{t}{T}}dt,
\end{equation}
 and the power concentrated at frequency $k$ is defined to be $|\hat f(k)|^2$.
With Neuroscience applications in mind, one might want to deal with frequencies in hertz.
In this case, a frequency $k$ on a $[0,T]$ period gives a frequency of $n=\frac{k}{T}$ in hertz. One can then adopt the following notation:
\begin{equation}
\label{eq:TF-n-1}
\hat{f}(nT)=\frac{1}{T}\int_0^Tf(t)e^{-2i\pi nT \frac{t}{T}}dt.
\end{equation}
which gives:
\begin{equation}
\label{eq:TF-n-2}
\hat{f}(nT)=\frac{1}{T}\int_0^Tf(t)e^{-2i\pi n t}dt.
\end{equation}
Since we deal with discrete signals with time step $\Delta t$, the computation turns into
\begin{equation}
\label{eq:TF-k-d}
\hat{f}(k)=\frac{1}{N}\sum_{j=0}^{N-1}f(j\Delta t)e^{-2i\pi k \frac{j\Delta t}{T}} , \quad
k \in \{0,1,...,N-1\}, \, N=\frac{T}{\Delta t}
\end{equation}
or
\begin{equation}
\label{eq:TF-k-d}
\hat{f}(nT)=\frac{1}{N}\sum_{j=0}^{N-1}f(j\Delta t)e^{-2i\pi n j\Delta t}, \quad
n \in \{0,1/T,...,(N-1)/T\}\ .
\end{equation}

 A plot of the mean value of the modulus of the Fourier coefficients according to formula \eqref{eq:TF-k-d}, over 2 seconds, for the signals generated by the system Eq. \eqref{eq:a-y}
is given in Fig 5A. The broad-band nature of the signal is quite evident. In Panel B, we show
the value of those coefficient of the signal as a function of time, instead of averaging over the samples
taken at different time points as is shown in Panel A. This plot shows the wandering nature of
the frequencies and power of the signal from Eq. \eqref{eq:a-y}. It captures well the phenomenon
depicted in \cite{Xin-2012} and \cite{Cha-2018}.

\section{Summary and Discussion}

Our larger aims are (i) to help make biology quantitative and (ii) to bring biological
ideas to mathematics, to dynamical systems  in particular.

With regard to (i), our challenge was to simulate gamma-band activity. We went beyond
that to simulate the E and I-conductances of typical neurons, which according to the
REI mechanism is responsible for producing gamma rhythms in a driven
local population. We  proposed a two dimensional model inspired by the FHN and Leslie-Gower
models, providing both theoretical and numerical analyses that give a panorama of the
expected dynamical behavior in a range of parameters. We then introduced some
stochasticity in the model that allows variation in amplitude, phase and frequencies,
in a way that is consistent with real electric signals of the brain. Because solutions of
our stochastic model are realistic, they can be used in future theoretical studies to simulate
input currents to neurons under drive. More generally, the techniques demonstrated
are generalizable to simulations of other rhythmic activity in biological systems.

With regard to (ii), we seek to contribute to the repertoire of dynamical systems theory by
bringing to the community's attention phenomena that are important that occur
frequently in biology.
Dynamical systems theory is well positioned to study periodic or quasi-periodic behavior (as in
KAM theory); there is also a fairly well developed theory of chaotic systems. Rhythms of the brain
are produced entirely naturally, and they are neither periodic nor chaotic, but somewhere in between.
They occur as a result of a push-and-pull mechanism that is typical of the interaction between
 opposing groups of agents in biology. A dynamical theory of biological rhythms, or more
 generally of the interplay between opposing groups of agents, will enrich the subject.

\section{Appendix}
In this appendix, we give details about the analysis near the point $C$. As discussed in section $2$, for $\epsilon$ very small, and the range of parameters under consideration,  the slow-fast analysis allows to locate the limit cycle: trajectories follow the attractive parts of the critical manifolds, exit the attractive cubic part of the critical manifold  through the jump point, and follow the fast fibers outside of the critical manifold. This gives the main picture. However, the point $C$ at which the trajectory exits the part of the critical manifold $u=0$ is significantly below the point $(0,-Ka_1a_2)$ at which the critical manifold becomes unstable. The trajectory follows the unstable part of the critical manifold along a significant distance, which is referred to the so called canard phenomenon. The aim of this appendix is to detail quite simple  computations estimating this distance and providing the coordinates of $C$. The final result is given by the theorem \ref{th:eps} at the end of the appendix. To facilitate the reading, we prepare the theorem with a few propositions.   In the first proposition, we rewrite system \eqref{eq:a-y} around the fold point $(0,-Ka_1a_2)$.

\begin{proposition}
Around the fold point $C=(0,-Ka_1a_2)$, system \eqref{eq:a-y} rewrites:
\begin{equation}
\label{eq:a-chv}
\left\{\begin{array}{rcl}
\dot{x}&=&x(Kx(a_1+a_2)+x^2-y)\\
\dot{y}&=&\epsilon\big(-Ka_1a_2(c+Ka_1a_2)-bKa_1a_2x+(c+2Ka_1a_2)y+bxy-y^2\big)\\
\end{array}
\right.
\end{equation}
\end{proposition}
\begin{proof}
The result follows from the change of variables
\[u=x,\, v=-Ka_1a_2+y.\]

\end{proof}
Next, we apply the change of variables:
\begin{equation}
\label{eq:chv2}
x=rx_2,\, y=r^2y_2,\, \epsilon=r^3.
\end{equation}
The following proposition holds:
\begin{proposition}
After change of variables  equation \eqref{eq:a-chv} writes
\begin{equation}
\label{eq:a-chv2}
\left\{\begin{array}{rcl}
\dot{x}_2&=&K(a_1+a_2)x_2^2+rx^3_2-ry_2\\
\dot{y}_2&=&-Ka_1a_2(c+Ka_1a_2)-bKa_1a_2rx_2+(c+2Ka_1a_2)r^2y_2+br^3x_2y_2-r^4y_2^4\\

\end{array}
\right.
\end{equation}
\end{proposition}
\begin{proof}
The result follows from the change variables \eqref{eq:chv2} and change of time $\tau=rt$.
\end{proof}

Next, consider equation \cref{eq:a-chv2} with $r=0$, $i.e.$
\begin{equation}
\label{eq:K2r=0}
\left\{\begin{array}{rcl}
\dot{x}_2&=&K(a_1+a_2)x_2^2\\
\dot{y}_2&=&-Ka_1a_2(c+Ka_1a_2)\\
\end{array}
\right.
\end{equation}

 The following result, which follows from explicit computation, holds:

\begin{proposition}
 The solution of system \eqref{eq:K2r=0} is:
\begin{equation}
\label{eq:sol:K2r=0}
\begin{array}{rcl}
x_2(t)&=&\frac{1}{x_2^{-1}(0)-K(a_1+a_2)t}\\
y_2(t)&=&y_2(0)-Ka_1a_2(c+Ka_1a_2)t\\
\end{array}
\end{equation}
i.e.
\begin{equation*}
x_2(t)=\frac{1}{\frac{1}{x_2(0)}+\frac{a_1+a_2}{a_1a_2(c+Ka_1a_2)}(y_2(t)-y_2(0))}
\end{equation*}
or
\begin{equation*}
y_2(t)=y_2(0)-\frac{a_1a_2(c+Ka_1a_2)}{a_1+a_2}(\frac{1}{x_2(0)}-\frac{1}{x_2(t)}
\end{equation*}
It follows that orbits solutions of \cref{eq:sol:K2r=0} have the following properties:
\begin{enumerate}
\item Every orbit has a horizontal asymptote given by
\[ y_2 = y_2(0)-\frac{a_1a_2(c+Ka_1a_2)}{a_1+a_2}\frac{1}{x_2(0)}\]
\item  Every orbit has a  vertical asymptote $x_2= 0^{+}$.

\end{enumerate}

\end{proposition}
Assume that a trajectory crosses the line $y=0$ $(v=-Ka_1a_2)$ at
\[x(0)=k\epsilon +o(\epsilon)\]
Let
\[\bar{x}=\frac{a_2}{2},\]
and let $\bar{y}$ the ordinate at which the trajectory crosses the line $x=\bar{x}$.
The following theorem holds.

\begin{theorem}
\label{th:eps}
Assume that the trajectory crosses the line $y=0$ at
\[x(0)=k\epsilon +o(\epsilon)\]
then
\[\lim_{\epsilon \rightarrow 0}\bar{y}=-\frac{a_1a_2(c+Ka_1a_2)}{k(a_1+a_2)}\]
\end{theorem}
\begin{proof}
We have
\begin{equation*}
\begin{array}{rcl}
\bar{y}&=&\epsilon^{\frac{2}{3}}\bar{y}_2\\[0.2cm]
&=&\epsilon^{\frac{2}{3}}\big(\frac{a_1a_2(c+Ka_1a_2)}{a_1+a_2}(\frac{1}{\bar{x}_2}-\frac{1}{x_2(0)})+\epsilon^{\frac{1}{3}}O(1) \big)\\[0.2cm]
&=&\epsilon^{\frac{2}{3}}\big(\frac{a_1a_2(c+Ka_1a_2)}{a_1+a_2}(\frac{\epsilon^{\frac{1}{3}}}{\bar{x}}-\frac{\epsilon^{\frac{1}{3}}}{k\epsilon +o(\epsilon)})+\epsilon^{\frac{1}{3}}O(1) \big)\\[0.2cm]
&=&\frac{a_1a_2(c+Ka_1a_2)}{a_1+a_2}(\frac{\epsilon}{\bar{x}}-\frac{1}{k +h(\epsilon)})+\epsilon O(1)  \mbox{ where } \lim_{\epsilon\rightarrow 0} h(\epsilon)=0\\
\end{array}
\end{equation*}
which gives the result.
\end{proof}

\section*{Acknowledgments} BA would like to thank R\'egion Normandie France and ERDF (European Regional Development Fund) XTERM, and CNRS IEA for funding.

\bibliographystyle{bmc-mathphys}
\bibliography{references}


TITLE: What are the elements in the subgroup generated by $\langle(1234)\rangle$ in $S_4$?
QUESTION [0 upvotes]: What are the elements in the subgroup generated by $\langle(1234)\rangle$ in $S_4$?

REPLY [2 votes]: Consider the element $(1234)$. It has order $4$, hence $<(1234)>$ has 4 elements : the identity, $(1234)$, $(1234)^2=(13)(24)$, $(1234)^3=(1432)$.


TITLE: Convolution product of arithmetic functions
QUESTION [0 upvotes]: An arithmetic function is a real-valued function whose domain is the set of positive integers. Define the convolution product of two arithmetic functions $f$ and $g$ to be the arithmetic function $f \star g$, where $$(f \star g)(n) = \sum_{ij = n}f(i)g(j), \text{ and } f^{\star k} = f \star f \star \cdots \star f \text{ } (k \text{ times}).$$ We say that two arithmetic functions $f$ and $g$ are dependent if there exists a nontrivial polynomial of two variables $P(x,y) = \sum_{i,j} a_{ij}x^iy^j$ with real coefficients such that $$P(f,g) = \sum_{i,j}a_{ij}f^{\star i} \star g^{\star j} = 0$$ and say that they are independent if they are not dependent. Let $p$ and $q$ be two distinct primes and set $$f_1(n) = \begin{cases}1 \quad \text{if } n = p,\\0 \quad \text{otherwise}; \end{cases} \qquad \qquad f_2(n) = \begin{cases}1 \quad \text{if } n = q,\\0 \quad \text{otherwise}.\end{cases}$$ Prove that $f_1$ and $f_2$ are independent.

My book says that $$f_{1}^{\star i}\star f_{2}^{\star j}(n) = \begin{cases}1 \,\, \text{if} \,\, n=p^{i}q^{j} \\ 0 \,\, \text{otherwise} \end{cases}$$ and the result follows from that, but how do they get that and how does the result follow from that?

REPLY [1 votes]: Let us write the convolution $$f_1\star f_1 (n) = \sum_{i\cdot j =n}f_1(i)f_1(j).$$
The function $f_1(i)\ne 0\iff i=p$, hence the only way to get a nonzero product inside the sum above is to put $i=j=p$. This implies that
$$f_1^{\star 2}(n)=\begin{cases}1,&n=p^2,\\0,&n\ne p^2.\end{cases}$$
By using the same technique, you will get that
$$f_1^{\star k}(n)=\begin{cases}1,&n=p^k,\\0,&n\ne p^k.\end{cases}$$
Can you now apply the same technique to $f_1^{\star k}\star f_2^{\star p}$?


TITLE: Showing that $(y-c)(1-F(y)) \underset{y \rightarrow \infty}{\rightarrow} 0$ for $F$ the cdf of an integrable random variable?
QUESTION [0 upvotes]: I was reading the comments in this question and it was suggested that this limit is $0$ by dominated convergence theorem. I'm not able to see how. If anyone sees it and could give some hints?
Thank you.

REPLY [1 votes]: It suffices to show that $yP(Y>y)\to 0$ as $y\to\infty$ as $cP(Y>y)\to 0$ is clear. To see this note that for $y>0$
$$
0\leq yI(Y>y)\leq YI(Y>y).
$$
Take expectations to yield that
$$
yP(Y>y)\leq E(YI(Y>y))\to 0
$$
as $y\to \infty$ by the dominated convergence theorem.


TITLE: A cone with diameter $12$cm and height $8$cm. Find the volume of the inscribed Sphere.
QUESTION [1 upvotes]: A cone with diameter $12\ cm$ and height $8\ cm$. Find the volume of the inscribed Sphere. Can someone help me solve this maths problem?

REPLY [1 votes]: The cone has radius 6, height 8 and slant height 10. I leave it to you to draw a diagram involving Pythagorean triplet proportions of sides 3:4:5 ;
Now the centre of circle touching two lines (slanted & base) is to be connected to base corner. It must be a bisector... whose property is ... proportion of segments on either side of the bisector to be same:
$$ \frac {6} {10} = \frac {x}  {8-x} $$
which fixes radius of base as 3 units. Now you can take it further to find volume.


TITLE: $(a+b\sqrt{2})^n =(c+d\sqrt{3})^m =>ab=cd=0$.
QUESTION [1 upvotes]: If $a,b,c,d \in Q, m,n \in N^*$ and $(a+b\sqrt{2})^n =(c+d\sqrt{3})^m$ then to show that
$ab=cd=0$.
An idea to solve it by contradiction but ...

REPLY [1 votes]: Assume that $AB\neq0$.
Hint: $\mathbb{Q}[\sqrt{2}] \cap \mathbb{Q} [\sqrt{3}] = \mathbb{Q}$.

 This tells us that $(a+b\sqrt{2})^n = p \in \mathbb{Q}$.

Hint: Consider conjugates.

 Taking the conjugates, we get that $(a-b\sqrt{2})^n = p$. Taking the difference, we get that $ 2b \sum {n \choose 2i} a^{n-2i}b^{i}2^i = 0$.

If $n$ is even, we have a non-zero ($2b$) multiplied by a (positive) sum of squares, which can't be 0. This is a contradiction.

 If $n$ is odd, we have $2ab \sum{n \choose 2i} a^{n-2i-1} b^i 2^i = 0 $. Once again, we have a non-zero ($2ab$) multiplied by a (positive) sum of squares, which can't be 0. This is a contradiction.

Hence, $AB=0$. Similarly, $CD=0$.


TITLE: Show that $\mathbb{Q}(\sqrt{1+\sqrt{3}})\cap \mathbb{Q}(\sqrt{1-\sqrt{3}})=\mathbb{Q}(\sqrt{3})$.
QUESTION [2 upvotes]: Show that $\mathbb{Q}(\sqrt{1+\sqrt{3}})\cap \mathbb{Q}(\sqrt{1-\sqrt{3}})=\mathbb{Q}(\sqrt{3})$.

I know that by closure, $\sqrt{3}\in \mathbb{Q}(\sqrt{1+\sqrt{3}})\cap \mathbb{Q}(\sqrt{1-\sqrt{3}})$, but what about the other containment? I want to use the order of the extensions. I know the degrees of each one are $4$ over the field $\mathbb{Q}$, and I know that $\mathbb{Q}(\sqrt{1+\sqrt{3}})$ and $ \mathbb{Q}(\sqrt{1-\sqrt{3}})$ are not equal to each other. But where can I go from here? How do I show that $[\mathbb{Q}(\sqrt{1+\sqrt{3}})\cap \mathbb{Q}(\sqrt{1-\sqrt{3}}):\mathbb{Q}(\sqrt{3}]=1$?

REPLY [1 votes]: Let $\alpha:=\sqrt{1+\sqrt{3}}$ and $\beta:=\sqrt{1-\sqrt{3}}$. You know that
$$[\Bbb{Q}(\alpha):\Bbb{Q}]=[\Bbb{Q}(\beta):\Bbb{Q}]=4,$$
and that $\Bbb{Q}(\sqrt{3})\subset\Bbb{Q}(\alpha)\cap\Bbb{Q}(\beta)$. Because the degree is multiplicative over towers of field extensions
$$[\Bbb{Q}(\alpha):\Bbb{Q}]=[\Bbb{Q}(\alpha):\Bbb{Q}(\alpha)\cap\Bbb{Q}(\beta)]\cdot[\Bbb{Q}(\alpha)\cap\Bbb{Q}(\beta):\Bbb{Q}(\sqrt{3})]\cdot[\Bbb{Q}(\sqrt{3}):\Bbb{Q}],$$
and because you know the two fields $\Bbb{Q}(\alpha)$ and $\Bbb{Q}(\beta)$ are not the same, you have
$$[\Bbb{Q}(\alpha):\Bbb{Q}(\alpha)\cap\Bbb{Q}(\beta)]>1,$$
and clearly $[\Bbb{Q}(\sqrt{3}):\Bbb{Q}]=2$. Because $[\Bbb{Q}(\alpha):\Bbb{Q}]=4$ it follows that
$$[\Bbb{Q}(\alpha):\Bbb{Q}(\alpha)\cap\Bbb{Q}(\beta)]=2
\qquad\text{ and }\qquad
[\Bbb{Q}(\alpha)\cap\Bbb{Q}(\beta):\Bbb{Q}(\sqrt{3})]=1,$$
which means that $\Bbb{Q}(\alpha)\cap\Bbb{Q}(\beta)=\Bbb{Q}(\sqrt{3})$.