# On the $L^{r}$-differentiability of Two Lusin Classes and a Full Descriptive
Characterization of the $HK_{r}$-integral

PAUL MUSIAL  Chicago State University
9501 South King Drive, Chicago, Illinois, 60628 USA<EMAIL_ADDRESS>,
VALENTIN SKVORTSOV  Lomonosov Moscow State University, Mathematics Department
and Moscow Center of Fundamental and Applied Mathematics
Moscow, 119991 Russia<EMAIL_ADDRESS>,  PIOTR SWOROWSKI  Casimirus the
Great University
Powstańców Wielkopolskich 2, 85-090 Bydgoszcz, Poland<EMAIL_ADDRESS>and  FRANCESCO TULONE  University of Palermo
via Archirafi 34, 90123 Palermo, Italy<EMAIL_ADDRESS>
###### Abstract.

It is proved that any function of a Lusin-type class, the class of
$ACG_{r}$-functions, is differentiable almost everywhere in the sense of a
derivative defined in the space $L^{r}$, $1\leq r<\infty$. This leads to
obtaining a full descriptive characterization of a Henstock-Kurzweil-type
integral, the $HK_{r}$-integral, which serves to recover functions from their
$L^{r}$-derivatives. The class $ACG_{r}$ is compared with the classical Lusin
class $ACG$ and it is shown that a continuous $ACG$-function can fail to be
$L^{r}$-differentiable almost everywhere.

###### 2020 Mathematics Subject Classification:

Primary 26A39, 26A46

##  1\. Introduction

The well-known classical definition of the Lebesgue integral in terms of
absolute continuity of the indefinite integral is an example of a so-called
full descriptive characterization of an integral. Definitions of this type are
known also for a number of non-absolute generalizations of the Lebesgue
integral. The most well-known among them is Lusin’s definition of Denjoy
integral (see [15]).

In general, we say that a descriptive characterization of a constructive
integration process is obtained if a class of functions is found so that it
coincides with the class of the indefinite integrals of all functions
integrable in the sense of this process, and if functions of the class are
differentiable almost everywhere in a certain sense, the derivative being
equal to the integrand almost everywhere. If the differentiability (in the
appropriate sense, corresponding to differentiability property of the
indefinite integral) of functions of the class is an additional assumption put
on the class, then we have only a partial descriptive characterization of the
considered integral. If on the other hand the corresponding differentiability
property is a consequence of the description of the class and not an
additional requirement, such a characterization is said to be a full
descriptive characterization. For example, each continuous function from
Lusin’s class $ACG$ is differentiable almost everywhere only in the sense of
the approximate derivative, and this class gives a full descriptive
characterization of the wide Denjoy integral (the $D$-integral, see [15]). At
the same time if we restrict this class by considering only continuous and
differentiable (in the usual sense) $ACG$-functions then we obtain a partial
descriptive characterization of the so-called Khintchine integral (see [5]).
Descriptive characterizations of other non-absolute integrals are discussed in
[2, 13, 14, 17, 18, 19].

In this paper we obtain a full descriptive characterization of a Henstock-
Kurzweil-type integral, the $HK_{r}$-integral (defined in 2004 by Musial and
Sagher [6]), which serves to integrate a derivative defined in the space
$L^{r}$, $1\leq r<\infty$. This $L^{r}$-derivative was introduced in [1] by
Calderón and Zygmund to be used in some estimates for solutions of elliptic
partial differential equations. The $HK_{r}$-integral was defined as an
extension of a Perron-type integral, the $P_{r}$-integral, which was defined
earlier by L. Gordon [3] and which also recovers a function from its
$L^{r}$-derivative. The $HK_{r}$-integral turned out to be strictly wider than
the $P_{r}$-integral (see [7]). It was also shown in [6] that the indefinite
$HK_{r}$-integral is $L^{r}$-differentiable almost everywhere and belongs to a
Lusin-type class of $ACG_{r}$-functions. Some other properties of these
integrals are investigated in [9, 10, 11, 12].

Another class, also giving a descriptive characterization of the
$HK_{r}$-integral, was obtained in [8] in terms of absolute continuity of so-
called $L^{r}$-variational measure generated by a function belonging to
$L^{r}$. Both classes, the class of $ACG_{r}$-functions and the class of
functions generating absolutely continuous $L^{r}$-variational measure, were
known to coincide with the class of the indefinite $HK_{r}$-integrals, but
only under the additional assumption that the functions are
$L^{r}$-differentiable almost everywhere (see [8]). The problem of
$L^{r}$-differentiability almost everywhere of functions of those classes was
left open in [8]. So in fact only partial descriptive characterizations of the
$HK_{r}$-integral were obtained in [6] and [8].

The main aim of the present paper is to prove that any $ACG_{r}$-function is
$L^{r}$-differentiable almost everywhere and thereby obtain a full descriptive
characterization of the $HK_{r}$-integral (see Theorem 13 in Section 3). In
Section 4, we construct an example of a function from a Lusin-type class,
$[ACG]$, which, even under an additional assumption of continuity, fails to
have the property of being $L^{r}$-differentiable almost everywhere for any
$r$. Thereby we show that the class $[ACG]$ is not contained in the class
$ACG_{r}$ for any $r$. As a byproduct of this example we get a new proof,
based on the $L^{r}$-derivative, of a result, previously obtained in [16],
stating that the $HK_{r}$-integral does not cover the wide Denjoy integral.

##  2\. Preliminaries

We recall here the main definitions and facts related to the notion of the
$HK_{r}$-integral (see [3, 6]).

We work in a fixed closed interval $[a,b]$. An interval $I$ is a non-
degenerate closed subinterval of $[a,b]$. We assume in all these definitions
that the functions $F$ belong to the class $L^{r}=L^{r}[a,b]$, $1\leq
r<\infty$. We start with $L^{r}$-derivates and the $L^{r}$-derivative.

######  Definition 1.

The upper-right $L^{r}$-derivate of $F$ at $x$, denoted by $D_{r}^{+}F(x)$, is
defined as the infimum of all numbers $\alpha$ such that

(1) $\left(\frac{1}{h}\int_{0}^{h}[F(x+t)-F(x)-\alpha t]_{+}^{r}\mkern
1.5mudt\right)^{\\!1/r}\\!=\,o(h)\quad\text{as}\ h\rightarrow 0^{+}.$

If no real number $\alpha$ satisfies (1), we set $D_{r}^{+}F(x)=+\infty$. The
lower-right, upper-left, and lower-left $L^{r}$-derivates of $F$ at $x$ are
defined correspondingly as

$\displaystyle D_{+,r}F(x)=-D_{r}^{+}(-F)(x),\quad
D_{r}^{-}F(x)=D_{r}^{+}(-F)(-x),$ $\displaystyle\text{and}\quad
D_{-,r}F(x)=-D_{r}^{+}F(-x).$

If all four $L^{r}$-derivates of $F$ at $x$ are equal, the common value,
denoted by $F^{\prime}_{r}(x)$, is the $L^{r}$-derivative of $F$ at $x$. If
$F^{\prime}_{r}(x)$ is finite we say that $F$ is $L^{r}$-differentiable at $x$
and the value of $F^{\prime}_{r}(x)$, say $\alpha$, is uniquely defined by

(2) $\left(\frac{1}{h}\int_{-h}^{h}|F(x+t)-F(x)-\alpha t|^{r}\mkern
1.5mudt\right)^{\\!1/r}\\!=\,o(h)\quad\text{as}\ h\to 0^{+}.$

We need the following property of $L^{r}$-derivates (see [3, Corollary of
Theorem 4])

######  Proposition 2.

If the two upper (or the two lower) $L^{r}$-derivates of a function $F$ are
less than $+\infty$ (resp. greater than $-\infty$) on a set $E$, then almost
everywhere on $E$ the $L^{r}$-derivative of $F$ exists and is finite.

For completeness we give the following definitions; see [4] for a full
treatment.

######  Definition 3.

Let $F\colon[a,b]\to\mathbb{R}$ and let $c\in[a,b]$. The function $F$ is said
to be approximately continuous at $c$ if there exists a measurable set $E$
containing $c$ such that $c$ is a point of density of $E$ and $F|_{E}$ is
continuous at $c$.

######  Definition 4.

Let $F\colon[a,b]\to\mathbb{R}$ and let $c\in[a,b]$. The function $F$ is said
to be approximately differentiable at $c$ if there exists a measurable set $E$
containing $c$ such that $c$ is a point of density of $E$ and

$\lim_{x\rightarrow c,x\in E}\frac{F(x)-F(c)}{x-c}$

is finite. The value of this limit is the approximate derivative of $F$ at
$c$, and we will denote it by $F^{\prime}_{\textup{ap}}(c)$. Note that the
value of $F^{\prime}_{\textup{ap}}(c)$ does not depend on the choice of $E$.

To define the $HK_{r}$-integral [6] we need some auxiliary notions. A tagged
interval is a pair $(I,x)$ where a point $x\in I$ is a tag. We say that tagged
intervals $(I^{\prime},x^{\prime})$ and $(I^{\prime\prime},x^{\prime\prime})$
are nonoverlapping if intervals $I^{\prime}$ and $I^{\prime\prime}$ are
nonoverlapping, i.e., they have no inner points in common. A partition is any
finite collection $\pi$ of pairwise nonoverlapping tagged intervals. A gauge
is a strictly positive function $\delta$ on $[a,b]$ (or on a subset of
$[a,b]$). We say that a tagged interval $(I,x)$ is $\delta$-fine if
$I\subset(x-\delta(x),x+\delta(x))$.

A partition is $\delta$-fine if all its elements are $\delta$-fine. A
partition $\pi$ is tagged in a set $E\subset[a,b]$ if $x\in E$ for each
element $(I,x)$ of $\pi$.

The Lebesgue measure on $[a,b]$ will be denoted by $\mu$. We recall the
definition of the $L^{r}$-Henstock–Kurzweil integral given in [6].

######  Definition 5.

A function $f\colon[a,b]\rightarrow\mathbb{R}$ is $L^{r}$-Henstock–Kurzweil
integrable (HKr-integrable) on $[a,b]$ if there exists a function $F\in L^{r}$
such that for any $\varepsilon>0$ there exists a gauge $\delta$ such that for
any $\delta$-fine partition $\\{([c_{i},d_{i}],x_{i})\\}_{i}$ we have

$\sum\limits_{i}\left(\dfrac{1}{d_{i}-c_{i}}\int_{c_{i}}^{d_{i}}\lvert
F(y)-F(x_{i})-f(x_{i})(y-x_{i})\rvert^{r}\,dy\right)^{1/r}<\varepsilon.$

######  Definition 6.

[6] Let $E\subset[a,b]$ and $F\in L^{r}$. We say that $F\in AC_{r}(E)$ if for
each $\varepsilon>0$ there exist $\eta>0$ and a gauge $\delta$ defined on $E$
such that for any $\delta$-fine partition $\\{([c_{i},d_{i}],x_{i})\\}_{i}$
tagged in $E$ and such that

(3) $\sum_{i}(d_{i}-c_{i})<\eta$

we have

(4)
$\sum_{i}\left(\frac{1}{d_{i}-c_{i}}\int_{c_{i}}^{d_{i}}|F(y)-F(x_{i})|^{r}\,dy\right)^{\\!1/r}<\,\varepsilon.$

We say that $F$ is an $ACG_{r}$-function if
$[a,b]=\bigcup_{n=1}^{\infty}E_{n}$ where $F\in AC_{r}(E_{n})$ for all $n$.

######  Definition 7.

Let $E\subset[a,b]$. We say that $F\in AC(E)$ if for each $\varepsilon>0$
there exist $\eta>0$ such that for any finite collection of nonoverlapping
intervals $\\{[c_{i},d_{i}]\\}_{i}$ having endpoints in $E$ and such that

$\sum_{i}(d_{i}-c_{i})<\eta$

we have

$\sum_{i}|F(d_{i})-F(c_{i})|<\,\varepsilon.$

We say that $F$ is an ACG-function if $[a,b]=\bigcup_{n=1}^{\infty}E_{n}$,
where $F\in AC(E_{n})$ for all $n$. If, moreover, all $E_{n}$ can be chosen
closed, we say $F$ is an [ACG]-function.

######  Definition 8.

A function $f\colon[a,b]\to\mathbb{R}$ is said to be D-integrable (integrable
in the wide Denjoy sense) if there exists a continuous $ACG$-function
$F\colon[a,b]\to\mathbb{R}$ such that $F^{\prime}_{\textup{ap}}(x)=f(x)$ at
almost all $x\in[a,b]$. We then define $\int_{a}^{b}f=F(b)-F(a)$.

######  Definition 9.

A function $f\colon[a,b]\to\mathbb{R}$ is said to be Kubota integrable if
there exists an approximately continuous $[ACG]$-function
$F\colon[a,b]\to\mathbb{R}$ such that $F^{\prime}_{\textup{ap}}(x)=f(x)$ at
almost all $x\in[a,b]$. We then define $\int_{a}^{b}f=F(b)-F(a)$.

We will use the following property of $ACG_{r}$-functions.

######  Proposition 10.

[6, Corollary 1] If $F$ is an $ACG_{r}$-function, then $[a,b]$ can be
represented as $[a,b]=\bigcup_{n=1}^{\infty}E_{n}$ with $F\in AC(E_{n})$ for
each $n$.

Recall that $F\in L^{r}$ and so is measurable. Hence the above proposition
implies that any $ACG_{r}$-function is a measurable $VBG$-function, and so by
the Denjoy-Khintchine Theorem (see [15, Chapter VII, Section 4, Theorem 4.3])
we obtain:

######  Theorem 11.

Any $ACG_{r}$-function is approximately differentiable a.e. on $[a,b]$.

##  3\. $L_{r}$-differentiability of $ACG_{r}$-functions

We now state and prove the main theorem of this paper.

######  Theorem 12.

Any $ACG_{r}$-function $F$ is $L_{r}$-differentiable a.e. on $[a,b]$.

###### Proof.

According to Theorem 11, $F^{\prime}_{\textup{ap}}$ exists and is finite
almost everywhere on $[a,b]$. We will show that

(5) $F_{r}^{\prime}(x)=F_{\textup{ap}}^{\prime}(x)\ \ \text{a.e.}$

By Lusin’s theorem on $C$-property it is enough to prove (5) on a closed set
$E$ of positive measure such that $F$ restricted to $E$ is continuous. We can
suppose that $F$ is approximately differentiable at each point of $E$. Note
that the approximate derivative and all $L_{r}$-derivates are measurable (see
[15, 3]). To simplify computation we consider first the case where
$F^{\prime}_{\textup{ap}}$ is non-negative on $E$. We start by proving that in
this case

(6) $D_{r}^{+}F(x)=F_{\textup{ap}}^{\prime}(x)\ \ \text{a.e.\ on $E$}.$

If not, then having in mind the relation between $L^{r}$-derivates and
approximate derivates (see [3, Theorem 2]) we get

(7) $D_{r}^{+}F(x)>F_{\textup{ap}}^{\prime}(x)\geq 0$

on a set $E^{\prime}\subset E$ of positive measure. We can suppose that
$E^{\prime}$ is closed. By Definition 4 we can find a set
$E^{\prime\prime}\subset E^{\prime}$ of positive outer measure where $F\in
AC_{r}(E^{\prime\prime})$. By Definition 6 for $\varepsilon=1$ there exist on
$E^{\prime\prime}$ a gauge $\delta$ and a positive number $\eta$ such that
inequality (4) holds for each $\delta$-fine partition
$\\{([c_{i},d_{i}],x_{i})\\}_{i}$ tagged in $E^{\prime\prime}$ and satisfying
(3). Now having fixed $\delta$ corresponding to the chosen $\varepsilon$ we
can find a set $T\subset E^{\prime\prime}$ of positive outer measure on which
$\delta(x)\geq\gamma$ for some constant $\gamma>0$. Let $\overline{T}$ be the
closure of $T$ and let $W$ be the set of two-sided limit points of
$\overline{T}$. We then have $\overline{T}=W\cup C$ where $C$ is countable. As
the set $E^{\prime}$ is closed we have $W\subset\overline{T}\subset
E^{\prime}\subset E$ and $\mu(W)>0$.

We show now that the inequality (4) with $\varepsilon=2$, i.e., the inequality

(8)
$\sum_{i}\left(\frac{1}{d_{i}-c_{i}}\int_{c_{i}}^{d_{i}}|F(y)-F(x_{i})|^{r}\,dy\right)^{\\!1/r}<\,2$

holds for any $\gamma$-fine partition $\\{([c_{i},d_{i}],x_{i})\\}_{i}$ tagged
in $W$, and satisfying (3) with the same $\eta$ which was found for
$\varepsilon=1$ and for $\delta$-fine partitions tagged in $E^{\prime\prime}$.
Indeed, consider such a $\gamma$-fine partition. By the definition of $W$ and
by the continuity of $F$ with respect to $E$, we can choose, for each $i$, a
point $x_{i}^{\prime}\in T$ so close to $x_{i}$ inside of $[c_{i},d_{i}]$ that
$([c_{i},d_{i}],x_{i}^{\prime})$ is also $\gamma$-fine and
$\sum_{i}|F(x_{i})-F(x^{\prime}_{i})|<1$. Then

$\displaystyle\sum_{i}\left(\frac{1}{d_{i}-c_{i}}\int_{c_{i}}^{d_{i}}|F(y)-F(x_{i})|^{r}\mkern
1.5mudy\right)^{\\!1/r}$
$\displaystyle\qquad\leq\sum_{i}\left(\frac{1}{d_{i}-c_{i}}\int_{c_{i}}^{d_{i}}|F(y)-F(x_{i}^{\prime})|^{r}\mkern
1.5mudy\right)^{\\!1/r}$
$\displaystyle\qquad\qquad+\sum_{i}\left(\frac{1}{d_{i}-c_{i}}\int_{c_{i}}^{d_{i}}|F(x_{i}^{\prime})-F(x_{i})|^{r}\mkern
1.5mudy\right)^{\\!1/r}\\!<1+1=2.$

So (8) must be true for any $\gamma$-fine partition tagged in $W$ and
satisfying (3). We now construct such a $\gamma$-fine partition for which (8)
fails to be true. This will give us a contradiction proving (6).

Let $P\subset W$ be the set of density points of $W$. We have $\mu(P)=\mu(W)$.
Let $x$ be any point of $P$. There exists a set $S_{x}$ having $x$ as a
density point, such that

(9) $\lim_{\begin{subarray}{c}t\to 0\\\ x+t\in
S_{x}\end{subarray}}\\!\frac{F(x+t)-F(x)}{t}=F_{\textup{ap}}^{\prime}(x)\geq
0.$

Note that $x$ is also a density point for the set $W\cap S_{x}$. Moreover, we
can find a closed set $A_{x}\subset W\cap S_{x}$ for which $x$ is still a
density point. Let the open set $B_{x}$ be the complement of $A_{x}$ in
$[a,b]$. Choose the number $l$ so that
$D^{+}_{r}(x)>l>F_{\textup{ap}}^{\prime}(x)$. Then by Definition 1 there
exists a constant $c_{x}>0$ and a sequence $h_{k}\to 0^{+}$ such that for each
$k$

(10) $\left(\frac{1}{h_{k}}\int_{0}^{h_{k}}[F(x+t)-F(x)-lt]_{+}^{r}\mkern
1.5mudt\right)^{\\!1/r}\\!>\,c_{x}h_{k}.$

We now put some restrictions on the size of $h_{k}$. We cover $P$ by an open
set $G\subset(a,b)$ such that $\mu(G)<2\mu(P)$, and we assume that all
intervals $[x,x+h_{k}]$ are subsets of $G$. As
$A_{x}\cap[x,x+h]\subset\\{x+t:F(x+t)-F(x)<lt\\}$ for small enough $h$, we can
also assume that for each $h_{k}$

(11) $\int_{0}^{h_{k}}[F(x+t)-F(x)-lt]_{+}^{r}\mkern 1.5mudt\\\
=\int_{B_{x}\cap[x,x+h_{k}]}[F(u)-F(x)-l(u-x)]_{+}^{r}\mkern 1.5mudu.$

Note that this equality in particular implies that in our further
consideration we can assume that for all $k$,

(12) $\mu(B_{x}\cap[x,x+h_{k}])>0.$

Indeed, otherwise

$\int_{0}^{h_{k}}[F(x+t)-F(x)-lt]_{+}^{r}\,dt=0$

contradicting the claim that $D^{+}_{r}F(x)>l$.

Fix some natural $n>\max{\\{2\mu(P)/\eta,2/\mu(P)\\}}$. As $x$ is a point of
dispersion for the set $B_{x}$ we can find $h_{x}\leq\gamma$ such that

(13)
$\frac{\mu(B_{x}\cap[x,x+h])}{h}<\min{\biggl{\\{}\frac{1}{n},\frac{c_{x}^{r}}{n^{r}}\biggr{\\}}}$

holds for $0<h\leq h_{x}$. So we can assume that $h_{k}<h_{x}\leq\gamma$ for
all $k$. Having in mind (10), (11), (12), and (13), we obtain:

$\displaystyle\left(\frac{1}{\mu(B_{x}\cap[x,x+h_{k}])}\int_{B_{x}\cap[x,x+h_{k}]}[F(u)-F(x)-l(u-x)]_{+}^{r}\mkern
1.5mudu\right)^{\\!1/r}$
$\displaystyle\quad=\left(\frac{h_{k}}{\mu(B_{x}\cap[x,x+h_{k}])}\right)^{\\!1/r}\\!\cdot\,\left(\frac{1}{h_{k}}\int_{0}^{h_{k}}[F(x+t)-F(x)-lt]_{+}^{r}\mkern
1.5mudt\right)^{\\!1/r}$ (14)
$\displaystyle\qquad\qquad\geq\biggl{(}\frac{n^{r}}{c^{r}_{x}}\biggr{)}^{\\!1/r}\,\cdot
c_{x}h_{k}=nh_{k}.$

As $B_{x}$ is an open set we have the representation

$B_{x}\cap(x,x+h_{k})=\bigcup_{j}\,(a_{j},b_{j}).$

Note that $a_{j}\in A_{x}\subset W$ and $b_{j}-a_{j}<h_{x}\leq\gamma$ for each
$j$. So $([a_{j},b_{j}],a_{j})$ are $\gamma$-fine intervals tagged in $W$.
Using the following obvious inequality for positive numbers

$\left(\frac{\sum_{k}a_{k}}{\sum_{k}b_{k}}\right)^{\\!1/r}\leq\left(\sum_{k}\frac{a_{k}}{b_{k}}\right)^{\\!1/r}\leq\sum_{k}\biggl{(}\frac{a_{k}}{b_{k}}\biggr{)}^{\\!1/r}$

(this is true for infinite sums provided the series are convergent), we get

(15)
$\displaystyle\begin{split}&\left(\frac{1}{\mu(B_{x}\cap[x,x+h_{k}])}\int_{B_{x}\cap[x,x+h_{k}]}[F(u)-F(x)-l(u-x)]_{+}^{r}\mkern
1.5mudu\right)^{\\!1/r}\\\
&\qquad\leq\left(\sum_{j}\frac{1}{b_{j}-a_{j}}\int_{a_{j}}^{b_{j}}[F(u)-F(x)-l(u-x)]_{+}^{r}\mkern
1.5mudu\right)^{\\!1/r}\\\
&\qquad\leq\left(\sum_{j}\frac{1}{b_{j}-a_{j}}\int_{a_{j}}^{b_{j}}|F(u)-F(a_{j})|^{r}\mkern
1.4mudu\right)^{\\!1/r}\\\
&\qquad\leq\sum_{j}\left(\frac{1}{b_{j}-a_{j}}\int_{a_{j}}^{b_{j}}|F(u)-F(a_{j})|^{r}\mkern
1.5mudu\right)^{\\!1/r}.\end{split}$

(We have taken into account that $F(a_{j})-F(x)<l(a_{j}-x)<l(u-x)$ for
$u\in(a_{j},b_{j})$). Combining (15) with (14) we get

(16)
$\sum_{j}\left(\frac{1}{b_{j}-a_{j}}\int_{a_{j}}^{b_{j}}|F(u)-F(a_{j})|^{r}\mkern
1.5mudu\right)^{\\!1/r}\geq nh_{k}.$

The family of intervals $\\{[x,x+h_{k}]\\}$ for all $x\in P$ and all $k$ forms
a Vitali cover of $P$. Applying the Vitali covering theorem we choose a
countable sequence of nonoverlapping intervals $[x_{j},x_{j}+h_{k_{j}}]$ which
are subsets of $G$ and cover $P$ up to a set of measure zero. Hence
$\mu(P)\leq\,\sum_{j}h_{k_{j}}\leq\mu(G)<2\mu(P)$. We can apply (16) to each
of these intervals and to the corresponding set
$B_{x_{j}}\cap[x_{j},x_{j}+h_{k_{j}}]$. Collecting all intervals constituting
$B_{x_{j}}\cap[x_{j},x_{j}+h_{k_{j}}]$ for all $j$ we obtain a family
$\\{(I_{i},z_{i})\\}_{i}$ of nonoverlapping $\gamma$-fine intervals tagged in
$P$, and so in $W,$ the total length of which is estimated, due to (13), by

$\sum_{j}\mu(B_{x_{j}}\cap[x_{j},x_{j}+h_{k_{j}}])<\sum_{j}\frac{h_{k_{j}}}{n}\leq\frac{2\mu(P)}{n}<\eta.$

So for any finite subfamily of $\\{(I_{i},z_{i})\\}_{i}$ the inequality (8)
should hold. At the same time by (16) and by the choice of $n$ we have

$\sum_{i}\left(\frac{1}{|I_{i}|}\int_{I_{i}}|F(u)-F(z_{i})|^{r}\,du\right)^{\\!1/r}>n\mu(P)>2,$

which must hold also for a certain finite subfamily of
$\\{(I_{i},z_{i})\\}_{i}$ (i.e., a $\gamma$-fine partition tagged in $W$).
This contradiction proves that our assumption (7) is false and (6) is true for
the case of non-negative $F^{\prime}_{\textup{ap}}$.

In the same way using the definition of the upper-left $L^{r}$-derivate, we
obtain a contradiction to the assumption that

(17) $D_{r}^{-}F(x)>F_{\textup{ap}}^{\prime}(x)\geq 0$

holds on a set of positive measure. In this case, instead of equality (11),
the equality

(18) $\int_{0}^{h_{k}}[F(x-t)-F(x)+lt]_{-}^{r}\mkern 1.5mudt\\\
=\int_{B_{x}\cup[x-h_{k},x]}[F(u)-F(x)-l(u-x)]_{-}^{r}\mkern 1.5mudu$

is used for respective $h_{k}$. This time the intervals
$\\{(a_{i},b_{i})\\}_{i}$ constituting $B_{x}\cap[x-h_{k},x]$ give rise to a
$\gamma$-fine partition $\\{(I_{i},z_{i})\\}_{i}$ tagged in $W$. Eventually we
obtain that:

(19) $D_{r}^{-}F(x)=F_{\textup{ap}}^{\prime}(x)\quad\text{a.e.\ on $E$}.$

(again in the case of non-negative $F^{\prime}_{\textup{ap}}$). Having (6) and
(19) we apply Proposition 2 and prove the statement of the theorem for the
case of non-negative $F^{\prime}_{\textup{ap}}$.

If $F^{\prime}_{\textup{ap}}$ is negative on a set of positive measure on
which (5) is false we apply the argument above to the function $-F$. ∎

As we have already noted, the above theorem shows that the class of
$ACG_{r}$-functions coincides with the class of indefinite integrals of all
$HK_{r}$-integrable functions, and we obtain a full descriptive
characterization of the $HK_{r}$-integral:

######  Theorem 13.

A function $f$ is HKr-integrable on $[a,b]$ if and only if there exists an
$ACG_{r}$-function $F$ such that $F_{r}^{\prime}=f$ a.e.; the function
$F(x)-F(a)$ being the indefinite HKr-integral of $f$.

##  4\. The $ACG$ property and $L_{r}$-differentiability

The following statement is proved in [16, Corollary 13]:

######  Proposition 14.

Each indefinite HKr-integral is an [ACG]-function.

Using this proposition and equality (5) for the indefinite $HK_{r}$-integral
together with its approximate continuity (see [6]), we get:

######  Theorem 15.

Every HKr-integrable function on $[a,b]$ is Kubota integrable on $[a,b]$ and
the values of these integrals coincide.

Once again using Proposition 14 and Theorem 13 we obtain a slight improvement
of Proposition 10.

######  Corollary 16.

Each $ACG_{r}$-function is an [ACG]-function.

We now show that an $ACG$-function, even a continuous one, may fail to be
$L_{r}$-differentiable a.e. This, in particular, will imply that the inclusion
stated in Corollary 16 is proper.

######  Theorem 17.

There exists a continuous $ACG$-function which is not $L_{r}$-differentiable
on a set of positive measure for any $r$.

###### Proof.

We let $r_{0,1}=[0,1]$. We remove from $r_{0,1}$ the open interval $u_{0,1}$
centered at $1/2$ having length $u_{0}=1/6$, leaving two closed intervals
$r_{1,1}$ and $r_{1,2}$ each having length $r_{1}=5/12$. We continue this
process in such a way that at level $n$ we are left with $2^{n}$ closed
intervals $\\{r_{n,k}\\}_{1\leq k\leq 2^{n}}$ each having length
$r_{n}:=(n+4)/(2^{n+1}(n+2))$. In order to accomplish this, from each interval
$r_{n,k}$ the concentric open interval $u_{n,k}$ of length
$u_{n}:=1/(2^{n}(n+2)(n+3))$ is removed. We have

$\left|\bigcup_{k=1}^{2^{n}}r_{n,k}\right|=\frac{1}{2}+\frac{1}{n+2}$

Let

$P=\bigcap_{n=1}^{\infty}\bigcup_{k=1}^{n}r_{n,k}$

$P$ is a symmetric, perfect, Cantor-like subset of $[0,1]$ having measure
$1/2$. Let $v_{n,k}$ be the open interval concentric to $u_{n,k}$ having
length $v_{n}:=u_{n}/2$. Define a function $F$ taking the value $0$ on $P$ and
$(-1)^{n}/n$ on $v_{n,k}$. Extend $F$ linearly on the subintervals of
$u_{n,k}$ adjacent to $v_{n,k}$.

$F$ is obviously continuous. It is $AC(P)$ and $AC(u_{n,k})$ for each $n$ and
$k$, and so it is $ACG$ on the whole interval $[0,1]$.

We will now show that $F$ fails to be $L_{r}$-differentiable at any point of
$P$ which is not right-hand isolated in $P$. Let $x$ be such a point. There
are arbitrarily large $n\in\mathbb{N}$ such that for some $h>0$ and $k$, $x+h$
is the midpoint of $v_{n,k}$. Let $l>0$ be such that $x+l$ is the left
endpoint of such $v_{n,k}$ and let $\alpha\in\mathbb{R}$. One can estimate

$\frac{1}{h}\left(\frac{1}{h}\int_{0}^{h}|F(x+t)-F(x)-\alpha
t|^{r}\,dt\right)^{\\!1/r}\\\
>\frac{1}{h}\left(\frac{1}{h}\int_{l}^{h}|F(x+t)-\alpha
t|^{r}\,dt\right)^{\\!1/r}\geq\frac{1}{h}\left(\frac{1}{h}\int_{l}^{h}|F(x+t)|^{r}\,dt\right)^{\\!1/r}$

when $\alpha\geq 0$ and $n$ is odd or $\alpha\leq 0$ and $n$ is even. Continue
the estimate from below:

$\displaystyle\geq\frac{2}{r_{n}}\left(\frac{2}{r_{n}}\cdot\frac{v_{n}}{2n^{r}}\right)^{\\!1/r}=\frac{2^{1-1/r}}{nr_{n}}\left(\frac{u_{n}}{r_{n}}\right)^{\\!1/r}$
$\displaystyle=2^{n+2}\cdot\frac{n+2}{n(n+4)^{1+1/r}(n+3)^{1/r}},$

which is unbounded in $n$. This proves that the function $F$ fails to be
$L_{r}$-differentiable on a set of measure $1/2$. ∎

As we have already mentioned the constructed example gives another proof of
the following result obtained in [16, Theorem 14]:

######  Corollary 18.

There exists a function which is D-integrable but which is HKr-integrable for
no $r$.

Indeed, for $F$ from the proof of Theorem 17 suppose
$f=F^{\prime}_{\textup{ap}}$ is $HK_{r}$-integrable. Then the indefinite
$HK_{r}$-integral $G=\int f$ is an $[ACG]$-function (Theorem 15) with
$G^{\prime}_{r}(x)=f(x)$ a.e. in $[0,1]$ (Theorem 13). Since
$G^{\prime}_{\textup{ap}}(x)=f(x)$ a.e. in $[0,1]$, due to monotonicity
property of the class of approximately continuous $[ACG]$-functions, $F-G$ is
constant, so that $F$, like $G$, is a.e. $L^{r}$-differentiable, giving a
contradiction.

## References

  * [1] A. P. Calderón, A. Zygmund, Local properties of solutions of elliptic partial differential equations. Studia Math. 20, 171–225 (1961)
  * [2] V. Ene, Real Functions – Current Topics, Lecture Notes in Mathematics 1603, Springer-Verlag, 1995.
  * [3] L. Gordon, Perron’s integral for derivatives in $L^{r}$. Studia Math. 28, 295–316 (1966/1967)
  * [4] R. A. Gordon, The Integrals of Lebesgue, Denjoy, Perron, and Henstock, Grad. Stud. Math. 4, Amer. Math. Soc., 1994.
  * [5] A. J. Khintchine, Sur le procede d’integration de M. Denjoy. Matematiczeskij Sbornik 30, 548–557 (1918)
  * [6] P. Musial, Y. Sagher, The $L^{r}$ Henstock–Kurzweil integral. Studia Math. 160(1), 53–81 (2004)
  * [7] P. Musial, V. Skvortsov, F. Tulone, The HKr-integral is not contained in the Pr-integral. Proc. Amer. Math. Soc. 150(5), 2107–2114 (2022)
  * [8] P. Musial, V. Skvortsov, F. Tulone, On descriptive characterization of an integral recovering a function from its $L^{r}$-derivative. Mathematical Notes 111(3), 414–422 (2022)
  * [9] P. Musial, V. Skvortsov, F. Tulone, Comparison of the $P_{r}$-integral with Burkill’s integrals and some applications to trigonometric series. Journal of Mathematical Analysis and Applications 523(1), article number 127019 (2023)
  * [10] P. Musial, F. Tulone, Integration by parts for the $L^{r}$ Henstock–Kurzweil integral. Electronic Journal of Differential Equations 2015(44), 1–7 (2015)
  * [11] P. Musial, F. Tulone, Dual of the class of $HK_{r}$-integrable functions. Minimax Theory and its Applications 4(1), 151–160 (2019)
  * [12] P. Musial, F. Tulone, The $L^{r}$-variational integral. Mediterranean Journal of Mathematics 19(3), article number 96 (2022)
  * [13] W. F. Pfeffer, The Lebesgue and Denjoy–Perron integrals from a descriptive point of view. Ricerche Mat. 48(2), 211–223 (1999)
  * [14] D. Preiss, B. Thomson, The approximate symmetric integral. Canadian Journal of Mathematics 41(3), 508–555 (1989)
  * [15] S. Saks, Theory of the Intergal, Dover Publications, 2005.
  * [16] V. Skvortsov, P. Sworowski, On the relation between Denjoy–Khintchine and HKr-integrals, arXiv:2211.09269, November 2022.
  * [17] V. A. Skvortsov, F. Tulone, Multidimensional dyadic Kurzweil–Henstock- and Perron-type integrals in the theory of Haar and Walsh series, Journal of Mathematical Analysis and Applications, 421(2), 1502–1518 (2015)
  * [18] V. A. Skvortsov, F. Tulone, Multidimensional P-adic integrals in some problems of harmonic analysis, Minimax Theory and its Applications, 2(1), 153–174 (2017)
  * [19] B. Thomson, On VBG Functions and the Denjoy–Khintchine Integral, Real Analysis Exchange 41, 1–54 (2015)


eurm10 msam10 definition[theorem]Definition Networks of piecewise linear
neural mass models–References

# Networks of piecewise linear neural mass models

S.COOMBES  Y.M.LAI  M. ŞAYLI  R.THUL.  Centre for Mathematical Medicine and
Biology, School of Mathematical Sciences, University of Nottingham,
Nottingham, NG7 2RD, UK

(2016)

###### Abstract

Neural mass models are ubiquitous in large scale brain modelling. At the node
level they are written in terms of a set of ordinary differential equations
with a nonlinearity that is typically a sigmoidal shape. Using structural data
from brain atlases they may be connected into a network to investigate the
emergence of functional dynamic states, such as synchrony. With the simple
restriction of the classic sigmoidal nonlinearity to a piecewise linear
caricature we show that the famous Wilson-Cowan neural mass model can be
explicitly analysed at both the node and network level. The construction of
periodic orbits at the node level is achieved by patching together matrix
exponential solutions, and stability is determined using Floquet theory. For
networks with interactions described by circulant matrices, we show that the
stability of the synchronous state can be determined in terms of a low-
dimensional Floquet problem parameterised by the eigenvalues of the
interaction matrix. Moreover, this network Floquet problem is readily solved
using linear algebra, to predict the onset of spatio-temporal network patterns
arising from a synchronous instability. We further consider the case of a
discontinuous choice for the node nonlinearity, namely the replacement of the
sigmoid by a Heaviside nonlinearity. This gives rise to a continuous-time
switching network. At the node level this allows for the existence of unstable
sliding periodic orbits, which we explicitly construct. The stability of a
periodic orbit is now treated with a modification of Floquet theory to treat
the evolution of small perturbations through switching manifolds via the use
of saltation matrices. At the network level the stability analysis of the
synchronous state is considerably more challenging. Here we report on the use
of ideas originally developed for the study of Glass networks to treat the
stability of periodic network states in neural mass models with discontinuous
interactions.

###### keywords:

General applied mathematics; Synchronisation; Nonsmooth equations; Complex
networks; Neural networks.

††volume: 000

##  1 Introduction

The Wilson-Cowan model [1, 2] is one of the most well-known neural mass models
for modelling the activity of cortex, and for a historical perspective see
[3]. Neural mass models generate brain rhythms using the notion of population
firing rates, aiming to side-step the need for large scale simulations of more
realistic networks of spiking neurons. Although they are not derived from
detailed conductance based models they can be motivated by a number of
phenomenological arguments [4], and typically take the form of systems of
nonlinear ordinary differential equations (ODEs). The Wilson-Cowan neural mass
model describes the dynamics of two interacting populations of neurons, one of
which is excitatory and the other inhibitory. Interactions are mediated
between the populations with the use of a nonlinear sigmoidal firing rate
function. In its most simple incarnation it consists of two nonlinear ODEs,
and as such has been widely studied using techniques from phase-plane analysis
and numerical bifurcation theory. At the network level the model can either be
posed on a graph or a continuous space, and since the 1970s there has been a
large amount of attention devoted to the analysis of these models and their
application in neuroscience [5]. Recent examples of their use include
reconciling information from anatomical and functional data [6], understanding
phase-amplitude coupling (whereby the amplitude of a higher frequency brain
rhythm is modulated by the phase of lower frequency activity) [7], modelling
epilepsy [8], and understanding the emergence of cortical resonant frequencies
[9]. Indeed there are many variants of the Wilson-Cowan neural mass model now
in use for interpreting neuroimaging data [10], including those of Zetterberg
et al. [11], Jansen and Rit [12], and Liley et al. [13]. Neural mass models
are a key component of the Virtual Brain project that aims to deliver the
first simulation of the human brain based on individual large-scale
connectivity [14]. Such large-scale brain network models are especially
relevant to understanding resting state networks [15], whereby different
regions of the brain’s sensorimotor system oscillate slowly and synchronously
in the absence of any explicit task.

However, at heart it is well to note that from a mathematical modelling
perspective all neural mass models to date are essentially low dimensional
coupled ODEs with a sigmoidal firing rate nonlinearity, exemplified by the
Wilson-Cowan model. Using extensions of the techniques originally developed by
Amari [16], the continuum or neural field [4] Wilson-Cowan model has been
analysed when the choice of this firing rate nonlinearity is a Heaviside
function. This has been possible because of a smoothing of the firing rate
with a spatial kernel representing anatomical connectivity. However, when
posed on a graph, representing a network of interacting neural populations, no
such smoothing arises. Surprisingly there are hardly any mathematical results
for such networks, as opposed to their continuum counterparts for which there
are now a plethora ranging from the properties of localised states through to
travelling waves, as reviewed in [17]. This discrepancy is really a reflection
of the fact that there are many more techniques for studying smooth dynamical
systems as opposed to nonsmooth. However, the body of mathematical work in
this area is rapidly growing, driven in part by its importance to engineering
[18, 19]. Given their relevance to large scale brain dynamics it is highly
desirable to develop mathematical techniques for the analysis of Wilson-Cowan
style neural mass models at the network level. Here we advocate for the
replacement of smooth sigmoidal nonlinearities in neural mass models by more
tractable functions, including piecewise linear (PWL) and piecewise constant
functions. A PWL continuous choice has been used in several previous studies,
including those of Hansel and Sompolinsky [20], and Kilpatrick and Bressloff
[21], whilst the discontinuous Heaviside (piecewise constant) choice has
proven especially popular since the seminal work of Amari [16]. In these
instances this has facilitated the construction of certain types of localised
states in continuum neural field models. However, in a discrete neural network
context there is a major mathematical difference in the analysis of network
states for the case of continuous vs. discontinuous firing rates. As well as
introducing a simple methodology to treat the construction of periodic orbits
in idealised Wilson-Cowan networks, this is one of the major topics we wish to
address in this paper.

First in §2 we introduce the model for an isolated Wilson-Cowan node with a
PWL firing rate. The description of dynamical states with reference to
switching manifolds becomes very useful. We show how matrix exponentials can
be used to patch together a periodic orbit, and that Floquet theory simplifies
considerably to yield explicit formulas for determining solution stability.
Next in §3 we consider a network of PWL Wilson-Cowan nodes, with nodes
arranged along a ring with distance-dependent interactions. This particular
choice of coupling guarantees the existence of the synchronous state. We then
develop a linear stability analysis of this state and show that this leads to
a tractable variational problem of a very similar type to that for the single
node, albeit now parameterised by the eigenvalues of the connectivity matrix.
We use this to determine instabilities that can lead to the formation of
spatio-temporal network patterns. Next in §4 we consider the case that the
firing rate is a Heaviside function, for which the techniques developed for
studying PWL systems break down. Once again periodic orbits can be constructed
using matrix exponentials, although standard Floquet theory must be now
augmented to cope with the evolution of linearised perturbations through the
switching manifolds. This is most readily achieved with the use of saltation
matrices that have commonly been used for the study of nonsmooth mechanical
systems [22]. However, at the network level the stability of the synchronous
state is much harder to determine than for the continuous model. Here we show
that ideas from the study of Glass networks developed by Edwards [23] are
particularly useful, and that stability is strongly influenced by the temporal
order in which network components cross switching manifolds, and that this in
turn is determined by the choice of initial perturbation. Finally in §6 we
conclude with an overview of the new results about synchrony in networks of
neural mass models, and discuss the natural extension of this work to treat
non-synchronous states.

##  2 The Wilson-Cowan model and a piecewise linear reduction

For their activity-based neural mass model Wilson and Cowan [1, 2]
distinguished between excitatory and inhibitory sub-populations. This seminal
(space-clamped) model can be written succinctly in terms of the pair of
coupled differential equations:

$\frac{{\rm d}u}{{\rm d}t}=-u+F(I_{u}+w^{uu}u-w^{vu}v),\qquad\tau\frac{{\rm
d}v}{{\rm d}t}=-v+F(I_{v}+w^{uv}u-w^{vv}v),$ (1)

Here $u=u(t)$ is a temporal coarse-grained variable describing the proportion
of excitatory cells firing per unit time at the instant $t$. Similarly the
variable $v$ represents the activity of an inhibitory population of cells. The
constants $w^{\alpha\beta}$, $\alpha,\beta\in\\{u,v\\}$, describe the weight
of all synapses from the $\alpha$th population to cells of the $\beta$th
population, and $\tau$ is a relative time-scale. The nonlinear function $F$
describes the expected proportion of neurons in population $\alpha$ receiving
at least threshold excitation per unit time, and is often taken to have a
sigmoidal form. Here the terms $I_{\alpha}$ represent external inputs (that
could be time varying). For a historical perspective on the Wilson-Cowan model
see [5], and for a more recent reflection by Cowan see [24]. To reduce the
model to a mathematically tractable form we consider the choice of a PWL
firing rate function given by

$F(x)=\begin{cases}0&x\leq 0\\\ \epsilon^{-1}x&0<x<\epsilon\\\
1&x\geq\epsilon\end{cases}.$ (2)

For appropriate choices of parameters the Wilson-Cowan model, with the firing
rate given by (2), can support stable oscillations. An example is shown in
Fig. 1, where we also plot the four switching manifolds defined by the
condition that arguments to the function $F$ in (1) take on the values zero
and $\epsilon$.

Figure 1: Phase plane for the Wilson-Cowan network with a PWL firing rate,
showing a stable periodic orbit (light blue). Parameters: $\epsilon=0.04$,
$\tau=0.6$, $I_{u}=-0.05$, $I_{v}=-0.3$, $w^{uu}=1$, $w^{vu}=2$, $w^{uv}=1$,
and $w^{vv}=0.25$. The straight lines in red and green show the switching
manifolds, where $I_{u}+w^{uu}u-w^{vu}v=0,\epsilon$ and
$I_{v}+w^{uv}u-w^{vv}v=0,\epsilon$ respectively.

Away from the switching manifolds the dynamics governing the evolution of
trajectories is linear, and may be constructed using matrix exponentials. To
simplify further analysis it is first convenient to introduce new variables
$(U,V)$ such that $u=(w^{vu}(V-I_{v})-w^{vv}(U-I_{u}))/|W|$, where $|W|=\det
W$, and $v=(w^{uu}(V-I_{v})-w^{uv}(U-I_{u}))/|W|$, as well as the matrices

$W=\begin{bmatrix}w^{uu}&-w^{vu}\\\ w^{uv}&-w^{vv}\end{bmatrix},\qquad
J=\begin{bmatrix}1&0\\\ 0&1/\tau\end{bmatrix}\qquad A=-WJW^{-1}.$ (3)

With these choices (1) transforms to

$\frac{{\rm d}}{{\rm d}t}\begin{bmatrix}U\\\
V\end{bmatrix}=A\begin{bmatrix}U-I_{u}\\\
V-I_{v}\end{bmatrix}+WJ\begin{bmatrix}F(U)\\\ F(V)\end{bmatrix}.$ (4)

In the representation (4) we see that the four switching manifolds are simply
defined by $U=0$, $U=\epsilon$, $V=0$, and $V=\epsilon$. The periodic orbit
shown in Fig. 1 (encircling an unstable fixed point) crosses each of these
manifolds twice, so that the periodic trajectory is naturally decomposed into
eight separate pieces. On each piece we shall denote the time-of-flight for a
trajectory to travel from one switching manifold to another by $\Delta_{i}$,
$i=1,\ldots,8$, so that the period of the orbit is given by
$\Delta=\sum_{i=1}^{8}\Delta_{i}$. As an explicit example of how to construct
a trajectory between two switching manifolds, consider the region where $0\leq
U\leq\epsilon$ and $V<0$. In this case the solution of (4) is given by

$\begin{bmatrix}U(t)\\\ V(t)\end{bmatrix}={\rm
e}^{A_{+}(\epsilon)t}\begin{bmatrix}U(0)\\\
V(0)\end{bmatrix}-{A}_{+}^{-1}(\epsilon)\left({\rm
e}^{A_{+}(\epsilon)t}-I_{2}\right)A\begin{bmatrix}I_{u}\\\
I_{v}\end{bmatrix},\qquad t\geq 0,$ (5)

where

$A_{+}(\epsilon)=\left(A+\epsilon^{-1}WJ\begin{bmatrix}1&0\\\
0&0\end{bmatrix}\right).$ (6)

It is a simple matter to write down the trajectories in each of the remaining
regions of phase space visited by a periodic orbit. We may then use these
matrix exponential formulas to patch together solutions, setting the origin of
time in each region such that initial data in one region comes from final data
from a trajectory in a neighbouring region. We shall denote the periodic orbit
by $(\overline{U},\overline{V})$ such that
$(\overline{U}(t),\overline{V}(t))=(\overline{U}(t+\Delta),\overline{V}(t+\Delta))$.
If we consider initial data with $(\overline{U}(0),\overline{V}(0))=(U_{0},0)$
then the eight times-of-flight and the unknown $U_{0}$ are determined self-
consistently by the nine equations $\overline{V}(\Delta_{1})=\epsilon$,
$\overline{U}(\Delta_{2})=\epsilon$, $\overline{U}(\Delta_{3})=0$,
$\overline{V}(\Delta_{4})=\epsilon$, $\overline{V}(\Delta_{5})=0$,
$\overline{U}(\Delta_{6})=0$, $\overline{U}(\Delta_{7})=\epsilon$,
$\overline{V}(\Delta_{8})=0$, and $\overline{U}(\Delta_{8})=U_{0}$. The
numerical solution of this nonlinear algebraic system of equations can be used
to construct periodic orbits such as the one shown in Fig. 1. Note that the
construction of periodic orbits that do not cross all of the switching
manifolds can similarly be performed (requiring the simultaneous solution of
fewer equations). To determine stability we can turn directly to Floquet
theory for planar systems which tells us that the non-trivial Floquet exponent
is given by

$\sigma=\frac{1}{\Delta}\int_{0}^{\Delta}\operatorname{Tr}D(s){\rm d}s,$ (7)

where $D(s)$ denotes the Jacobian of the system evaluated along the periodic
orbit. In general this is a hard quantity to evaluate for systems where the
periodic orbit is not available in closed form. However, for the PWL Wilson-
Cowan model the Jacobian is piecewise constant and we have that

$\sigma=\frac{1}{\Delta}\sum_{i=1}^{8}\Delta_{i}\operatorname{Tr}A_{i},$ (8)

where $A_{2}=A_{4}=A_{6}=A_{8}=A$, $A_{3}=A_{7}=A_{+}(\epsilon)$, and
$A_{1}=A_{5}=A_{-}(\epsilon)$, where

$A_{-}(\epsilon)=\left(A+\epsilon^{-1}WJ\begin{bmatrix}0&0\\\
0&1\end{bmatrix}\right).$ (9)

Thus a periodic orbit is stable if $\sigma<0$. In Fig. 2 we present a plot of
$\sigma$ as a function of $\tau$, to show that the periodic solution in Fig. 1
is stable.

Figure 2: A plot of the non-trivial Floquet exponent for the PWL Wilson-Cowan
model (left axis), as a function of the relative time-scale $\tau$, with the
period of the orbit also shown (right axis). Parameters as in Fig. 1. Periodic
orbits emerge via a supercritical Hopf bifurcation as $\tau$ increases through
$\tau_{\text{Hopf}}=(w^{vv}+\epsilon)/(w^{uu}-\epsilon)\sim 0.3$. We see that
the branch of periodic orbits shown is stable, with stability decreasing to
zero as the solution is lost with increasing $\tau$. This loss of existence
occurs because of a grazing bifurcation (coincident with a saddle-node
bifurcation of periodic orbits) at $\tau_{\text{graze}}\sim 0.6$ whereby part
of the trajectory develops a point of inflection on the switching manifold
$v=(I_{u}+w^{uu}u)/w^{vu}$ (red solid line in Fig. 1), such that beyond
bifurcation the trajectory does not cross the switching manifold and instead
is attracted to the stable fixed point at $(u,v)=(0,0)$.

Given the above method to construct and determine the stability of a periodic
orbit, we next show how to extend this aproach to treat synchronous solutions
in networks of Wilson-Cowan oscillators.

##  3 A piecewise linear Wilson-Cowan network

The study of coupled oscillator networks in biology, physics, and engineering
is now commonplace. Two particularly well known tools for studying patterns of
phase-locked states and their instabilities are the theory of weakly coupled
oscillators [25], and the master stability function (MSF) [26]. The reduction
of a coupled limit cycle network to a phase oscillator network has proven very
useful for gaining insight into phenomena ranging from the synchronisation of
flashing fireflies [27] to behaviours in social networks [28], and for a
recent review see [29]. However, there is an obvious limitation to such an
approach, namely the restriction to weak interaction (and near identical
oscillators). The MSF approach (for identical oscillators) does not require
any such restriction on coupling strength, and can be used to determine the
stability of the synchronous state in terms of the eigen-structure of the
network connectivity matrix. However, the numerical evolution of a system of
dynamical equations, arising from a Floquet variational problem, must be
performed. Importantly the MSF approach can be combined with group theoretical
techniques used in the study of symmetric dynamical systems to analyse the
stability of cluster states within symmetric networks of dynamical units [30,
31]. Here we favour the MSF approach and show it simplifies considerably for a
PWL choice of firing rate function. This allows us to improve upon previous
mathematical studies of Wilson-Cowan networks, such as those by Campbell and
Wang [32] (who treated networks with nearest neighbour coupling and
established the condition for synchrony), Ueta and Chen [33] (who performed a
numerical bifurcation analysis for small networks), and Ahmadizadeh et al.
[34] (who used perturbation techniques and numerics to study synchrony in
networks with diffusive coupling).

We now consider a network of Wilson-Cowan nodes given by

$\displaystyle\frac{{\rm d}u_{i}}{{\rm d}t}$
$\displaystyle=-u_{i}+F\left(I_{u}+\sum_{j=1}^{N}\mathcal{W}^{uu}_{ij}u_{j}-\sum_{j=1}^{N}\mathcal{W}^{vu}_{ij}v_{j}\right),$
(10) $\displaystyle\tau\frac{{\rm d}v_{i}}{{\rm d}t}$
$\displaystyle=-v_{i}+F\left(I_{v}+\sum_{j=1}^{N}\mathcal{W}^{uv}_{ij}u_{j}-\sum_{j=1}^{N}\mathcal{W}^{vv}_{ij}v_{j}\right),\qquad
i=1,\ldots,N,$ (11)

subject to the constraints $\sum_{j=1}^{N}\mathcal{W}^{uu}_{ij}=w^{uu}$,
$\sum_{j=1}^{N}\mathcal{W}^{vu}_{ij}=w^{vu}$,
$\sum_{j=1}^{N}\mathcal{W}^{uv}_{ij}=w^{uv}$, and
$\sum_{j=1}^{N}\mathcal{W}^{vv}_{ij}=w^{vv}$ for all $i$. These row-sum
constraints are natural for networks arranged on a ring, and guarantee the
existence of a synchronous orbit $(u_{i}(t),v_{i}(t))=(u(t),v(t))$ for all
$i=1,\ldots,N$, where $(u(t),v(t))$ is given by the solution of (1).

It is now convenient to introduce a vector notation with
$X=(u_{1},v_{1},u_{2},v_{2},\ldots,u_{N},v_{N})\in\RSet^{2N}$ and consider a
change of variables $Y=\mathcal{W}X+C$, where
$C=\mathbf{1}_{N}\otimes(I_{u},I_{v})$, and

$\mathcal{W}=\mathcal{W}^{uu}\otimes\begin{bmatrix}1&0\\\
0&0\end{bmatrix}-\mathcal{W}^{vu}\otimes\begin{bmatrix}0&1\\\
0&0\end{bmatrix}+\mathcal{W}^{uv}\otimes\begin{bmatrix}0&0\\\
1&0\end{bmatrix}-\mathcal{W}^{vv}\otimes\begin{bmatrix}0&0\\\
0&1\end{bmatrix}.$ (12)

Here the symbol $\otimes$ denotes the usual tensor product for matrices, and
$\mathbf{1}_{N}$ is an $N$-dimensional vector with all entries equal to unity.
This means that the switching manifolds can be succinctly described by
$Y_{i}=0$ and $Y_{i}=\epsilon$, and the dynamics takes the form

$\frac{{\rm d}}{{\rm d}t}Y=\mathcal{A}(Y-C)+\mathcal{W}\mathcal{J}F(Y),$ (13)

where

$\mathcal{J}=I_{N}\otimes
J,\qquad\mathcal{A}=-\mathcal{W}\mathcal{J}\mathcal{W}^{-1},$ (14)

where $J$ is given by (3) and $I_{N}$ is the $N\times N$ identity matrix. If
we denote the synchronous solution by
$\overline{Y}(t)=(\overline{U}(t),\overline{V}(t),\overline{U}(t),\overline{V}(t),\ldots,\overline{U}(t),\overline{V}(t))$
and consider small perturbations such that $Y=\overline{Y}+\delta Y$, then
these evolve according to

$\frac{{\rm d}}{{\rm d}t}\delta Y=\mathcal{A}\delta
Y+\mathcal{W}\mathcal{J}DF(\overline{Y})\delta Y,$ (15)

where $DF(\overline{Y})$ is the Jacobian of $F$ evaluated along the periodic
orbit.

Given the constraints on the matrices $\mathcal{W}^{\alpha\beta}$, with
$\alpha,\beta\in\\{u,v\\}$ it is natural to take these to be circulant
matrices with
$\mathcal{W}^{\alpha\beta}_{ij}=\mathcal{W}^{\alpha\beta}_{|i-j|}$. In this
case the normalised eigenvectors of $\mathcal{W}^{\alpha\beta}$ are given by
$e_{p}=(1,\omega_{p},\omega_{p}^{2},\ldots,\omega_{p}^{N-1})/\sqrt{N}$, where
$p=0,\ldots,N-1$, and $\omega_{p}=\exp(2\pi ip/N)$ are the $N$th roots of
unity. The corresponding complex eigenvalues are given by
$\nu^{\alpha\beta}=\nu^{\alpha\beta}(p)$ where

$\nu^{\alpha\beta}(p)=\sum_{\mu=0}^{N-1}\mathcal{W}^{\alpha\beta}_{\mu}\omega_{p}^{\mu}.$
(16)

If we introduce the matrix of eigenvectors
$P=[e_{0}~{}e_{1}~{}\ldots~{}e_{N-1}]$, then we have that

$\displaystyle(P\otimes I_{2})^{-1}\mathcal{W}(P\otimes I_{2})$
$\displaystyle=\Lambda^{uu}\otimes\begin{bmatrix}1&0\\\
0&0\end{bmatrix}-\Lambda^{vu}\otimes\begin{bmatrix}0&1\\\
0&0\end{bmatrix}+\Lambda^{uv}\otimes\begin{bmatrix}0&0\\\
1&0\end{bmatrix}-\Lambda^{vv}\otimes\begin{bmatrix}0&0\\\ 0&1\end{bmatrix},$
$\displaystyle=\operatorname{diag}(\Lambda(0),\Lambda(1),\ldots,\Lambda(N-1))\equiv\Lambda,$
(17)

where
$\Lambda^{\alpha\beta}=\operatorname{diag}(\nu^{\alpha\beta}(0),\nu^{\alpha\beta}(1),\ldots,\nu^{\alpha\beta}(N-1))$,
and

$\Lambda(p)=\begin{bmatrix}\nu^{uu}(p)&-\nu^{vu}(p)\\\
\nu^{uv}(p)&-\nu^{vv}(p)\end{bmatrix},\qquad p=0,1,\ldots,N-1.$ (18)

Moreover, it is easy to establish that in the above notation $(P\otimes
I_{2})^{-1}\mathcal{A}(P\otimes I_{2})=-\Lambda(I_{N}\otimes J)\Lambda^{-1}$.

If we now consider perturbations of the form $\delta Z=(P\otimes
I_{2})^{-1}\delta Y$ then from (15) we find that the linearised dynamics is
described by the system

$\frac{{\rm d}}{{\rm d}t}\delta Z=\Lambda(I_{N}\otimes
J)\left[-\Lambda^{-1}+(I_{N}\otimes D)\right]\delta Z,$ (19)

where $D\in\RSet^{2\times 2}$ is the Jacobian of
$(F(\overline{U}),F(\overline{V}))$, and is a piecewise constant matrix that
is only non-zero if $0<\overline{U}(t)<\epsilon$ or
$0<\overline{V}(t)<\epsilon$. In the former case $[DF]_{11}=\epsilon^{-1}$
with all other entries zero, and in the latter case $[DF]_{22}=\epsilon^{-1}$
with all other entries zero. We see that (19) has a block structure where the
dynamics in each of $N$ $2\times 2$ blocks is given by

$\frac{{\rm d}}{{\rm d}t}\xi=[A(p)+\Lambda(p)JD]\xi,\qquad
p=0,\ldots,N-1,\qquad\xi\in\RSet^{2},$ (20)

with $A(p)=-\Lambda(p)J\Lambda^{-1}(p)$. Thus, comparing to (4), we see that
the variational equation for the network is identical to that for a single
Wilson-Cowan unit with $W$ replaced by $\Lambda(p)$. We note that for $p=0$
the variational problem is identical to that for an isolated node since
$\Lambda(0)=W$ (using
$\nu^{\alpha\beta}(0)=\sum_{\mu=0}^{N-1}\mathcal{W}_{\mu}^{\alpha\beta}=w^{\alpha\beta}$).
Thus to determine the stability of the synchronous state we only have to
consider a set of $N$ two dimensional variational problems. Exploiting the
fact that between switching manifolds the variational problem defined by (20)
is time-independent we may construct a solution in a piecewise fashion from
matrix exponentials and write $\xi(t)=\exp[(A(p)+\Lambda(p)JD)t]\xi(0)$. We
may then build up a perturbed trajectory over one period of oscillation in the
form $\xi(\Delta)=\Gamma(p)\xi(0)$, where $\Gamma(p)\in\RSet^{2\times 2}$ is
given by

$\Gamma(p)={\rm e}^{A(p)\Delta_{8}}{\rm e}^{A_{+}(p;\epsilon)\Delta_{7}}{\rm
e}^{A(p)\Delta_{6}}{\rm e}^{A_{-}(p;\epsilon)\Delta_{5}}{\rm
e}^{A(p)\Delta_{4}}{\rm e}^{A_{+}(p;\epsilon)\Delta_{3}}{\rm
e}^{A(p)\Delta_{2}}{\rm e}^{A_{-}(p;\epsilon)\Delta_{1}},$ (21)

where

$A_{+}(p;\epsilon)=\left(A(p)+\epsilon^{-1}\Lambda(p)J\begin{bmatrix}1&0\\\
0&0\end{bmatrix}\right),\quad
A_{-}(p;\epsilon)=\left(A(p)+\epsilon^{-1}\Lambda(p)J\begin{bmatrix}0&0\\\
0&1\end{bmatrix}\right).$ (22)

Thus if a periodic orbit of an isolated Wilson-Cowan node is stable then the
synchronous network solution will be stable provided all the eigenvalues of
$\Gamma(p)$, for $p=0,\ldots,N-1$, lie in the unit disc (excluding the one
that arises from time-translation invariance, with a value $+1$). For a fixed
value of $p$ one of three bifurcations is possible, namely a tangent
instability defined by $\det(\Gamma(p)-I_{2})=0$, a period-doubling
instability defined by $\det(\Gamma(p)+I_{2})=0$, and a Neimark-Sacker
bifurcation defined by $\det\Gamma(p)=1$. If there is a $p=p_{c}$ such that
one of these instabilities occurs then the excited network state will
correspond to the eigenvector $\text{Re}\,e_{p_{c}}$.

###  3.1 Example: a ring network

By way of illustration of the above theory let us consider a network of
Wilson-Cowan nodes arranged on a ring with an odd number of nodes. Introducing
a distance between nodes indexed by $i$ and $j$ as
$\operatorname{dist}(i,j)=\min(|i-j|,N-|i-j|)$, we can define a set of
exponentially decaying connectivity matrices according to

$\mathcal{W}^{\alpha\beta}_{ij}=w^{\alpha\beta}\frac{{\rm
e}^{-\operatorname{dist}(i,j)/\sigma_{\alpha\beta}}}{\sum_{j=0}^{N-1}{\rm
e}^{-\operatorname{dist}(0,j)/\sigma_{\alpha\beta}}}.$ (23)

Thus we have a set of four circulant matrices parametrised by the four spatial
scales $\sigma_{\alpha\beta}$ that respect the row-sum constraints
$\sum_{j=1}^{N}\mathcal{W}^{\alpha\beta}_{ij}=w^{\alpha\beta}$. In Fig. 3 we
show a plot of the eigenvalues of $\Gamma(p)$ for $p=0,\ldots,N-1$ for two
different parameter choices. In one case all of the eigenvalues (excluding the
one arising from time-translation invariance) lie within the unit disc, whilst
in the other one leaves the unit disc along the negative real axis.

Figure 3: Spectral plots in the complex plane for a Wilson-Cowan ring network,
with spatial scales $\sigma_{\alpha\beta}=\sigma$ for all $\alpha,\beta$, and
$N=31$. Other parameters as in Fig. 1. Left: $\sigma=0.15$, and the
synchronous solution is predicted to be linearly stable. Right:
$\sigma=0.191$, and the synchronous solution is predicted to be linearly
unstable.

This latter scenario predicts an instability of the synchronous state, and is
consistent with direct numerical simulations. Moreover, by studying the
spectrum under parameter variation we can find the value of $p=p_{c}$ which
goes unstable first. In Fig. 4 we show time courses (obtained by direct
numerical simulation) for the components $u_{i}(t)$ of the emergent network
state just beyond the point of instability, as well as a plot of the real part
of the spatial eigenvector $e_{p_{c}}$. We see that the spatial pattern of the
network state is well predicted by $e_{p_{c}}$, suggesting that the
bifurcation is supercritical.

Figure 4: Direct numerical simulation of a Wilson-Cowan ring network, with
$N=31$, just beyond the point of synchronous instability where $\sigma=0.191$.
Other parameters as in Fig. 1. Here we plot the components $u_{i}(t)$ in a
space-time plot. The shape of the unstable mode $e_{p_{c}}$, with $p_{c}=16$
(and also $p_{c}=17$ because of a degeneracy) is depicted in blue at the top
of the figure. The bifurcation point of the linear instability is found to be
in excellent agreement with simulations, with the spatial pattern of the
emergent network state predicted by $e_{p_{c}}$.

##  4 The Heaviside world

In a recent paper Harris and Ermentrout [35] considered a single Wilson-Cowan
population with a Heaviside nonlinearity, where the firing rate in (1) takes
the form $F(x)=H(x)$, where $H(x)=0$ for $x<0$ and $H(x)=1$ for $x>0$. The
choice of a Heaviside firing rate has been very popular in mathematical
neuroscience ever since the seminal work of Amari (for neural field models),
as nicely exemplified by his recent article on the “Heaviside World” [36]. A
case in point is the work of Laing and Chow [37] for understanding binocular
rivalry. They considered a neural mass network model with recurrent
excitation, cross-inhibition, adaptation, and synaptic depression and showed
that the use of a Heaviside nonlinearity allowed the explicit calculation of
the dominance durations of perceptions. A more recent use of the Heaviside
firing rate has been by McCleney and Kilpatrick [38] for neural activity
models with spike rate adaptation to understand the dynamics of up-down
states. Using techniques from Filippov systems and differential inclusions
Harris and Ermentrout made a study of periodic orbits for a Heaviside firing
rate using a boundary value problem approach. Here we show that we can recover
their results using the matrix exponential approach of §2. Moreover, we also
extend their work on a single node by showing how to determine the stability
of periodic orbits using a nonsmooth version of Floquet theory.

In the representation (4), with $F=H$, we see that the there are two switching
manifolds defined by $U=0$ and $V=0$. If we introduce the indicator functions
$h_{1}(U,V)=U$ and $h_{2}(U,V)=V$ then we can define these manifolds (lines in
this case) as

$\Sigma_{i}=\left\\{(U,V)\in\RSet^{2}~{}|~{}h_{i}(U,V)=0\right\\}.$ (24)

These switching manifolds naturally divide the plane into four sets. We denote
these by $D_{++}=\\{(U,V)~{}|U\geq 0,V\geq 0\\}$, $D_{+-}=\\{(U,V)~{}|U\geq
0,V\leq 0\\}$, $D_{--}=\\{(U,V)~{}|U\leq 0,V\leq 0\\}$, and
$D_{-+}=\\{(U,V)~{}|U\leq 0,V\geq 0\\}$. If we denote the elements of $A$ by
$A_{ij}$, $i=1,2$ and $j=1,2$, where

$A=-\frac{1}{|W|}\begin{bmatrix}w^{vu}w^{uv}/\tau-w^{uu}w^{vv}&w^{uu}w^{vu}(1-1/\tau)\\\
w^{vv}w^{uv}(1/\tau-1)&w^{uv}w^{vu}-w^{uu}w^{vv}/\tau\end{bmatrix},\qquad|W|=w^{vu}w^{uv}-w^{uu}w^{vv},$
(25)

then the $U$-nullclines are given by

$V=I_{v}-\frac{A_{11}(U-I_{u})}{A_{12}}+\frac{1}{A_{12}}\begin{cases}-w^{uu}+w^{vu}/\tau&(U,V)\in
D_{++}\\\ -w^{uu}&(U,V)\in D_{+-}\\\ 0&(U,V)\in D_{--}\\\ w^{vu}/\tau&(U,V)\in
D_{-+}\end{cases},$ (26)

and the $V$-nullclines are given by

$V=I_{v}-\frac{A_{21}(U-I_{u})}{A_{22}}+\frac{1}{A_{22}}\begin{cases}-w^{uv}+w^{vv}/\tau&(U,V)\in
D_{++}\\\ -w^{uv}&(U,V)\in D_{+-}\\\ 0&(U,V)\in D_{--}\\\ w^{vv}/\tau&(U,V)\in
D_{-+}\end{cases}.$ (27)

An example set of nullclines is shown in Fig. 5.

Figure 5: Phase plane for a Wilson-Cowan node with a Heaviside firing rate
(transformed coordinates), showing the $U$-nullclines (red) and $V$-nullclines
(green), as well as a stable periodic orbit (blue), and an unstable periodic
sliding orbit (dashed magenta). Parameters (excluding $\epsilon$) as in Fig.
1.

To discuss fixed points and their stability it is first necessary to complete
the description of the dynamics on the switching manifolds. We do this using
the Filippov convex method [39] and extend our discontinuous system into a
convex differential inclusion. The Filippov extension of (4) is then

$\frac{{\rm d}}{{\rm d}t}\begin{bmatrix}U\\\ V\end{bmatrix}\in
F(U,V)=\begin{cases}F_{++}(U,V)&(U,V)\in D_{++}\\\
\overline{\text{co}}\left(\\{F_{++},F_{+-}\\},\kappa_{1}\right)&(U,V)\in
D_{++}\cap D_{+-}\\\ F_{+-}(U,V)&(U,V)\in D_{+-}\\\
\overline{\text{co}}\left(\\{F_{+-},F_{--}\\},\kappa_{2}\right)&(U,V)\in
D_{+-}\cap D_{--}\\\ F_{--}(U,V)&(U,V)\in D_{--}\\\
\overline{\text{co}}\left(\\{F_{--},F_{-+}\\},\kappa_{3}\right)&(U,V)\in
D_{--}\cap D_{-+}\\\ F_{-+}(U,V)&(U,V)\in D_{-+}\\\
\overline{\text{co}}\left(\\{F_{-+},F_{++}\\},\kappa_{4}\right)&(U,V)\in
D_{-+}\cap D_{++}\end{cases},$ (28)

where $F_{\alpha\beta}(U,V)=A[U-I_{u},V-I_{v}]^{T}+b_{\alpha\beta}$ for
$\alpha,\beta\in\\{+,-\\}$ and

$b_{++}=\begin{bmatrix}w^{uu}-w^{vu}/\tau\\\
w^{uv}-w^{vv}/\tau\end{bmatrix},\quad b_{+-}=\begin{bmatrix}w^{uu}\\\
w^{uv}\end{bmatrix},\quad b_{--}=\begin{bmatrix}0\\\ 0\end{bmatrix},\quad
b_{-+}=\begin{bmatrix}-w^{vu}/\tau\\\ -w^{vv}/\tau\end{bmatrix}.$ (29)

Here $\overline{\text{co}}(\\{f,g\\},\kappa)=\kappa f+(1-\kappa)g$ with
$\kappa\in[0,1]$ is the closed convex hull of all values between $f$ and $g$.
A sliding solution may exist along a switching manifold such that
$\dot{h}_{i}=\nabla h_{i}\cdot F=0$. The functions $\kappa_{j}$,
$j=1,\ldots,4$, are chosen to ensure that $\dot{h}_{i}=0$ along any switching
manifold. For example if a sliding solution exists along the line $U=0$ for
$V<0$ then we would construct $\kappa_{2}$ using $\nabla h_{1}=(1,0)$ and
$F(0,V)=\kappa_{2}F_{+-}(0,V)+(1-\kappa_{2})F_{--}(0,V)$ yielding

$\kappa_{2}=\frac{(1,0)\cdot
F_{--}(0,V)}{(1,0)\cdot(F_{--}(0,V)-F_{+-}(0,V))}.$ (30)

As illustrated in Fig. 5 it is possible for two nullclines to intersect and
create a fixed point $(U_{ss},V_{ss})$. In the example shown this occurs for
$U<0$ and $V<0$, so that $(U_{ss},V_{ss})=(I_{u},I_{v})$. Linear stability
analysis shows that this is a stable node (with eigenvalues of $A$, namely
$-1$ and $-1/\tau$). Moreover, this system also supports pseudo equilibria
where either a nullcline touches a switching manifold, or two switching
manifolds intersect. A thorough exploration of the pseudo equilibria of (1)
can be found in [35]. Here we shall simply focus on the pseudo equilibrium at
$(U_{ss},V_{ss})=(0,0)$, and characterise its stability by considering
trajectories around this point. In fact given the PWL nature of the dynamics
it is sensible to consider the construction of periodic orbits, and determine
the stability of the pseudo equilibrium in terms of the stability of
encircling small amplitude orbits.

###  4.1 Periodic orbits and their stability

A non-sliding periodic orbit around $(0,0)$ can be constructed in terms of the
times-of-flight in each region $D_{\alpha\beta}$. If we denote these four
times by the symbols $\Delta_{\alpha\beta}$ then the period of the orbit is
given by $\Delta=\Delta_{++}+\Delta_{-+}+\Delta_{--}+\Delta_{+-}$. We may then
use a matrix exponential solution:

$\begin{bmatrix}U(t)\\\ V(t)\end{bmatrix}={\rm e}^{At}\begin{bmatrix}U(0)\\\
V(0)\end{bmatrix}+(I_{2}-{\rm e}^{At})\left[\begin{bmatrix}I_{u}\\\
I_{v}\end{bmatrix}-A^{-1}WJ\begin{bmatrix}H(U)\\\
H(V)\end{bmatrix}\right],\qquad t\geq 0.$ (31)

to patch together solutions, setting the origin of time in each region such
that initial data in one region comes from final data from a trajectory in a
neighbouring region. We shall denote the periodic orbit by
$(\overline{U},\overline{V})$ such that
$(\overline{U}(t),\overline{V}(t))=(\overline{U}(t+\Delta),\overline{V}(t+\Delta))$.To
indicate which region we are considering we shall simply add $\alpha\beta$
subscripts to the formula in (31). In this way a periodic orbit that visits
all four regions in turn can be parameterised by the five unknowns
$\overline{U}_{++}(0)$, $\overline{V}_{++}(\Delta_{++})$,
$\overline{U}_{-+}(\Delta_{-+})$, $\overline{V}_{--}(\Delta_{--})$,
$\overline{U}_{+-}(\Delta_{+-})$, and $\Delta_{\alpha\beta}$. These are
determined self-consistently by the five equations
$\overline{U}_{++}(\Delta_{++})=0$, $\overline{V}_{-+}(\Delta_{-+})=0$,
$\overline{U}_{--}(\Delta_{--})=0$, $\overline{V}_{+-}(\Delta_{+-})=0$, and
$\overline{U}_{+-}(\Delta_{+-})=\overline{U}_{++}(0)$. To determine the
stability of such an orbit we may use the nonsmooth Floquet theory described
in [40]. In essence this treats the propagation of perturbations through a
switching manifold using a saltation matrix, such that
$Y(T^{+})=\lim_{\epsilon\searrow 0}Y(T+\epsilon)=KY(T)$, where $Y=(U,V)$
denotes the vector state of the system and $K\in\RSet^{2\times 2}$ is the
saltation matrix that acts at time $T$. Saltation matrices can be derived in a
number of ways, with a general prescription in terms of an indicator function
$h$ as [22]

$K=I_{2}+\frac{\left[\dot{Y}(T^{+})-\dot{Y}(T)\right]\left[\vphantom{\dot{Y}(T)}\nabla_{Y}h(Y(T))\right]^{\mathsf{T}}}{\nabla_{Y}h(Y(T))\cdot\dot{Y}(T)}.$
(32)

Alternatively, in the context of the PWL model discussed in §2, we can obtain
the relevant saltation matrices by considering the approximation
$H(x)=\lim_{\epsilon\rightarrow 0}F(x)$. To see this we introduce the vector
$\overline{Y}(t)=(\overline{U}(t),\overline{V}(t))$ and linearise the
equations of motion (4) by considering $Y(t)=\overline{Y}(t)+\delta Y(t)$, for
small perturbations $\delta Y(t)=(\delta U,\delta V)$. The linearised
equations of motion are given by

$\frac{{\rm d}}{{\rm d}t}\delta Y=\left[A+WJDF(\overline{Y}(t))\right]\delta
Y.$ (33)

Here $DF(\overline{Y}(t))$ is the piecewise constant matrix described after
(19). Consider for example the time of flight, $t_{1}(\epsilon)$, between
$U=\epsilon$ and $U=0$. For small $\epsilon$ we may estimate $t_{1}(\epsilon)$
using the result that $U(t)\simeq
U(t_{0})+\left.\dot{U}\right|_{t=t_{0}}(t-t_{0})$, giving
$t_{1}(\epsilon)=-\epsilon/\left.\dot{U}\right|_{t=\Delta_{++}}$. The
corresponding change in state across this small time interval can be obtained
by integrating (33) to give

$\delta Y(T^{+})-\delta Y(T)=\lim_{\epsilon\rightarrow
0}\int_{T}^{T+t_{1}(\epsilon)}WJ\begin{bmatrix}\epsilon^{-1}&0\\\
0&0\end{bmatrix}\delta Y(t){\rm d}t.$ (34)

Thus we obtain $\delta Y(T^{+})=K_{1}\delta Y^{-}$, with the saltation matrix
$K_{1}$ given by

$K_{1}=I_{2}-\frac{1}{\dot{U}(\Delta_{++})}WJ\begin{bmatrix}1&0\\\
0&0\end{bmatrix}.$ (35)

The other saltation matrices (describing the passage through
$\epsilon$-neighbourhoods of $U=0$ and $V=0$) are constructed in a similar
fashion, and found to be

$\displaystyle K_{2}$
$\displaystyle=I_{2}-\frac{1}{\dot{V}(\Delta_{-+})}WJ\begin{bmatrix}0&0\\\
0&1\end{bmatrix},$ $\displaystyle K_{3}$
$\displaystyle=I_{2}+\frac{1}{\dot{U}(\Delta_{--})}WJ\begin{bmatrix}1&0\\\
0&0\end{bmatrix},$ $\displaystyle K_{4}$
$\displaystyle=I_{2}+\frac{1}{\dot{V}(\Delta_{+-})}WJ\begin{bmatrix}0&0\\\
0&1\end{bmatrix}.$ (36)

It is straightforward to check that the saltation matrices (35)-(36) are
equivalent to those defined by (32). We now pass to the limit $\epsilon=0$, to
treat the Heaviside firing rate. Between switching events the perturbations
evolve according to $\exp(A(t-T))\delta Y(T^{+})$, for $t>T$, where $\delta
Y(T^{+})$ is the perturbation at the switching time. Thus after one period of
oscillation we may put this all together to obtain

$\delta Y(\Delta)=\Gamma\delta Y(0),\qquad\Gamma=K_{4}{\rm
e}^{A\Delta_{+-}}K_{3}{\rm e}^{A\Delta_{--}}K_{2}{\rm
e}^{A\Delta_{-+}}K_{1}{\rm e}^{A\Delta_{++}}.$ (37)

The periodic orbit will be stable if the eigenvalues of $\Gamma$ lie within
the unit disc. Note that one of the Floquet multipliers is equal to one,
corresponding to perturbations along the periodic orbit. Let us denote the
other eigenvalue by ${\rm e}^{\sigma\Delta}$ and use the result that
$\det\Gamma={\rm e}^{\sigma\Delta}\times 1$. Hence,

$\displaystyle{\rm e}^{\sigma\Delta}$ $\displaystyle=\left(\prod_{i=1}^{4}\det
K_{i}\right)\det{\rm e}^{A\Delta_{+-}}\det{\rm e}^{A\Delta_{--}}\det{\rm
e}^{A\Delta_{-+}}\det{\rm e}^{A\Delta_{++}}$
$\displaystyle=\frac{\dot{\overline{V}}(\Delta_{+-}^{+})}{\dot{\overline{V}}(\Delta_{+-})}\frac{\dot{\overline{U}}(\Delta_{--}^{+})}{\dot{\overline{U}}(\Delta_{--})}\frac{\dot{\overline{V}}(\Delta_{-+}^{+})}{\dot{\overline{V}}(\Delta_{-+})}\frac{\dot{\overline{U}}(\Delta_{++}^{+})}{\dot{\overline{U}}(\Delta_{++})}\det{\rm
e}^{A\Delta_{+-}}\det{\rm e}^{A\Delta_{--}}\det{\rm e}^{A\Delta_{-+}}\det{\rm
e}^{A\Delta_{++}}.$ (38)

Using the fact that $\det{\rm e}^{At}={\rm e}^{\operatorname{Tr}A\,t}$ we find

$\sigma=-\left(1+\frac{1}{\tau}\right)+\frac{1}{\Delta}\log\frac{\dot{\overline{V}}(\Delta_{+-}^{+})}{\dot{\overline{V}}(\Delta_{+-})}\frac{\dot{\overline{U}}(\Delta_{--}^{+})}{\dot{\overline{U}}(\Delta_{--})}\frac{\dot{\overline{V}}(\Delta_{-+}^{+})}{\dot{\overline{V}}(\Delta_{-+})}\frac{\dot{\overline{U}}(\Delta_{++}^{+})}{\dot{\overline{U}}(\Delta_{++})}.$
(39)

A periodic orbit will be stable provided $\sigma<0$. We shall say that the
pseudo-equilibria at $(0,0)$ is unstable (stable) if it is enclosed by a
stable (unstable) periodic orbit of arbitrarily small amplitude. We shall say
that there is a pseudo-Hopf bifurcation at $(0,0)$ when the pseudo-equilibrium
changes stability, namely when $\sigma=0$. A plot of $\sigma=\sigma(\tau)$
(not shown) for the parameters of Fig. 2, shows very similar behaviour as for
the steep PWL firing rate function. In essence we may regard the second term
on the right hand side of (39) as a correction term to standard Floquet theory
to cope with the nonsmooth nature of the Heaviside firing rate.

###  4.2 An unstable periodic sliding orbit

The Wilson-Cowan node can also support an unstable periodic orbit that has a
component which slides along the switching manifold $U=0$ for
$V\in[V_{1},V_{2}]$, as depicted in Fig. 5. The points $V_{1,2}$ are easily
calculated by determining the points at which the $U$-nullclines touch the
switching manifold where $U=0$, and are found to be
$V_{1}=(A_{11}I_{u}+A_{12}I_{V}-w^{uu})/A_{12}$ and
$V_{2}=V_{1}+w^{uu}/A_{12}$. In reverse time initial data close to a sliding
trajectory would be attracted to it. Thus we can think of constructing an
unstable periodoc sliding orbit, of the type shown in Fig. 5, by breaking it
into five pieces. All pieces of this orbit are constructed similarly to before
(see above), except the component that slides. Using the Filippov method and
equation (30) we find $\kappa_{2}=(A_{11}I_{u}-A_{12}V+A_{12}I_{v})/w^{uu}$,
with the sliding dynamics prescribed by

$\frac{{\rm d}}{{\rm d}t}\begin{bmatrix}U\\\
V\end{bmatrix}=\begin{bmatrix}0&0\\\
0&A_{22}-A_{11}w^{uv}/w^{uu}\end{bmatrix}\begin{bmatrix}U\\\
V\end{bmatrix}+\begin{bmatrix}0\\\ b_{s}\end{bmatrix},$ (40)

where $b_{s}=-A_{12}I_{u}-A_{22}I_{v}+(A_{11}I_{u}+A_{12}I_{v})w^{uv}/w^{uu}$.
In backward time the periodic sliding orbit shown in Fig. 5 would slide up
along $U=0$ until the point $V=V_{2}$, where it would leave the switching
manifold.

We now turn our attention to networks built from Wilson-Cowan nodes with a
Heaviside firing rate.

##  5 A network of Heaviside Wilson-Cowan nodes

As we have shown in §4 the replacement of a sigmoidal firing rate by a
Heaviside function can lead to highly tractable models for which substantial
analytical results can be obtained (with the use of matrix exponentials and
saltation matrices). However, at the network level the mathematical
differences between the treatment of smooth and nonsmooth firing rates are
considerably amplified relative to those at the single node level. At the node
level it is well known that regarding the Heaviside function as the steep
limit of a sigmoidal function can lead to arbitrarily many different non-
equivalent dynamical systems. This is simply due to the non-uniqueness of the
singular limits by which smooth functions may tend towards discontinuities.
For a recent perspective on this issue see the work of Jeffrey [41]. Thus
there is no reason to assume that taking the limit $\epsilon\rightarrow 0$ for
the PWL network considered in §3 will be relevant to a Wilson-Cowan network
with a Heaviside nonlinearity. Namely the approximation of a Heaviside
function by a continuous function such that $H(x)=\lim_{\epsilon\rightarrow
0}F(x)$, where $F(x)$ is given by (2), may have little utility given that
pointwise convergence need not imply distributional convergence.

We now return to the network introduced in §3, but replace the dynamics of
each node with the Heaviside limit studied in the previous section. For the
following analysis, it is convenient to rewrite (13) as

$\frac{{\rm d}}{{\rm
d}t}Y=\mathcal{A}(Y-\mathcal{F}(Y))\,,\qquad\mathcal{F}(Y)=C-\mathcal{A}^{-1}\mathcal{W}\mathcal{J}H(Y).$
(41)

The network model (13), with a Heaviside nonlinearity, is reminiscent of a so-
called Glass network originally introduced for the study of biochemical
networks that are dominated by switch-like behaviour [42, 43], though here the
model has two-time scales. For a nice survey of periodic and aperiodic
behaviour in Glass networks we recommend the article by Edwards [23], and for
the application to gene networks see Edwards and Glass [44].

The synchronous network state is given by (31) (remembering the row-sum
constraint on the network connections). To study its linear stability we
consider values of the perturbed network state $Y$ that are close to the
synchronous network state at the unperturbed crossing times. Let
$\overline{T}_{i}$ denote the time that the synchronous state moves between
one of the four quadrants (as illustated in Fig. 5). We then make the ansatz
that the perturbed network state $Y$ can be expressed with respect to the
synchronous orbit at one of the switching times $\overline{T}_{i}$ and write
$Y(t)=\overline{Y}(\overline{T}_{i})+\delta Y(t)$ with $t$ in the
neighbourhood of $\overline{T}_{i}$.

We first construct the saltation matrix through a switch, indexed by
$i=1,\ldots,4$. Suppose that the $k$th crossing occurs at a perturbed crossing
time $T_{i,k}$. The network states at two consecutive crossings are related
via

$Y(T_{i,k+1})={\rm
e}^{\mathcal{A}(T_{i,k+1}-T_{i,k})}Y(T_{i,k})+\left(I_{2N}-{\rm
e}^{\mathcal{A}(T_{i,k+1}-T_{i,k})}\right)\mathcal{F}(Y(T_{i,k}^{+}))\,.$ (42)

This equation is obtained by integrating (41) using the observation that
$\mathcal{F}$ is constant between crossings. By linearising (42) we can relate
the perturbations between crossing events as

$\delta Y(T_{i,k+1})=\delta Y(T_{i,k})+Y^{i,k}\delta T_{i,k}\,,$ (43)

where
$Y^{i,k}=\mathcal{A}(\overline{Y}(\overline{T}_{i})-\mathcal{F}(Y(T_{i,k}^{+})))$
and $\delta T_{i,k}=T_{i,k+1}-T_{i,k}$. For the node that crosses at
$T_{i,k+1}$, the corresponding component of $\delta Y(T_{i,k+1})$, say at
position $m$, vanishes, since
$Y_{m}(T_{i,k+1})=\overline{Y}_{m}(\overline{T}_{i})$ (namely the $m$th
component of the perturbed trajectory equals the $m$th component of the
synchronous orbit). Here, $m\in\\{1,3,\ldots,2N-1\\}$ or $m\in\\{2,4,\ldots
2N\\}$, depending on whether the crossing occurs along the $V$ or $U$ axis. We
then see from (43) that

$\delta T_{i,k}=-\frac{\delta Y_{m}(T_{i,k})}{Y^{i,k}_{m}}\,.$ (44)

At this point, $m$ is still unknown. However, since $m$ corresponds to the
node that crosses before any of the other remaining nodes do so, we find it by
minimising (44) over the possible values of $m$, and we denote it by $m_{k}$.
When we combine (43) and (44), we find that $\delta
Y(T_{i,k+1})=\Gamma_{i,k}\delta Y(T_{i,k})$ with

$\Gamma_{i,k}=I_{2N}-\frac{Y^{i,k}e_{m_{k}}^{\mathsf{T}}}{Y^{i,k}_{m_{k}}}\\\
,$ (45)

where $e_{m}$ is the $m$th canonical basis vector in $\mathbb{R}^{2N}$. The
saltation matrix for each of the four switches is then given by

$L_{i}=\Gamma_{i,N-1}\Gamma_{i,N-2}\cdots\Gamma_{i,1}\,,\qquad i=1,\ldots,4.$
(46)

The ordering of matrix multiplications in (46) is determined by the iterative
minimisation of the perturbations given by (44).

In the next step, we analyse how a perturbed network state is propagated
between saltation events. Let $T_{i}^{+}$ denote the time when the last node
crosses between quadrants. Here, the superscript makes explicit that all nodes
have crossed into the next quadrant. The next network event occurs when one of
the nodes crosses into the subsequent quadrant. This happens at a time
$T_{i+1}^{-}$, where the superscript indicates that only one node has crossed.
We will make the ansatz that $T_{i}^{+}=\overline{T}_{i}+\delta T_{i}^{+}$ and
$T_{i+1}^{-}=\overline{T}_{i+1}+\delta T_{i+1}^{-}$. We see from (41) that

$Y(T_{i+1}^{-})={\rm
e}^{\mathcal{A}(T_{i+1}^{-}-T_{i}^{+})}Y(T_{i}^{+})+\left(I_{2N}-{\rm
e}^{\mathcal{A}(T_{i+1}^{-}-T_{i}^{+})}\right)\mathcal{F}(Y(T_{i}^{+}))\,,$
(47)

from which we obtain after linearisation

$\delta Y(T_{i+1}^{-})={\rm e}^{\mathcal{A}\Delta_{i}}\left(\delta
Y(T_{i}^{+})-\overline{Y}^{\prime}(\overline{T}_{i}^{+})\delta
T_{i}^{+}\right)+\overline{Y}^{\prime}(\overline{T}_{i+1}^{-})\delta
T_{i+1}^{-}\,,\\\ $ (48)

where we have used the fact that
$\mathcal{F}(Y(T_{i}^{+}))=\mathcal{F}(Y(T_{i+1}^{-}))$, since $\mathcal{F}$
is constant between crossing events. Here $\overline{Y}^{\prime}(t)$ denotes
the differential of $\overline{Y}(t)$ with respect to $t$. As above, the
component of $\delta Y(T_{i+1}^{-})$ that corresponds to the node that
switches first, say at position $m$, vanishes. Taking the $m$th component of
(48) then yields an expression for the perturbation of the crossing time

$\delta
T_{i+1}^{-}=-\frac{f_{m}^{i}}{\overline{Y}_{m}^{\prime}(\overline{T}_{i+1}^{-})}\,,$
(49)

where the vector $f^{i}\in\RSet^{2N}$ is given by ${\rm
e}^{\mathcal{A}\Delta_{i}}\left(\delta
Y(T_{i}^{+})-\overline{Y}^{\prime}(\overline{T}_{i}^{+})\delta
T_{i}^{+}\right)$. We again find the value of $m$ by minimising (49) over all
admissible values of $m$ and refer to it as $m_{i}$. This leads to $\delta
Y(T_{i+1}^{-})=\Gamma_{i}\delta Y(T_{i}^{+})$ with

$\Gamma_{i}=\left(G_{i}-\frac{\overline{Y}^{\prime}(\overline{T}_{i+1}^{-})}{\overline{Y}_{m_{i}}^{\prime}(\overline{T}_{i+1}^{-})}e_{m_{i}}^{T}G_{i}\right)\,,\qquad$
(50)

and

$G_{i}={\rm
e}^{\mathcal{A}\Delta_{i}}\left(I_{2N}-\frac{\overline{Y}^{\prime}(\overline{T}_{i}^{+})e_{1}^{T}\delta
T_{i}^{+}}{\delta Y_{1}(T_{i}^{+})}\right)\,.$ (51)

Taken together, we obtain after one period

$\delta Y(T_{4}^{+})=\Psi\delta
Y(0)\,,\qquad\Psi=L_{4}\Gamma_{4}L_{3}\Gamma_{3}L_{2}\Gamma_{2}L_{1}\Gamma_{1}\,.$
(52)

The matrices $\Gamma_{i}$ act to propagate perturbations across a quadrant,
and the $L_{i}$ propagate perturbations through a switch. At first sight, the
definition of $G_{i}$ suggests that we have introduced a dependence of
$\Gamma_{i}$ on $\delta Y(0)$ through the inclusion of $\delta Y(T_{i}^{+})$.
This dependence can be avoided by noting that $\delta T_{i}^{+}=\delta
T_{i}^{-}+\sum_{k}\delta T_{i,k}$ and the repeated use of (44), (45) and (49).
The drawback of this approach is that the resultant operator does not lend
itself to an interpretation of successive propagations and saltations, nor is
it numerically advantageous. Moreover, this operator would only remove the
_explicit_ dependence of $\Psi$ on $\delta Y(0)$. The minimisation steps that
are necessary to determine the order in which nodes switch already leads to an
_implicit_ dependence of $\Psi$ on $\delta Y(0)$. Changing $\delta Y(0)$ can
lead to a different order of switching, and since matrix multiplication does
not commute, $\Psi$ can be different for different $\delta Y(0)$. This has
profound implications for asserting linear stability. The usual argument that
the eigenvalues of $\Psi$ determine linear stability does not hold anymore. To
see this, consider the propagation of $\delta Y(0)$ over multiple periods,
i.e.

$\delta Y^{(1)}=\Psi^{(0)}\delta Y(0)\,,\quad\delta Y^{(2)}=\Psi^{(1)}\delta
Y^{(1)}\,,\quad\delta Y^{(3)}=\Psi^{(2)}\delta Y^{(2)}\,,\quad\ldots$ (53)

so that

$\delta Y^{(m)}=\Psi^{(m-1)}\Psi^{(m-2)}\cdots\Psi^{(0)}\delta Y(0)\,.$ (54)

The eigenvalues of $\Psi^{(i)}$ and $\Psi^{(j)}$ can be different for $i\neq
j$. For some value of $i$ $\Psi^{(i)}$ may have all eigenvalues in the unit
disc, whilst for another value of $i$ there may be some eigenvalues outside
the unit disc. Over one period, perturbations can therefore grow or shrink.
This entails that for a product of operators as in (54), $\delta Y^{(m)}$ may
be smaller than $\delta Y^{(0)}$, although some $\Psi^{(i)}$ might have some
eigenvalues that lie outside the unit disc. Instead of looking at the
eigenvalues of individual $\Psi^{(i)}$, we could have studied the eigenvalues
of the product of operators in (54). We would have come to the same conclusion
since eigenvalues of the product operator move into and out of the unit disc
as we increase $m$.

Figures 6 and 7 illustrate the dependence of the spectra on random initial
conditions $\delta Y(0)$. In both figures, the left panel shows the spectra
for initial conditions when all eigenvalues of $\Psi^{(0)}$ lie within the
unit disc. The middle panel displays spectra with some eigenvalues outside the
unit disc, and the right panel is a blowup of the middle panel around the unit
disc. For Fig. 6, we chose a value of $\sigma$ such that the synchronous orbit
of the PWL network, with a small values of $\epsilon=0.001$, is linearly
stable. We observe that the eigenvalues of the Heaviside network cluster
around those of the PWL network. While it appears that the majority of
synchronous solutions are stable (for this parameter choice), some initial
conditions lead to eigenvalues outside the unit disc. When zooming into the
unit disc, we see some degree of clustering, although this is not as
pronounced as for the stable solutions.

Figure 6: Spectral plots for a Heaviside Wilson-Cowan ring network with
spatial scales $\sigma_{\alpha\beta}=0.215$ for all $\alpha,\beta$, and $N=5$.
We sampled $2000$ random initial conditions, and eigenvalues are shown as open
red circles. The filled blue circles are the eigenvalues of the PWL network
with the same parameter values and $\epsilon=0.001$. (Left) Spectra for
initial conditions that lead to eigenvalues that all fall into the unit disc.
(Middle) Spectra for initial conditions that lead to eigenvalues outside the
unit disc. (Right) Blow-up of the the middle panel around the unit disc. Other
parameters as in Fig. 1.

For larger values of $\sigma$, the synchronous state of the PWL network
becomes unstable (for small $\epsilon)$. The left panel of Fig. 7 shows that
the eigenvalues of the Heaviside network that all fall into the unit disc
exhibit only a weak association with the eigenvalues of the PWL network. In
addition, it seems that more initial conditions lead to unstable synchronous
solutions than stable ones. This mirrors the behaviour in Fig. 6, where the
majority of initial conditions gives rise to stable solutions. The blow-up in
the right panel of Fig. 7 illustrates that the eigenvalues of the Heaviside
network form clusters around those of the PWL network.

Figure 7: Spectral plots for a Heaviside Wilson-Cowan ring network with
spatial scales $\sigma_{\alpha\beta}=0.23$ for all $\alpha,\beta$, and $N=5$.
We sampled $2000$ random initial conditions, and eigenvalues are shown as open
red circles. The filled blue circles are the eigenvalues of the PWL network
with the same parameter values with $\epsilon=0.001$. (Left) Spectra for
initial conditions that lead to eigenvalues that all fall into the unit disc.
(Middle) Spectra for initial conditions that lead to eigenvalues outside the
unit disc. (Right) Blow-up of the the middle panel around the unit disc. Other
parameters as in Fig. 1.

While the notion of linear stability in terms of eigenvalues of the propagator
is lost for the Heaviside network, it appears that the clustering of these
eigenvalues reflects the stability of the PWL system, at least for small
values of $\epsilon$ (where the PWL firing rate becomes more switch like).

##  6 Conclusion

In this paper we have shown that the combination of two popular approaches in
dynamical systems, namely PWL modelling of low dimensional oscillators and the
MSF, can be combined to give insight into the behaviour of network states in
neural mass network models. This is natural for this type of system since the
sigmoidal nonlinearity, ubiquitous throughout neuroscience modelling of large
scale brain dynamics, is well caricatured by a PWL reduction. We have focused
here on the bifurcation of the synchronous network state, and shown how this
can be determined in terms of a set of low-dimensional Floquet problems, each
of which can be solved using simple linear algebra. In essence the PWL aspect
of the model allows the variational problem for stability to be solved without
recourse to the numerical solution of an ordinary differential equation.
Closed form solutions are patched together, and although this may appear
inelegant at first sight, it does lead to explicit formulas for Floquet
exponents at the single node level, and is easily cast into algorithmic form
for accurate numerical computations at the network level. This nicely
highlights the benefits of PWL modelling. Importantly the approach advocated
here is not just limited to the construction and stability of the synchronous
state. Pecora et al. [30] and Sorrentino et al. [31] have recently extended
the MSF approach to treat more exotic states making extensive use of tools
from computational group theory. Thus the work presented here is readily
extended to treat non-synchronous states, such as clusters, and for a further
discussion see [45]. From a neuroscience perspective it would also be
important to treat delays, arising from the finite propagation speed of action
potentials relaying signals between distinct brain regions [46]. In this case
we would hope to exploit the growing body of knowledge on PWL dynamics with
time delay, as exemplified by [47].

From a mathematical perspective we have also seen that there is an important
difference between the analysis of a high gain continuous PWL sigmoid and that
of a discontinuous switch-like Heaviside firing rate. Although this can be
facilitated with the use of saltation matrices (to propagate perturbations
through switching manifolds) there is no MSF style approach that reduces the
study of synchrony to a set of sub-network Floquet problems. Moreover, in
contrast to the linear stability analysis of continuous systems, there is now
a new challenge of addressing the temporal order in which perturbations to
network states pass through a switching manifold. To treat this we have made
use of ideas originally developed for Glass networks [23], though note that
similar issues of ordering also arise in the analysis of pulse-coupled systems
[48, 49, 50]. In essence the analysis of a Wilson-Cowan network with a
Heaviside firing rate must be performed carefully, and with non-standard
tools, as its behaviour can differ from that of a similar network with a high
gain PWL sigmoid.

###### Acknowledgements.

This work was supported by the Engineering and Physical Sciences Research
Council [grant number EP/P007031/1]. M. Şayli was supported by a grant from
TÜBITAK.

## References

  * [1] H R Wilson and J D Cowan.  Excitatory and inhibitory interactions in localized populations of model neurons.  Biophysical Journal, 12:1–24, 1972.
  * [2] H R Wilson and J D Cowan.  A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue.  Kybernetik, 13:55–80, 1973.
  * [3] J D Cowan.  History of concepts and techniques.  Intelligent Systems, 3:375–400, 2004.
  * [4] S Coombes, P beim Graben, R Potthast, and J Wright, editors.  Neural Fields: Theory and Applications.  Springer, 2014.
  * [5] A Destexhe T J Sejnowski.  The Wilson–Cowan model, 36 years later.  Biological Cybernetics, 101:1–2, 2009.
  * [6] M W Woolrich and K E Stephan.  Biophysical network models and the human connectome.  NeuroImage, 80:330–338, 2013.
  * [7] A C E Onslow, M W Jones, and R Bogacz.  A canonical circuit for generating phase-amplitude coupling.  PLoS ONE, 9:1–15, 2014.
  * [8] H G E Meijer, T L Eissa, B Kiewiet, J F Neuman, C A Schevon, R G Emerson, R R Goodman, G M McKhann, C J Marcuccilli, A K Tryba, J D Cowan, S A van Gils, and W van Drongelen.  Modeling focal epileptic activity in the Wilson–Cowan model with depolarization block.  The Journal of Mathematical Neuroscience, 5:7, 2015.
  * [9] C A Lea-Carnall, M A Montemurro, N J Trujillo-Barreto, L M Parkes, and W El-Deredy.  Cortical resonance frequencies emerge from network size and connectivity.  PLoS Computational Biology, 12:1–19, 2016.
  * [10] P A Valdes-Sosa, J M Sanchez-Bornot, R C Sotero, Y Iturria-Medina, Y Aleman-Gomez, J Bosch-Bayard, F Carbonell, and T Ozaki.  Model driven EEG/fMRI fusion of brain oscillations.  Human Brain Mapping, 30:2701–21, 2009.
  * [11] L H Zetterberg, L Kristiansson, and K Mossberg.  Performance of a model for a local neuron population.  Biological Cybernetics, 31:15–26, 1978.
  * [12] B H Jansen and V G Rit.  Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns.  Biological Cybernetics, 73:357–366, 1995.
  * [13] D T J Liley, P J Cadusch, and M P Dafilis.  A spatially continuous mean field theory of electrocortical activity.  Network, 13:67–113, 2002.
  * [14] P Sanz-Leon, S A Knock, A Spiegler, and V K Jirsa.  Mathematical framework for large-scale brain network modeling in The Virtual Brain.  NeuroImage, 111:385–430, 2015.
  * [15] M Breakspear.  Dynamic models of large-scale brain activity.  Nature Neuroscience, 20:340–352, 2017.
  * [16] S Amari.  Dynamics of pattern formation in lateral-inhibition type neural fields.  Biological Cybernetics, 27:77–87, 1977.
  * [17] S Coombes.  Waves, bumps, and patterns in neural field theories.  Biological Cybernetics, 93:91–108, 2005.
  * [18] M di Bernardo, C Budd, A R Champneys, and P Kowalczyk.  Piecewise-smooth Dynamical Systems: Theory and Applications.  Applied Mathematical Sciences. Springer, 2008.
  * [19] A Colombo, M di Bernardo, S J Hogan, and M R Jeffrey.  Bifurcations of piecewise smooth flows: Perspectives, methodologies and open problems.  Physica D, 241:1845–1860, 2012.  Dynamics and Bifurcations of Nonsmooth Systems.
  * [20] D Hansel and H Sompolinsky.  Methods in Neuronal Modeling, chapter Modeling Feauture Selectivity in Local Cortical Circuits.  MIT Press, second edition, 1998.
  * [21] Z P Kilpatrick and P C Bressloff.  Effects of synaptic depression and adaptation on spatiotemporal dynamics of an excitatory neuronal network.  Physica D, 239:547–560, 2010.
  * [22] R Leine and H Nijmeijer.  Dynamics and Bifurcations of Non-Smooth Mechanical Systems, volume 18 of Lecture Notes in Applied and Computational Mechanics book series.  Springer-Verlag, Berlin, 2004.
  * [23] R Edwards.  Analysis of continuus-time switching networks.  Physica D, 146:165–199, 2000.
  * [24] J D Cowan, J Neuman, and W van Drongelen.  Wilson-Cowan equations for neocortical dynamics.  The Journal of Mathematical Neuroscience, 6(1), 2016.
  * [25] Y Kuramoto.  Chemical Oscillations, Waves and Turbulence.  Springer-Verlag, Heidelberg, 1984.
  * [26] L M Pecora and T L Carroll.  Master stability functions for synchronized coupled systems.  Physical Review Letters, 80:2109–2112, 1998.
  * [27] G B Ermentrout.  An adaptive model for synchrony in the firefly Pteroptyx malaccae.  Journal of Mathematical Biology, 29:571–585, 1991.
  * [28] S Assenza, R Gutiérrez, J Gómez-Gardenes, V Latora, and S Boccaletti.  Emergence of structural patterns out of synchronization in networks with competitive interactions.  Scientific Reports, 1, 2011.
  * [29] F Dörfler and F Bullo.  Synchronization in complex networks of phase oscillators: A survey.  Automatica, 50:1539–1564, 2014.
  * [30] L M Pecora, F Sorrentino, A M Hagerstrom, T E Murphy, and R Roy.  Cluster synchronization and isolated desynchronization in complex networks with symmetries.  Nature Communications, 5, 2014.
  * [31] F Sorrentino, L M Pecora, A M Hagerstrom, T E Murphy, and R Roy.  Complete characterization of the stability of cluster synchronization in complex dynamical networks.  Science Advances, 2:e1501737–e1501737, 2016.
  * [32] S Campbell and D Wang.  Synchronization and desynchronization in a network of locally coupled Wilson-Cowan oscillators.  IEEE Transactions on Neural Networks, 7:541–554, 1996.
  * [33] T Ueta and G Chen.  On synchronization and control of coupled Wilson Cowan neural oscillators.  International Journal of Bifurcation and Chaos, 13:163–175, 2003\.
  * [34] S Ahmadizadeh, D Nesic, D B Grayden, and D R Freestone.  Analytic synchronization conditions for a network of Wilson and Cowan oscillators.  IEEE 54th Annual Conference on Decision and Control, 2015.
  * [35] J Harris and B Ermentrout.  Bifurcations in the Wilson–Cowan Equations with Nonsmooth Firing Rate.  SIAM Journal on Applied Dynamical Systems, 14:43–72, 2015.
  * [36] S-I Amari.  Neural Field Theory, chapter Heaviside World: Excitation and Self-Organization of Neural Fields.  Springer, 2013.
  * [37] C R Laing and C C Chow.  A spiking neuron model for binocular rivalry.  Journal of Computational Neuroscience, 12:39–53, 2002.
  * [38] Z T McCleney and Z P Kilpatrick.  Entrainment in up and down states of neural populations: non-smooth and stochastic models.  Journal of Mathematical Biology, 73:1131–1160, 2016.
  * [39] A F Filippov.  Differential Equations with Discontinuous Righthand Sides.  Kluwer Academic Publishers, Norwell, 1988.
  * [40] S Coombes and R Thul.  Synchrony in networks of coupled non-smooth dynamical systems: Extending the master stability function.  European Journal of Applied Mathematics, 27:904–922, 2016.
  * [41] M R Jeffrey.  Smoothing tautologies, hidden dynamics, and sigmoid asymptotics for piecewise smooth systems.  Chaos, 25(10):103125, 2015.
  * [42] L Glass and S A Kauffman.  the logical analysis of continuous non-linear biochemical control networks.  Journal of Theoretical Biology, 39:103–129, 1971.
  * [43] L Glass.  Classification of biological networks by their qualitative dynamics.  Journal of Theoretical Biology, 54:85–107, 1975.
  * [44] R Edwards and L Glass.  Combinatorial explosion in model gene networks.  Chaos, 10:691–704, 2000.
  * [45] R Nicks, L Chambon, and S Coombes.  Clusters in nonsmooth oscillator networks.  Physical Review X, submitted, 2017.
  * [46] G Deco, V K Jirsa, A R McIntosh, O Sporns, and R Kötter.  Key role of coupling, delay, and noise in resting brain fluctuations.  Proceedings of the National Academy of Sciences, 106:10302–10307, 2009.
  * [47] D V Senthilkumar and M. M Lakshmanan.  Bifurcations and chaos in time delayed piecewise linear dynamical systems.  International Journal of Bifurcation and Chaos, 15:2895–2912, 2005\.
  * [48] M Timme, F Wolf, and T Geisel.  Coexistence of regular and irregular dynamics in complex networks of pulse-coupled oscillators.  Physical Review Letters, 89:258701, 2002.
  * [49] H Kielblock, C Kirst, and M Timme.  Breakdown of order preservation in symmetric oscillator networks with pulse-coupling.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 21:025113, 2011.
  * [50] P Goel and B Ermentrout.  Synchrony, stability, and firing patterns in pulse-coupled oscillators.  Physica D, 163:191 – 216, 2002.