# Paradoxical phase response of gamma rhythms facilitates their entrainment in
heterogeneous networks

Xize Xu, Hermann Riecke

###### Abstract

The synchronization of different $\gamma$-rhythms arising in different brain
areas has been implicated in various cognitive functions. Here, we focus on
the effect of the ubiquitous neuronal heterogeneity on the synchronization of
PING (pyramidal-interneuronal network gamma) and ING (interneuronal network
gamma) rhythms. The synchronization properties of rhythms depends on the
response of their collective phase to external input. We therefore determined
the macroscopic phase-response curve for finite-amplitude perturbations
(fmPRC), using numerical simulation of all-to-all coupled networks of
integrate-and-fire (IF) neurons exhibiting either PING or ING rhythms. We show
that the intrinsic neuronal heterogeneity can qualitatively modify the fmPRC.
While the phase-response curve for the individual IF-neurons is strictly
positive (type I), the fmPRC can be biphasic and exhibit both signs (type II).
Thus, for PING rhythms, an external excitation to the excitatory cells can, in
fact, delay the collective oscillation of the network, even though the same
excitation would lead to an advance when applied to uncoupled neurons. This
paradoxical delay arises when the external excitation modifies the internal
dynamics of the network by causing additional spikes of inhibitory neurons,
whose delaying within-network inhibition outweighs the immediate advance
caused by the external excitation. These results explain how intrinsic
heterogeneity allows the PING rhythm to become synchronized with a periodic
forcing or another PING rhythm for a wider range in the mismatch of their
frequencies. We demonstrate a similar mechanism for the synchronization of ING
rhythms. Our results identify a potential function of neuronal heterogeneity
in the synchronization of coupled $\gamma$-rhythms, which may play a role in
neural information transfer via communication through coherence.

Department of Engineering Science and Applied Mathematics,
Northwestern University, Evanston, IL 60208, USA

Author Summary

The interaction of a large number of oscillating units can lead to the
emergence of a collective, macroscopic oscillation in which many units
oscillate in near-unison or near-synchrony. This has been exploited
technologically, e.g., to combine many coherently interacting, individual
lasers to form a single powerful laser. Collective oscillations are also
important in biology. For instance, the circadian rhythm of animals is
controlled by the near-synchronous dynamics of a large number of individually
oscillating cells. In animals and humans brain rhythms reflect the coherent
dynamics of a large number of neurons and are surmised to play an important
role in the communication between different brain areas. To be functionally
relevant, these rhythms have to respond to external inputs and have to be able
to synchronize with each other. We show that the ubiquitous heterogeneity in
the properties of the individual neurons in a network can contribute to that
ability. It can allow the external inputs to modify the internal network
dynamics such that the network can follow these inputs over a wider range of
frequencies. Paradoxically, while an external perturbation may delay
individual neurons, their ensuing within-network interaction can
overcompensate this delay, leading to an overall advance of the rhythm.

##  1 Introduction

Collective oscillations or rhythms representing the coherent dynamics of a
large number of coupled oscillators play a significant role in many systems.
In the technological realm they range from laser arrays and Josephson
junctions to micromechanical oscillators [1, 2]. Among the important
biological examples are the heart rhythm, the circadian rhythm generated by
the suprachiasmatic nucleus [3], the segmentation clock controlling the somite
formation during development [4], and brain waves [5]. One prominent brain
rhythm is the widely observed $\gamma$-rhythm with frequencies in the range
30-100Hz. The coherent spiking of the neurons underlying this rhythm likely
enhances the downstream impact of the neurons participating in the rhythm. The
rhythmic alternation of low and high activity has been suggested to play a
significant role in the communication between different brain areas [6, 7].
That communication has also been proposed to be controled by the coherence of
the rhythms in the participating brain areas [8, 9, 10, 11, 12, 13].

For collective oscillations or rhythms to play a constructive role in a system
they need to respond adequately to external perturbations and stimuli. For
instance, for the circadian rhythm it is essential that it can be reliably
entrained by light and phase-lock to its daily variation. Similarly, if
rhythms are to play a significant role in the communication between different
brain areas, their response to input from other areas represents a significant
determinant of their function. Moreover, the stimulation and entrainment of
$\gamma$-rhythms by periodic sensory input is being considered as a
therapeutic approach for some neurodegenerative diseases [14].

Even small perturbations can affect oscillations significantly in that they
can advance or delay the oscillations, i.e. they can change the phase of the
oscillators. This change typically depends not only on the strength of the
perturbation but, importantly, also on the timing of the perturbations and is
expressed in terms of the phase response curve (PRC), which has been studied
extensively for individual oscillators [15]. For infinitesimal perturbations
the PRC can be determined elegantly using the adjoint method [16].

If the collective oscillation of a network of interacting oscillators is
sufficiently coherent, that system can be thought of as a single effective
oscillator. Consequently, the response of the macroscopic phase of the
collective oscillation to external perturbations and the mutual interaction of
multiple collective oscillations is of interest. The macroscopic phase-
response curve (mPRC) has been obtained in various configurations, including
noise-less heterogeneous phase oscillators [17, 18], noisy identical phase
oscillators [19, 20], noisy excitable elements [21], and noisy oscillators
described by the theta-model [22], which is equivalent to the quadratic
integrate-fire model for spiking neurons. Recent work has used the reduction
of networks of quadratic integrate-fire neurons to two coupled differential
equations for the firing rate and the mean voltage [23], which is related to
the Ott-Antonsen theory [24, 25], to develop a method to obtain the
infinitesimal macroscopic PRC (imPRC) for excitatory-inhibitory spiking
networks [26, 27].

A key difference between the response of an individual oscillator to a
perturbation and that of a collective oscillation is the fact that the degree
of synchrony of the collective oscillation can change as a result of the
perturbation, reflecting a change in the relations between the individual
oscillators. Thus, the phase response of a collective oscillation to a brief
perturbation consists not only of the immediate change in the phases of the
individual oscillators caused by the perturbation, but includes also a change
in the collective phase that can result from the subsequent convergence back
to the phase relationship between the oscillators corresponding to the
synchronized state, which is likely to have been changed by the perturbation
[18]. Interestingly, it has been observed that the infinitesimal macroscopic
phase response can be qualitatively different from the phase response of the
individual elements. Thus, even if the individual oscillators have a type-I
PRC, i.e. a PRC that is strictly positive or negative, the mPRC of the
collective oscillation can be of type II, i.e. it can exhibit a sign change as
a function of the phase [21, 22, 28].

Here we investigate the interplay between external perturbations and the
internal interactions among neurons in inhibitory and in excitatory-inhibitory
networks exhibiting $\gamma$-rhythms of the ING- and of the PING-type. We
focus on networks comprised of neurons that are not identical, leading to a
spread in their individual phases and a reduction in the degree of their
synchrony. How does this phase dispersion affect the response of the
macroscopic phase of the rhythm to perturbations? Does it modify the ability
of the network to follow a periodic perturbation ?

We show that the dispersion in the phase together with the within-network
interactions among the neurons can be the cause of a paradoxical phase
response: an external perturbation that _delays_ each individual neuron can
_advance_ the macroscopic rhythm. We identify the following mechanism
underlying this paradoxical response: external perturbations that delay
individual neurons sufficiently allow the within-network inhibition generated
by early-spiking neurons to suppress the spiking of less excited neurons. This
results in a reduced within-network inhibition, which reduces the time to the
next spike volley, speeding up the rhythm. This paradoxical phase response
increases with the neuronal heterogeneity and allows the network to phase-lock
to periodic external perturbations over a wider range of detuning. Thus, the
desynchronization within the network enhances its synchronizabilty with other
networks. The mechanism is closely related to that underlying the enhancement
of synchronization of collective oscillations by uncorrelated noise [29] and
the enhanced entrainment of the rhythm of a homogeneous network to periodic
input if that input exhibits phase dispersion across the network [30, 31]. We
demonstrate and analyze these behaviors for networks of inhibitory neurons
(ING-rhythm) and for networks comprised of excitatory and inhibitory neurons
(PING-rhythm).

##  2 Results

We investigated the impact of neuronal heterogeneity on the response of the
phase of $\gamma$-rhythms to brief external perturbations and the resulting
ability of rhythms to synchronize to periodic input. As described in the
Methods, we used networks comprised of minimal integrate-fire neurons that
interact with each other through synaptic pulses modeled via delayed double-
exponentials. To study ING-rhythms all neurons were inhibitory, while for the
PING-rhythms we used excitatory-inhibitory networks. In both cases, the
coupling within each population was all-to-all. Throughout, we implemented the
neuronal heterogeneity by injecting a different steady bias current $I_{bias}$
into each neuron. Our analysis suggests that the origin of the neuronal
heterogeneity plays only a minor role as long as it leads to a dispersion of
their spike times [29].

### Paradoxical Phase Response of Heterogeneous Networks: ING-Rhythm

Figure 1: ING-rhythm can be advanced by inhibition while individual neurons
are delayed. (A) Top: spike raster of neurons spiking sequentially in the
order of their input strength (increasing with neuron index). Bottom: mean
voltage across the network (LFP). (B) External inhibition advanced the rhythm.
Top: raster plot of spikes without (black) and with (red) external inhibitory
pulse. Bottom: Average of the total inhibitory current each neuron received
from the other neurons within the network.  $I^{(I)}=$ 20.4 pA, $C_{V}^{(I)}$=
0.15, $f_{network}=47$ Hz. In (B), perturbations were made with a square-wave
inhibitory current pulse with duration 0.1 ms and amplitude 3200 pA to each
neuron, resulting in a 4 mV rapid hyperpolarization.

In the absence of external perturbations the all-to-all inhibition among the
neurons lead to rhythmic firing of the neurons. Due to their heterogeneity
they did not spike synchronously but sequentially, as shown in Fig.1A, where
the neurons are ordered by the strength of their bias current. The dependence
of the phase dispersion on the coefficient of variation of the heterogeneity
in the bias current (CV) is shown in Suppl. Figure S1. For sufficiently large
heterogeneity some neurons never spiked: while the weak bias current they
received would have been sufficient to induce a spike eventually, the strong
inhibition that was generated by the neurons spiking earlier in the cycle
suppressed those late spikes. Neurons with strong bias current could spike
multiple times.

A brief, inhibitory external input delivered to all neurons (green dashed line
in Fig.1B) delayed each neuron. The degree of this individual delay depended
on the timing of the input, as is reflected in the PRC of the individual
neurons. If the perturbation was applied during the time between the spike
volleys, the delay of each neuron had no further consequence and the overall
rhythm was delayed. However, if the same inhibitory perturbation arrived
during a spike volley (dashed green line in Fig.1B), it could advance the
overall rhythm. As illustrated in Fig.1B, only the spiking of the late neurons
was delayed by the perturbation. Importantly, with this delay some neurons did
not spike before the within-network inhibition triggered by the early-spiking
neurons (dashed blue line in Fig.1B) became strong enough to suppress the
spiking of the late neurons altogether. With fewer neurons spiking, the all-
to-all inhibition within the network was reduced, allowing all neurons to
recover earlier, which lead to a shorter time to the next spike volley. If the
speed-up was larger than the immediate delay induced by the external
inhibition, the overall phase of the rhythm was advanced by the delaying
inhibition.

Figure 2: fmPRC of heterogeneous ING network. (A) Phase shift in response to
inhibition for different neuronal heterogeneity but fixed natural frequency
($f_{network}=40$Hz). The paradoxical phase advance increased with neuronal
heterogeneity. (B) fmPRC changed qualitatively with the amplitude of the
perturbation. Left: fmPRC for three different perturbation amplitudes. Right:
raster plot of spikes without (black) and with (red) external inhibition. Top:
strong perturbation advanced the network. Bottom: weak perturbation applied at
the same time as in the top figure. The network was delayed. (C) Maximal phase
advance as a function of neuronal heterogeneity and external inhibition
strength. The threshold of the inhibition amplitude to obtain an advance
decreased with heterogeneity (white line). $f_{network}$ was kept constant
($f_{network}=40$Hz). In (A)-(C), perturbations were made with a square-wave
inhibitory current pulse with duration 0.1 ms to each interneuron. In (A), the
amplitude of the current was 1600 pA, resulting in a 2 mV rapid
hyperpolarization.

As the example in Fig.1B shows, the paradoxical phase response requires proper
timing of the perturbation. We therefore determined quantitatively the
macroscopic phase-response curve (PRC) of the rhythm. To do so we measured
computationally the amount a brief current injection shifted the phase of the
rhythm (Fig.2A). We defined the phase as the normalized time since the first
spike in the most recent volley of spikes. Reflecting the strictly positive
PRC of the individual integrate-fire neurons, without heterogeneity ($CV=0$)
external inhibition always delayed the rhythm, independent of the timing of
the pulse. In contrast, in heterogeneous networks the rhythm could be advanced
if the same inhibitory perturbation was applied shortly after the first spikes
in the spike volley ($\phi_{inh}>0$). Increasing the neuronal heterogeneity
enhanced this phase advance, since it shifted the within-network inhibition
driven by the leading neurons to earlier times, while it delayed the lagging
neurons. As a result, for the same external perturbation, a larger fraction of
neurons that would spike in the absence of the external inhibition was
sufficiently delayed to have their spikes be suppressed by the within-network
inhibition (cf. Fig.1B), reducing the within-network inhibition and with it
the time to the next spike volley. To keep the frequency of the unperturbed
network fixed in Fig.2A, we reduced the tonic input with increasing
heterogeneity, which enhanced the phase advance. However, even if the tonic
input was kept fixed, the phase advance increased with heterogeneity (Suppl.
Fig.S2).

For weak heterogeneity the paradoxical phase response occurred only for
sufficiently strong perturbations, i.e. it did not arise in the infinitesimal
macroscopic PRC (imPRC). Thus, the phase response changed qualitatively as the
amplitude of the perturbation was strong enough to delay the spikes of
sufficiently many slow neurons until the self-inhibition of the network set in
and suppressed their spikes (Fig.2B). As the CV of the neurons was increased,
the dispersion was large enough that the spikes of the lagging neurons were
suppressed by the self-inhibition of the network even in the absence of an
external perturbation. Above that threshold value of CV the paradoxical phase
response occurred even for infinitesimal perturbations (Fig.2C).

Figure 3: The paradoxical phase response of a heterogeneous ING network is
robust. (A) The phase advance of the fmPRC decreased with the natural
frequency ($CV^{(I)}=0.15$). (B) The fmPRC did not depend sensitively on the
within-network coupling strength $W$ ($CV^{(I)}=0.15$, $I^{(I)}=$15.8 pA). (C)
The Fourier spectrum of the LFP as a function of the coupling strength $W$.
Parameters as in (B). (D) Paradoxical phase response in the absence of an
explicit delay, $\tau_{d}=0$, for different effective synaptic delays due to
different synaptic rise times $\tau_{1}^{I}$ ($CV^{(I)}=0.05$, $I^{(I)}=$15.8
pA). For low $\tau_{1}^{I}$ (blue curve), the shift alternated in subsequent
cycles reflecting the subhamornic nature of the rhythm. (E) The Fourier
spectrum of the LFP as a function of the effective synaptic delay (synaptic
time constant of rise $\tau_{1}^{I}$). With decreasing $\tau_{1}^{I}$, a
subharmonic peak emerged and eventually the rhythm disintegrated. Parameters
as in D. In (A), (B) and (D), perturbations were made with a square-wave
inhibitory current pulse with duration 0.1 ms to each interneuron. In (A) and
(B), the amplitude of the current was 1600 pA, resulting in a 2 mV rapid
hyperpolarization. In (D), the amplitude of the current was 400 pA, resulting
in a 0.5 mV rapid hyperpolarization.

The paradoxical phase response was robust with respect to changes in the
natural frequency of the network, the coupling strength, and the effective
synaptic delay, as long as the rhythm persisted. The paradoxical phase advance
increased with decreasing natural frequency of the network, since the
inhibition had a stronger effect for lower mean input strength (Fig.3A).
Changing the within-network coupling strength by a factor of 2 up or down did
not substantially affect the paradoxical phase response (Fig.3B) nor the
strength of the rhythm (Fig.3C). Even without explicit synaptic delay
($\tau_{d}=0$), the effective delay given by the double-exponential synaptic
interaction was sufficient to render a paradoxical response (Fig.3D). However,
when this effective delay was reduced by decreasing the rise time
$\tau_{1}^{I}$ of the synaptic current, the rhythm itself developed a strong
subharmonic component and eventually disintegrated (Fig.3E). In the
subharmonic regime the paradoxical phase advance alternated in consecutive
cycles of the rhythm.

In [27, 13] the exact reduction of all-to-all coupled heterogeneous networks
of quadratic integrate-fire neurons to 2 coupled ordinary differential
equations for each network [23] has been used to obtain the infinitesimal
macroscopic phase-response curve (imPRC) for ING and PING networks. They
obtained biphasic response only if the excitatory perturbation was applied to
the population of inhibitory neurons; for perturbations to the excitatory
neurons they found only monophasic response (type-I). This is presumably due
to the lack of a delay in the single-exponential synaptic interactions used in
[27, 13].

### Enhancing entrainment of ING-rhythms through network heterogeneity

In order to allow communication by coherence [32, 11], the rhythms in
different brain areas need to be sufficiently phase-locked with each other. As
a simplification of two interacting $\gamma$-rhythms, we therefore
investigated the ability of the rhythm in a network to be entrained by a
periodic external input, particularly focusing on the possibly facilitating
role of neuronal heterogeneity. Motivated by the paradoxical phase response
induced by the heterogeneity, we addressed, in particular, the question
whether an ING network can be sped up by inhibition to entrain it with a
faster network.

Figure 4: Sketch of computational models. (A) ING rhythm receives periodic
inhibitory input generated from another ‘clock’ ING rhythm. (B) PING rhythm
receives periodic excitatory input by its E-population generated from another
‘clock’ PING rhythm.

The network considered here was the same as that used to analyze the fmPRC.
The within-network interaction was an all-to-all inhibition with synaptic
delay $\tau_{d}$, resulting in a rhythm with natural frequency $f_{natural}$,
Each neuron received heterogeneous input $I_{bias}$ and inhibitory periodic
pulses with frequency $f_{clock}$ . The latter can be considered as the output
of another ING-network and were, in fact, generated that way (Fig.4). We refer
to this external input as the ‘clock’. The detuning $\Delta
f=f_{clock}-f_{natural}$ was a key control parameter.

For periodic input the fmPRC allows the definition of an iterated map
describing the response of the network. For periodic $\delta$-pulses that map
is shown in Fig.5A. For positive detuning, i.e. when the clock is faster than
the network, entrainment requires that the phase response is paradoxical in
order for the rhythm to be sped up by the inhibition. If the heterogeneity and
the resulting phase response are sufficiently large, the maximum of the
iterated map crosses the diagonal, generating a stable and an unstable fixed
point. The former is the desired entrained state.

As the detuning is increased the iterated map is shifted downward. This can
decrease the slope of the iterated map at the fixed point below -1,
destabilizing the fixed point in a period-doubling bifurcation. For periodic
pulses comprised of double-exponential inhibitory currents (cf. (2,3)) a rich
bifurcation scenario emerged (Fig.5B). Note that the full map is not
continuous and not unimodal (cf. first bottom panel of Fig.5B). Nevertheless,
for $\Delta f<7.17$ Hz the attractor remains near the unstable fixed point and
displays a period-doubling cascade to chaos and multiple periodic windows. For
$\Delta f>7.28$ Hz, however, the attractor includes points on both sides of
the discontinuity (cf. third bottom panel in Fig.5B).

Figure 5: Connection between fmPRC and the synchronization of
$\gamma$-rhythms. (A) Iterated map of $\Phi_{inh}$. The network can be
synchronized by faster periodic inhibiton under sufficiently large advancing
phase response. (B) Top: The bifurcation diagram of the iterated map of
$\Phi_{inh}$ with varying detuning $\Delta f$. To the right of the magenta
dashed line ($\Delta f=7.28$ Hz) the attractors involve points on both sides
of the discountinuity of the map. Bottom from left to right: iterated maps of
$\Phi_{inh}$ for $\Delta f=0,2.44,8.8$ Hz. The distance between the diagonal
and subdiagonal line represents the detuning between the network and periodic
input. In (A), the fmPRC was determined for a $\delta$-pulse perturbation, in
(B) for a double-exponential inhibitory current (cf. (2,3)) was used as in
Fig.6.  Figure 6: Network heterogeneity enhances synchronization and phase-
locking of periodically driven ING rhythm. (A) Synchronization quantified
using $f_{dom}:f_{clock}$ with $f_{dom}$ and $f_{clock}$ being the dominant
frequencies of the Fourier spectrum of the LFP of the network and the clock,
respectively. The neuronal heterogeneity enhanced the synchronization by
shifting $f_{dom}$ to $f_{clock}$. Example 1: Synchronized with 1:1 phase-
locking. Example 2: Synchronized with subharmonic response (period 4). Example
3: synchronized with subharmonic response (chaotic). Example 4: Not
synchronized. Squares and dashed lines in the map of $\Phi_{inh}$ indicate
clock cycles in which the network did not spike ($\Phi_{inh}$ was arbitrarily
set to 2). (B) Subharmonic response. Color hue and saturation indicate the
frequency ratio $f_{sub}:f_{clock}$ and the ratio of the Fourier power at
these two frequencies. $f_{sub}$ is the frequency of the dominant peak of the
network power spectrum that satisfies $f_{sub}<f_{clock}$. The power ratio is
capped at 1\. Dashed line marks the value of input heterogeneity used in
Fig.5B. (C) Spiking variability and $var(\Phi_{inh})$ as a function of
neuronal heterogeneity and detuning. Color hue indicates the fraction of clock
cycles without spikes in the network. In particular, red indicates that the
network spikes in every cycle. Color saturation indicates $var(\Phi_{inh})$.
The neuronal heterogeneity enhances the tightness of the phase-locking.

Having clarified the role of the fmPRC in the network’s synchronizability and
ability to phase-lock, we investigated the role of neuronal heterogeneity in
more detail (Fig.6). To do that, we adjusted for each value of the input
heterogeneity the mean input strength $I^{(I)}$ so as to keep the natural
frequency $f_{network}$ constant ($f_{network}=44$ Hz). Then we determined the
extent of synchronization and phase-locking of the network under the influence
of periodic inhibitory input as a function of the detuning $\Delta f$ and
network heterogeneity $CV$. As shown above, the fmPRC of a heterogeneous
network could be biphasic with the amplitude of the paradoxical phase response
increasing with neuronal heterogeneity. Expecting that for sufficiently large
heterogeneity an ING-rhythm could be accelerated by a faster periodic
inhibition, we tested phase-locking predominantly for positive detuning,
corresponding to $f_{clock}>f_{network}$. We first investigated how neuronal
heterogeneity affected the synchronization by comparing the dominant frequency
$f_{dom}$ in the Fourier spectrum of the network’s LFP with $f_{clock}$. In
Fig.6A, the color hue indicates the ratio $f_{dom}:f_{clock}$. For small
heterogeneity, $f_{dom}$ was a rational multiple of $f_{clock}$ that depended
on the detuning, while for sufficiently large CV the network became
synchronized in the sense that $f_{dom}=f_{clock}$ (yellow). The range of
$\Delta f$ allowing synchronization became wider with increasing neuronal
heterogeneity, implying that the neuronal heterogeneity enhanced the
synchronization of the ING-rhythm. However, note that $f_{dom}=f_{clock}$ did
not imply a perfectly synchronized or a 1:1 phase-locked state. In fact,
various different subharmonic responses arose: example 2 shows a period-4
state, while in example 3 the dynamics were actually chaotic (Fig.6A) even
though $f_{dom}=f_{clock}$. Motivated by these observations, we divided the
states into three types:

  * •

Type 1: $f_{dom}\neq f_{clock}$, not synchronized, not phase-locked (Fig.6
example 4).

  * •

Type 2: $f_{dom}=f_{clock}$ with subharmonic response, might be poorly phase-
locked (Fig.6 example 3) or displaying rational ratio phase-locking (Fig.6
example 2).

  * •

Type 3: $f_{dom}=f_{clock}$, no subharmonic response, 1-to-1 phase-locking
(Fig.6 example 1).

The phase diagram Fig.6A does not differentiate between types 2 and 3. It only
shows that neuronal heterogeneity enhanced the synchronization of the network
by shifting $f_{dom}$ to $f_{clock}$. Therefore, we studied whether neuronal
heterogeneity also enhanced the synchronization by weakening the subharmonic
response and changing the synchronized state from type 2 to type 3, as well as
whether the dynamics of the fmPRC shown in the bifurcation diagram Fig.5B
could predict the phase relationship between the network and the clock. Using
the same simulation setup as in Fig.6A, the subharmonic response is shown in
Fig.6B. The color hue indicates the ratio $f_{sub}:f_{clock}$, where $f_{sub}$
is the frequency of the dominant peak of the LFP power spectrum that satisfies
$f_{sub}<f_{clock}$. The color saturation gives the ratio of the powers at
$f_{sub}$ and $f_{clock}$ (capped at 1). Thus, over most of the range of
positive detuning and neuronal heterogeneity tested, the fading-away of the
color with increasing heterogeneity reveals that the neuronal heterogeneity
weakened the subharmonic response. Over a small range of positive detuning,
increasing neuronal heterogeneity from small values induced perfect
synchronization (type 3) by weakening the subharmonic response with frequency
ratio $f_{sub}:f_{clock}=1:2$; the system traversed a continuous period-
doubling bifurcation in reverse, with type 2 (red) giving way to type 3
(white). Together with Fig.6A, this showed that neuronal heterogeneity could
enhance the synchronization both by making $f_{dom}=f_{clock}$ (from type 1 to
type 2) and by weakening the subharmonic response (from type 2 to type 3). The
range of detuning where increasing heterogeneity induced a type 3
synchronization became wider for larger synaptic delay within the network
(Suppl. Figure S3). Note that the bifurcation diagram (Fig.5B) based on the
fmPRC agrees well with the subharmonic response marked along the dashed line
at $CV=0.1$ in Fig.6B, suggesting that the fmPRC can well predict the
subharmonic response and persistent phase response of the network.

In addition to enhancing the frequency synchonization, neuronal heterogeneity
was also able to increase the tightness of the phase-locking. Over most of the
parameter regime investigated, the variance of the phase of the network
relative to the periodic input ($var(\Phi_{inh})$) decreased with increasing
heterogeneity, as indicated by the decrease in the color saturation in Fig.6C.
In fact, for detuning between 0 Hz and 2 Hz the heterogeneity reduced
$var(\Phi_{inh})$ to 0 (white), corresponding to the 1:1 phase-locked state.
Even for the 1:2 phase-locked state (cf. the red area in Fig.6B)
$var(\Phi_{inh})$ was very small for a range of heterogeneity and detuning (2
Hz to 4 Hz), indicating tight phase locking. Except for type-3 synchronized
states the size of the spike volleys varied between clock cycles. In fact,
over wide ranges of the parameters the network did not spike in each of the
clock cycles, as indicated by the color hue in Fig.6C, which gives the
fraction of cycles with no network spikes (e.g. Fig.6 example 4).

### Paradoxical phase response and entrainment of PING rhythms

Many $\gamma$-rhythms involve not only inhibitory neurons, but arise from the
mutual interaction of excitatory (E) and inhibitory (I) neurons (PING rhythm)
[33]. The key elements to obtain a paradoxical phase response and the ensuing
enhanced synchronization are self-inhibition within the network, neuronal
heterogeneity and effective synaptic delay. Since in PING rhythms the
connections from E-cells to I-cells and back to the E-cells form an effective
self-inhibiting loop, we asked whether PING-rhythms can exhibit behavior
similar to the behavior we identified for ING-rhythms.

Figure 7: Network heterogeneity enhances the synchronization and the tightness
of phase-locking of the PING rhythm. (A) fmPRC of PING networks with constant
natural frequency ($f_{network}=41$ Hz) but different neuronal heterogeneity.
Only with neuronal heterogeneity the phase was delayed by the excitation. (B)
Non-monotonicity of the paradoxical delay with constant natural frequency
($f_{network}=41$ Hz). B2-4: Top: raster plot of spikes in E-population (input
strength increased with cell index). Bottom: mean inhibitory synaptic
conductance within the PING network. The titles show the absolute and relative
increase in spike number (B2: $CV=0.05$, B3: $CV=0.1$, B4: $CV=0.2$). (C)
Subharmonic response of the PING rhythm with periodic excitation as function
of neuronal heterogeneity and detuning. $f_{network}$ was fixed at 41 Hz.
Color hue and saturation indicate the frequency ratio and power ratio at the
frequencies $f_{super}$ and $f_{clock}$ of the E-population’s LFP. $f_{super}$
was the frequency of the dominant peak of the LFP power spectrum that
satisfies $f_{super}>f_{clock}$. The power ratio was capped at 1. Generally,
the neuronal heterogeneity enhanced the synchronization of the PING rhythm by
weakening subharmonic response. (D) The tightness of the phase-locking
($var(\Phi_{exc})$) as a function of neuronal heterogeneity and detuning. The
neuronal heterogeneity enhanced the tightness of the phase-locking. For
$\Delta f\in[-22\mbox{Hz},\,-17.4\mbox{Hz}]$ the clock was twice as fast as
the network, resulting in vanishing $var(\Phi_{exc})$.

Considering a PING-rhythm generated by an E-I network comprised of integrate-
fire neurons, we first studied its fmPRC. To avoid that all I-cells receive
identical input and therefore spike as a single unit, the I-cells received, in
addition to the excitation from the E-cells, heterogeneous, tonic, Gaussian-
distributed subthreshold input with mean $I^{(I)}=36$ pA and $CV^{(I)}$=
0.167. The phase response of the network was probed by applying an identical
external excitatory perturbation to all E-cells and recording the resulting
phase shift of the LFP (cf. eqs.(7,8)) of the E-population, averaged across
500 realizations of the subthreshold input to the I-cells (Fig.7A). More
specifically, the perturbations consisted of a square-wave excitatory current
pulse with amplitude 76 pA and duration 0.1 ms to each E-cell, resulting in a
2 mV rapid depolarization. Without neuronal heterogeneity the external
excitation always advanced the phase of the rhythm resulting in an fmPRC that
was strictly positive. In the heterogeneous case, however, the PING rhythm
exhibited a paradoxical phase response, whereby the collective rhythm was
delayed while the individual neurons were advanced by the excitation. The
delay was caused by the increase of self-inhibition within the network that
was generated by the additional spikes in the E-population, which in turn
drove additional spikes in the I-population. In contrast to the fmPRC of the
ING-rhythm, this paradoxical phase response was not monotonic in the
heterogeneity. While weak heterogeneity resulted in strong delay, the delay
decreased with increasing intermediate CV-values and only increased again for
larger CV (Fig.7B left top). This non-monotonicity arose because we kept the
frequency of the network constant as we increased its heterogeneity. This
required a decrease in the tonic input to the E-cells with increasing
heterogeneity. For the stronger tonic input used for weak heterogeneity the
same external perturbation elicited more additional spikes than it did for
strong heterogeneity where the tonic input was weaker (cf. titles of subpanels
of Fig.7B). The total number of spikes occurring in each cycle of the
unperturbed network also decreased with increasing heterogeneity.
Consequently, the relative change in the number of spikes induced by the
perturbation was non-monotonic in the heterogeneity. As a result, the relative
change in the inhibitory synaptic conductance resulting from the perturbation
and with it the phase delay was also non-monotonic.

As for the ING rhythm, we investigated the role of neuronal heterogeneity in
the synchonizability and the ability of phase-locking of coupled PING rhythms.
In analogy to the ING-case, we considered the case of the E-population of a
PING network receiving periodic excitation generated by a clock PING network
(Fig.4B). As before, we adjusted the tonic input strength to the E-population
to keep the natural frequency of the network constant as we changed its
heterogeneity ($f_{network}=41$Hz). To probe the impact of the paradoxical
phase response on the synchronization we focused on negative detuning for
which the periodic external excitation needed to slow down the network in
order to achieve phase-locking. Indeed, with increasing heterogeneity the
network could become synchronized with the slower clock over a larger ranger
of the detuning as indicated by the fading saturation of the color in Fig.7C.
Here the color hue indicates the ratio $f_{super}:f_{clock}$, where
$f_{super}$ was determined as the frequency with the most power among the
frequencies higher than $f_{clock}$ in the Fourier spectrum of the
E-population’s LFP. The color saturation indicates the ratio of the power at
the frequencies $f_{super}$ and $f_{clock}$. Thus, a color hue closer to green
($f_{super}:f_{clock}=1:1$) or with a lower saturation implies better
synchronization. By observing how the width of the range of detuning allowing
synchronization varied with neuronal heterogeneity, we concluded that
generally, the synchronizability of PING rhythm was enhanced by the neuronal
heterogeneity by weakening subharmonic response. Note that for $CV\in[0,0.1]$
the synchonizability of the PING rhythm decreased slightly with neuronal
heterogeneity. This was consistent with the nonmonotonicity exhibited by the
fmPRC seen in Fig.7B. The neuronal heterogeneity played a similar role in the
tightness of the phase-locking as in the synchronizability (Fig.7D).

##  3 Discussion

In this paper we have analyzed the response of collective oscillations of
inhibitory and of excitatory-inhibitory networks of integrate-fire neurons to
external perturbations. For ING- and PING-rhythms we have shown that the
combination of neuronal heterogeneity and effective synaptic delay can
qualitatively change the phase response compared to the phase response of the
individual neurons generating the rhythm. Thus, perturbations that delay the
I-cells can paradoxically advance the ING-rhythm and perturbations that
advance the E-cells can delay the PING-rhythm. As a result, the macroscopic
phase-response curve for finite-amplitude perturbations (fmPRC) of the rhythm
can change sign as the phase of the perturbation is changed (type-II), even
though the PRC of all individual cells is strictly positive (type-I). This
change of the fmPRC enhances the ability of the $\gamma$-rhythm to synchronize
with other rhythms.

The key element of the mechanism driving the paradoxical phase response and
the enhanced synchronization is the cooperation of the external perturbation
and the effectively delayed within-network inhibition. In the ING-network a
suitably timed external perturbation delays the lagging - but not the early -
neurons sufficiently to allow the within-network inhibition triggered by the
early neurons to keep the lagging neurons from spiking. This reduces the
overall within-network inhibition and with it the duration of the cycle. Thus,
the perturbation modifies the internal dynamics of the rhythm, which leads to
changes in the phase of the rhythm that can dominate the immediate phase
change the perturbation induces. The situation is somewhat similar to that
investigated in [18]. There it had been pointed out that an external
perturbation of a collective oscillation can lead to changes in its phase in
two stages: i) an immediate change of the phases of all oscillators as a
direct result of the perturbation and ii) a subsequent slower change in the
collective phase resulting from the convergence of the disturbed phases back
to the synchronized state. That analysis was based on a network of phase
oscillators and could therefore not include a key element of our results,
which is the perturbation-induced change in the number of neurons that
actually spike and the resulting change in the within-network inhibition that
results in a change of the period of the rhythm. As discussed in [31, 34], for
ING-rhythms such a change in the number of spiking neurons underlies also the
enhanced phase-locking found in [30].

Going beyond ING-rhythms, we showed that PING-rhythms can also exibit a
paradoxical phase response _via_ a mechanism that is analogous to that of ING-
rhythms. For that analysis we have focused on excitatory-inhibitory networks
with only connections between but not within the excitatory and inhibitory
populations. For excitatory inputs to the excitatory cells to generate a
paradoxical phase response it is necessary that the additional spikes of the
excitatory neurons that are caused by the external perturbation induce
additional spikes of the inhibitory neurons. This behavior arises if the
inibitory population is also allowed to be heterogeneous. Moreover, the
within-network inhibition has to be strong enough to be able to suppress the
spiking of lagging excitatory neurons. This is, e.g., found in mice piriform
cortex, where principal neurons driven by sensory input from the olfactory
bulb arriving early during a sniff recruit inhibitory interneurons via long-
range recurrent connections, resulting in the global, transient suppression of
subsequent cortical activity [35]. A characteristic feature of the paradoxical
phase response of the PING rhythm is the extended cycle time following
enhanced activation of the excitatory cells. A strong such correlation between
the cycle time and the previous LFP amplitude has been observed for the
$\gamma$-rhythm in hippocampus [36]. To assess whether this rhythm exhibits
paradoxical phase response would require comparing the macroscopic phase
response [37] with that of indvidual participating neurons.

In order for the global perturbation to affect the various neurons
differently, they have to be at different phases in their cycle. Our analysis
suggests that the specific cause for this heterogeneity in the spike times
does not play an important role. Indeed, as shown in [29], even fluctuating
heterogeneities that are generated by noise rather than static heterogeneities
reflecting intrinsic properties of neurons can enhance the synchronization of
multiple $\gamma$-rhythms in interconnected networks of identical neurons.
Note that the noise driving this synchronization is uncorrelated across
neurons. The analysis of that system revealed the same mechanism at work as
the one identified here.

In various previous analytical and computational analyses it has been found
that the dynamics of the macroscopic phase of a collective oscillation can
qualitatively differ from that of the microscopic phase. Thus, for interacting
groups of noisy identical phase oscillators the macroscopic phases of the
groups can tend to lign up with each other, even if all pair-wise interactions
between individual oscillators prefer the antiphase state, and vice versa
[20]. An analogous result has been obtained for heterogeneous populations of
noiseless oscillators [17].

Qualitative changes have also been found in the macroscopic phase response of
rhythms in noisy homogeneous networks when the noise level was changed [21,
22, 28]. Using a Fokker-Planck approach for globally coupled excitable
neurons, a type-I mPRC was obtained for weak noise, where the rhythm emerges
through a SNIC bifurcation, while a type-II mPRC arose for strong noise that
led to a Hopf bifurcation [21]. A similar approach was used to obtain the mPRC
via the adjoint method for an extension of the theta-model that implements
conductance-based synaptic interactions. Again, although individual theta-
neurons have a type-I PRC, a type-II mPRC was obtained for the rhythm, which
arose in a Hopf bifurcation [22]. This was also the case in an extension to
networks of excitable and inhibitory neurons [28].

Thus, results reminiscent of those presented here have been obtained
previously. However, the mechanism underlying them was not addressed in detail
and remained poorly understood. We expect that our analysis will provide
insight into those systems. The key element of the mechanism discussed here is
that due to the dispersion of the spike times, which either results from
neuronal heterogeneity or noise, the external perturbation enables the within-
network inhibition to suppress the spiking of a larger number of neurons than
without it. In our system this was facilitated by the delay with which spikes
triggered the within-network inhibition, which allowed some neurons to escape
its impact in the absence of the external perturbation, but not in its
presence. Our analysis showed, however, that the explicit delay is not
necessary; the effective delay resulting from a double-exponential synaptic
interaction is sufficient. In fact, when reducing that effective delay the
paradoxical phase response did not disappear until the delay was so short that
the rhythm itself developed a strong subharmonic component and disintegrated.

In this paper we have focused on a specific, very simple neuronal model, the
leaky integrate-fire model with conductance-based pulsatile coupling. In
previous work on the enhanced synchronization among $\gamma$-rhythms _via_
noise-induced spiking heterogeneity it was demonstrated that this result does
not depend sensitively on the neuron type. Comparable results were obtained
also with Morris-Lecar neurons for parameters in which the periodic spiking
arises from a SNIC-bifurcation, resulting in a type-I PRC as is the case for
integrate-fire neurons, but also for parameters for which the spiking is due
to a Hopf bifurcation, resulting in a type-II PRC [29]. For networks comprised
of heterogeneous neurons with type-II PRC the fmPRC of the collective
oscillation is likely to be more complex, since the heterogeneity allows the
same input to induce phase shifts with opposite signs for different neurons.
However, we expect that the interplay between the within-network inhibition
and the external perturbation can again substantially and qualitatively modify
the fmPRC by changing the number of neurons participating in the rhythm.

In [29] the results were also found to be robust with respect to significant
changes in the network connectivity (random instead of all-to-all) as well as
the reversal potential of the inhibitory synapses, as long as the rhythm
itself persisted robustly (cf. [38]). In fact, the coupling did not even have
to be synaptic; collective oscillations of relaxation-type chemical
oscillators that were coupled diffusively were also shown to exhibit noise-
induced synchronization. These results suggest that the paradoxical phase
response found here arises in a much wider class of macroscopic collective
oscillations.

The strong paradoxical phase response that we demonstrated for heterogeneous
networks allows their rhythm to synchronize with a periodic external input
over a range of detuning that increases substantially with the neuronal
heterogeneity. This is reminiscent of computational results for anterior
cingulate cortex that investigated networks of excitatory neurons coupled via
a common population of inhibitory neurons. There heterogeneity was also found
to enhance the synchrony of rhythms, as measured in terms of coincident spikes
within 10ms bins [39].

The heterogeneity-enhanced synchrony we have identified suggests that the
coherence of $\gamma$-rhythms emerging in different interacting networks could
also be enhanced by neuronal heterogeneity. It has been proposed that the
coherence of different $\gamma$-rhythms, which has been observed to be
modified by attention [8], plays an important role in the communication
between the corresponding networks [32, 11]. Computational studies have shown
that the direction of information transfer between networks depends on the
relative phase of their rhythms [13, 12], which can be changed by switching
between different base states [40, 41]. Whether the enhanced synchrony
resulting from neuronal heterogeneity enhances this information transfer is
still an open question.

Disrupted $\gamma$-rhythms have been observed in multiple brain regions in
neurological diseases, especially Alzheimer’s disease. Optogenetic and sensory
periodic stimulation at $\gamma$-frequencies has been found to entrain the
$\gamma$-rhythm in hippocampus and visual cortex, respectively, and has
resulted in a significant reduction in total amyloid level [42]. Similar
neuro-protective effects of entrainment by external $\gamma$-stimulation have
also been found for other sensory modalities [43, 14]. This suggests that
$\gamma$-stimulation by sensory input might be a feasible therapeutic
approach. Our results suggest a potential role of neuronal heterogeneity in
this context.

From a functional perspective, it has been shown that the noise-induced
synchronization mentioned above can facilitate certain learning processes
[44]. Specifically, a read-out neuron was considered that received input from
neurons in two networks _via_ synapses that exhibited spike-timing dependent
plasticity. The two networks were interacting with each other and each of them
exhibited a $\gamma$-rhythm, albeit at different frequencies. For low noise
the two rhythms were not synchronized and the read-out neuron received inputs
from the two networks at uncorrelated times. These inputs drove the plasticity
inconsistently, leading only to a very slow overall evolution of the synaptic
weights, if any. However, for stronger noise the two networks were
synchronized, providing a more consistent spike timing that lead to
substantial changes in the synaptic weights. As a result, the read-out neuron
was eventually only driven by the network that had the larger natural
frequency in the absence of the coupling between the networks. It is expected
that synchrony by neuronal heterogeneity will have a similar impact.

##  4 Methods

Neuron model. Both E-cells and I-cells were modeled as leaky integrate-and-
fire neurons, each characterized by a membrane potential $V_{i}(t)$ satisfying

$\tau_{E,I}\frac{d}{dt}V_{i}=-(V_{i}-V_{rest})+\frac{I_{i}^{(syn)}}{g_{syn}}+\frac{I_{i}^{(ext)}}{g_{ext}}+\frac{I_{i}^{(bias)}}{g_{bias}}\,,$
(1)

where $V_{rest}$ is the resting potential and $\tau_{E,I}$ the membrane time
constants of the E- and I-cells, respectively. $I_{i}^{(syn)}(t)$ is the total
synaptic current that the neuron receives from the other neurons within the
network. $I_{i}^{(ext)}(t)$ is a time-dependent external input that represents
perturbations applied to determine the fmPRC or, in the study of
synchronization, the periodic input generated by the clock network.
$I_{i}^{(bias)}$ denotes a tonic, neuron-specific excitatory bias current that
implements the heterogeneity of the neuron properties. The corresponding
conductances are denoted $g_{syn}$, $,g_{ext}$, and $g_{bias}$. Upon the
$i^{\mbox{th}}$ neuron reaching the spiking threshold $V_{peak}$, the voltage
$V_{i}$ was reset to the fixed value $V_{reset}$. Parameters for the neuron
were kept fixed throughout all simulations (see Table 1). The local field
potential (LFP) of the network was approximated as the mean voltage across all
neurons $j=1,...N$ in the respective population.

Network model. We studied two types of networks: an ING network and a PING
network. The ING network was modeled as an all-to-all inhibitory network of
$N_{I}^{(ING)}$ interneurons. The PING network was modeled as a network of
$N_{I}^{(PING)}$ interneurons and $N_{E}^{(PING)}$ principal cells with all-
to-all interneuron-principal and principal-interneuron connections (i.e.,
without principal-principal and interneuron-interneuron connections). In PING,
only principal cells received external input $I^{ext}(t)$.

To gain insight into the interaction between two ING rhythms, we considered
the simplified situation in which all neurons in the network received strictly
periodic input $I^{(ext)}$, which was generated by another ING network
(‘clock’). Similarly, for PING rhythms, the E-cells of the PING network
received strictly periodic excitatory input $I^{(ext)}$ from another PING
network through all-to-all connection between their E populations.

Synaptic currents. We used delayed double-exponential conductance-based
currents to model the excitatory and the inhibitory synaptic inputs from
neuron $j$ to neuron $i$,

$I_{ij}^{(syn)}(t)=g_{syn}\frac{\tau_{E,I}}{\tau_{2}^{E,I}-\tau_{1}^{E,I}}\left(A_{ij}^{(2)}(t)-A_{ij}^{(1)}(t)\right)\,\left(V_{rev,j}-V_{i}(t)\right)\,,$
(2)

with the two exponentials $A_{ij}^{(1,2)}(t)$ satisfying

$\frac{d}{dt}A_{ij}^{(1,2)}(t)=-\frac{A_{ij}^{(1,2)}(t)}{\tau_{1,2}^{E,I}}+\sum_{k}W_{ij}\delta(t-t_{j}^{(k)}-\tau_{d})\,,$
(3)

where $V_{rev,j}$ is the synaptic reversal potential corresponding to the
synapse type, $W_{ij}$ the dimensionless synaptic strength, and $\delta$ the
Dirac $\delta$-function. All synapses of the same type (I-I, I-E, E-I) were
equally strong. The time constants of $A_{i}^{(1,2)}(t)$ satisfied
$\tau_{2}^{E,I}>\tau_{1}^{E,I}$. The synaptic current was normalized to render
the time integral independent of the synaptic time constants
$\tau_{1,2}^{E,I}$. The inhibitory synaptic currents had a slower decay than
the excitatory ones (cf. Table 1). We included an explicit synaptic delay
$\tau_{d}$ in the model. Every spike of the presynaptic neuron $j$ at time
$t_{j}^{(k)}$ triggered a jump in both $A_{ij}^{(1,2)}(t)$, making the
synaptic conductance rise continously after a synaptic delay $\tau_{d}$.

External periodic inputs were also modeled as double-exponential conductance-
based currents with $g_{syn}$ in (2,3) replaced by $g_{ext}$. The time
constants and delay were as for the within-network synaptic inputs.

Heterogeneous tonic input. The bias currents $I_{i}^{(bias)}$ of the ING
network were Gaussian distributed around $I_{mean}$ with a coefficient of
variation $CV$ and arranged in increasing order,
$I_{1}^{(bias)}>I_{2}^{(bias)}...>I_{N}^{(bias)}$. For the PING network, all
excitatory neurons received a heterogeneous bias $I_{E}^{(bias)}$ with mean
$I^{(E)}$ and a coefficient of variation $CV^{(E)}$. Similarly, the bias
currents $I_{I}^{(bias)}$ to the inhibitory neurons were characterized by
their mean $I^{(I)}$and their coefficient of variation $CV^{(I)}$. Without the
excitatory input from principal cells, the voltage of interneurons remained
below the spiking threshold. In our investigation of the impact of the
neuronal heterogeneity on the phase response and entrainment of the PING
rhythm we kept $CV^{(I)}$ fixed and varied $CV^{(E)}$.

Macroscopic Phase-response Curve for Finite-Amplitude Perturbations (fmPRC).

ING rhythm. For a single ING network, we applied a single inhibitory
$\delta$-pulse to each neuron $j=1,...N_{I}^{(ING)}$ at time $t_{inh}$ (dashed
green line in Fig.1B) and recorded the resulting phase shift $\Delta\varphi$.
The amplitude of the inhibitory perturbation to each neuron was the same. The
phase of the inhibition was defined as

$\phi_{inh}=\frac{t_{inh}-t_{firstspike}^{(unperturbed)}}{T},$ (4)

where $T$ was the period of the network LFP and
$t_{firstspike}^{(unperturbed)}$ the time of the first spike in the spike
volley of the unperturbed network that was closest to $t_{inh}$. The resulting
phase shift $\Delta\phi$ was given by

$\Delta\phi=\frac{\left(t_{firstspike}^{(unperturbed)}-t_{firstspike}^{(perturbed)}\right)}{T},$
(5)

where $t_{firstspike}^{(perturbed)}$ is the time of the first spike in the
corresponding volley in the perturbed network. $\Delta\phi$ and $\phi_{inh}$
were taken to be in the range $[-0.5\,\>0.5)$. Positive $\Delta\phi$ indicated
that the network was advanced by the perturbation, while negative indicated a
delay.

The periodic input (‘clock’) that was used to test the synchronizability of
the ING-rhythm was generated by a homogeneous ING network. The phase of the
periodic input in the $n^{\mbox{th}}$ clock cycle was defined by

$\Phi_{inh}^{(n)}=\frac{\left(t_{firstspike}^{(clock)(n)}+\tau_{d}-t_{firstspike}^{(network)(n)}\right)}{T},$
(6)

where $t_{firstspike}^{(network)(n)}$ was the time of the first spike in the
spike volley of the network in the $n^{\mbox{th}}$ cycle and
$t_{firstspike}^{(clock)(n)}$ the time of the spike of the clock. In contrast
to the definition of $\phi_{inh}$ in (4), the definition of $\Phi_{inh}^{(n)}$
included the delay $\tau_{d}$, since the inhibition generated by the clock
arrived with delay $\tau_{d}$ in the network.

PING rhythm. To probe the phase response of the PING network we used the same
approach as for the ING rhythm, except that we used excitatory instead of
inhibitory $\delta$-pulses and applied them only to the E-cells. The phase of
the excitation $\phi_{exc}$ and the resulting phase shift $\Delta\phi$ were
determined similarly as in the case of the ING rhythm,

$\phi_{exc}=\frac{t_{exc}-t_{firstspike}^{(unperturbed)}}{T},$ (7)

$\Delta\phi=\frac{(t_{firstspike}^{(unperturbed)}-t_{firstspike}^{(perturbed)})}{T},$
(8)

where $t_{firstspike}^{(perturbed)}$ and $t_{firstspike}^{(unperturbed)}$ were
the times of the first spike in the respective spike volleys of the
E-population.

Analogously to $\Phi_{inh}^{(n)}$, the phase of the periodic input in the
$n^{\mbox{th}}$ clock cycle was given by

$\Phi_{exc}^{(n)}=\frac{(t_{firstspike}^{(clock)(n)}+\tau_{d}-t_{firstspike}^{(network)(n)})}{T}.$
(9)

Throughout, the tonic, Gaussian distributed input to the interneurons in the
PING network was fixed: $I^{(I)}=36$ pA, $CV^{(I)}=0.167$.

ING network |
---|---
Parameter |  Value
$\tau_{I}$, membrane time constant |  20 ms
$u_{rest}$, resting potential |  -55 mV
$V_{peak}$, spiking threshold |  -50 mV
$V_{reset}$, reset voltage |  -60 mV
$\tau_{d}$, synaptic delay |  3 ms
$N_{I}^{(ING)}$, # of interneurons |  500
$W$, synaptic strength within the network |  $7.5\times 10^{-3}$
$W^{(ext)}$, synaptic strength for the input from the clock network |  $1.8\times 10^{-3}$

Synaptic currents Parameter Value $\tau_{1}^{E}$, time constant of rise in
excitatory synapse 0.5 ms $\tau_{2}^{E}$, time constant of decay in excitatory
synapse 2 ms $\tau_{1}^{I}$, time constant of rise in inhibitory synapse 0.5
ms $\tau_{2}^{I}$, time constant of decay in inhibitory synapse 5 ms
$V_{rev}^{I}$, reversal potential of inhibitory synapse -70 mV $V_{rev}^{E}$,
reversal potential of excitatory synapse 0 mV

PING network |
---|---
Parameter |  Value
$\tau_{E}$, membrane time constant of principal cells |  20 ms
$\tau_{I}$, membrane time constant of interneurons |  10 ms
$u_{rest}$, resting potential |  -70 mV
$V_{peak}$, spiking threshold |  -52 mV
$V_{reset}$, reset voltage |  -59 mV
$\tau_{d}$, synaptic delay |  1 ms
$N_{I}^{(PING)}$, # of interneurons |  200
$N_{E}^{(PING)}$, # of principal cells |  800
$W^{I}$, inhibitory synaptic strength within the network |  $5.4\times 10^{-3}$
$W^{E}$, excitatory synaptic strength within the network |  $1.67\times 10^{-3}$
$W^{(ext)}$, clock-network synaptic strength |  $1.6\times 10^{-3}$

Synaptic conductances Parameter Value Excitatory input on principal cells :
$g_{ext}^{(PING)}$, $g_{bias\>E}^{(PING)}$ 0.19 nS Excitatory input on
interneurons: $g_{bias}^{(ING)}$, $g_{syn\>EtoI}^{(PING)}$,
$g_{bias\>I}^{(PING)}$ 0.3 nS Inhibitory input on principal cells:
$g_{syn\>ItoE}^{(PING)}$ 2.5 nS Inhibitory input on interneurons:
$g_{ext}^{(ING)}$, $g_{syn}^{(ING)}$ 4 nS

Table 1: Parameters used in the model. Most parameters are based on [29, 45].

## References

  * [1] Bruesselbach H, Jones DC, Mangir MS, Minden M, Rogers JL.  Self-organized coherence in fiber laser arrays.  Optics Letters. 2005;30(11):1339–1341.  doi:10.1364/OL.30.001339.
  * [2] Wiesenfeld K, Colet P, Strogatz SH.  Synchronization Transitions in a Disordered Josephson Series Array.  Phys Rev Lett. 1996;76(3):404–407.
  * [3] Liu C, Weaver DR, Strogatz SH, Reppert SM.  Cellular Construction of a Circadian Clock: Period Determination in the Suprachiasmatic Nuclei.  Cell. 1997;91(6):855–860.
  * [4] Venzin OF, Oates AC.  What are you synching about? Emerging complexity of Notch signaling in the segmentation clock.  Developmental Biology. 2020;460(1):40–54.  doi:10.1016/j.ydbio.2019.06.024.
  * [5] Wang XJ.  Neurophysiological and Computational Principles of Cortical Rhythms in Cognition.  Physiol Rev. 2010;90(3):1195–1268.
  * [6] Börgers C, Epstein S, Kopell NJ.  Background gamma rhythmicity and attention in cortical local circuits: a computational study.  Proc Natl Acad Sci U S A. 2005;102(19):7002–7007.
  * [7] Börgers C, Kopell NJ.  Gamma oscillations and stimulus selection.  Neural Comput. 2008;20(2):383–414.
  * [8] Bosman CA, Schoffelen JM, Brunet N, Oostenveld R, Bastos AM, Womelsdorf T, et al.  Attentional stimulus selection through selective synchronization between monkey visual areas.  Neuron. 2012;75:875–888.  doi:10.1016/j.neuron.2012.06.037.
  * [9] Roberts MJ, Lowet E, Brunet NM, Ter Wal M, Tiesinga P, Fries P, et al.  Robust Gamma Coherence between Macaque V1 and V2 by Dynamic Frequency Matching.  Neuron. 2013;78(3):523–536.  doi:10.1016/j.neuron.2013.03.003.
  * [10] Buzsaki G, Schomburg EW.  What does gamma coherence tell us about inter-regional neural communication?  Nature Neuroscience. 2015;18(4):484–489.  doi:10.1038/nn.3952.
  * [11] Fries P.  Rhythms for Cognition: Communication through Coherence.  Neuron. 2015;88:220–235.  doi:10.1016/j.neuron.2015.09.034.
  * [12] Palmigiano A, Geisel T, Wolf F, Battaglia D.  Flexible information routing by transient synchrony.  Nature Neuroscience. 2017;20:1014–1022.  doi:10.1038/nn.4569.
  * [13] Dumont G, Gutkin B.  Macroscopic phase resetting-curves determine oscillatory coherence and signal transfer in inter-coupled neural circuits.  PLoS Computational Biology. 2019;15:e1007019.  doi:10.1371/journal.pcbi.1007019.
  * [14] Adaikkan C, Middleton SJ, Marco A, Pao PC, Mathys H, Kim DNW, et al.  Gamma Entrainment Binds Higher-Order Brain Regions and Offers Neuroprotection.  Neuron. 2019;102:929–943.e8.  doi:10.1016/j.neuron.2019.04.011.
  * [15] Schultheiss NW, Butera RJ.  Phase Response Curves in Neuroscience.  Springer; 2012.
  * [16] Brown E, Moehlis J, Holmes P.  On the Phase Reduction and Response Dynamics of Neural Oscillator Populations.  Neural Comput. 2004;16(4):673–715.
  * [17] Kawamura Y, Nakao H, Arai K, Kori H, Kuramoto Y.  Phase synchronization between collective rhythms of globally coupled oscillator groups: noiseless nonidentical case.  Chaos (Woodbury, NY). 2010;20:043110.  doi:10.1063/1.3491346.
  * [18] Levnajic Z, Pikovsky A.  Phase resetting of collective rhythm in ensembles of oscillators.  Physical review E, Statistical, nonlinear, and soft matter physics. 2010;82:056202.  doi:10.1103/PhysRevE.82.056202.
  * [19] Kawamura Y, Nakao H, Arai K, Kori H, Kuramoto Y.  Collective Phase Sensitivity.  Phys Rev Lett. 2008;101(2):024101.
  * [20] Kawamura Y, Nakao H, Arai K, Kori H, Kuramoto Y.  Phase Synchronization Between Collective Rhythms of Globally Coupled Oscillator Groups: Noisy Identical Case.  Chaos. 2010;20(4):043109.
  * [21] Kawamura Y, Nakao H, Kuramoto Y.  Collective phase description of globally coupled excitable elements.  Physical review E, Statistical, nonlinear, and soft matter physics. 2011;84:046211.  doi:10.1103/PhysRevE.84.046211.
  * [22] Kotani K, Yamaguchi I, Yoshida L, Jimbo Y, Ermentrout GB.  Population dynamics of the modified theta model: macroscopic phase reduction and bifurcation analysis link microscopic neuronal interactions to macroscopic gamma oscillation.  J R Soc Interface. 2014;11(95):20140058.  doi:10.1098/rsif.2014.0058.
  * [23] Montbrio E, Pazo D, Roxin A.  Macroscopic description for networks of spiking neurons.  Physical Review X. 2015;5:021028.
  * [24] Ott E, Antonsen TM.  Low dimensional behavior of large systems of globally coupled oscillators.  Chaos. 2008;18(3):037113.  doi:10.1063/1.2930766.
  * [25] Luke TB, Barreto E, So P.  Complete Classification of the Macroscopic Behavior of a Heterogeneous Network of Theta Neurons.  Neural Computation. 2013;25(12):3207–3234.  doi:10.1162/NECO_a_00525.
  * [26] Hannay KM, Booth V, Forger DB.  Collective phase response curves for heterogeneous coupled oscillators.  Physical Review E. 2015;92(2):022923.  doi:10.1103/PhysRevE.92.022923.
  * [27] Dumont G, Ermentrout GB, Gutkin B.  Macroscopic phase-resetting curves for spiking neural networks.  Physical Review E. 2017;96:042311.  doi:10.1103/PhysRevE.96.042311.
  * [28] Akao A, Ogawa Y, Jimbo Y, Ermentrout GB, Kotani K.  Relationship between the mechanisms of gamma rhythm generation and the magnitude of the macroscopic phase response function in a population of excitatory and inhibitory modified quadratic integrate-and-fire neurons.  Physical Review E. 2018;97(1):012209.  doi:10.1103/PhysRevE.97.012209.
  * [29] Meng JH, Riecke H.  Synchronization by uncorrelated noise: interacting rhythms in interconnected oscillator networks.  Scientific Reports. 2018;8:6949.  doi:10.1038/s41598-018-24670-y.
  * [30] White JA, Banks MI, Pearce RA, Kopell NJ.  Networks of interneurons with fast and slow gamma-aminobutyric acid type A (GABA(A)) kinetics provide substrate for mixed gamma-theta rhythm.  Proc Nat Acad Sci USA. 2000;97(14):8128–8133.  doi:10.1073/pnas.100124097.
  * [31] Serenevy AK, Kopell NJ.  Effects of Heterogeneous Periodic Forcing on Inhibitory Networks.  Siam Journal on Applied Dynamical Systems. 2013;12(3):1649–1684.  doi:10.1137/12089274X.
  * [32] Fries P.  A Mechanism for Cognitive Dynamics: Neuronal Communication Through Neuronal Coherence.  Trends Cogn Sci. 2005;9(10):474–480.
  * [33] Bartos M, Vida I, Jonas P.  Synaptic mechanisms of synchronized gamma oscillations in inhibitory interneuron networks.  Nature Reviews Neuroscience. 2007;8(1):45–56.  doi:10.1038/nrn2044.
  * [34] Cannon J, McCarthy MM, Lee S, Lee J, Börgers C, Whittington MA, et al.  Neurosystems: brain rhythms and cognitive processing.  European J Neuroscience. 2014;39:705–719.  doi:10.1111/ejn.12453.
  * [35] Bolding KA, Franks KM.  Recurrent cortical circuits implement concentration-invariant odor coding.  Science (New York, NY). 2018;361.  doi:10.1126/science.aat6904.
  * [36] Atallah BV, Scanziani M.  Instantaneous Modulation of Gamma Oscillation Frequency by Balancing Excitation with Inhibition.  Neuron. 2009;62(4):566–577.  doi:10.1016/j.neuron.2009.04.027.
  * [37] Akam T, Oren I, Mantoan L, Ferenczi E, Kullmann DM.  Oscillatory dynamics in the hippocampus support dentate gyrus-CA3 coupling.  Nature Neuroscience. 2012;15(5):763–768.  doi:10.1038/nn.3081.
  * [38] Börgers C, Krupa M, Gielen S.  The Response of a Classical Hodgkin-Huxley Neuron To an Inhibitory Input Pulse.  J Comput Neurosci. 2010;28(3):509–526.
  * [39] Adams NE, Sherfey JS, Kopell NJ, Whittington MA, LeBeau FEN.  Hetereogeneity in Neuronal Intrinsic Properties: A Possible Mechanism for Hub-Like Properties of the Rat Anterior Cingulate Cortex during Network Activity.  Eneuro. 2017;4(1):UNSP e0313–16.2017.  doi:10.1523/ENEURO.0313-16.2017.
  * [40] Witt A, Palmigiano A, Neef A, El Hady A, Wolf F, Battaglia D.  Controlling the oscillation phase through precisely timed closed-loop optogenetic stimulation: a computational study.  Front Neural Circuits. 2013;7:49.  doi:10.3389/fncir.2013.00049.
  * [41] Kirst C, Timme M, Battaglia D.  Dynamic information routing in complex networks.  Nature Communications. 2016;7:11061.  doi:10.1038/ncomms11061.
  * [42] Iaccarino HF, Singer AC, Martorell AJ, Rudenko A, Gao F, Gillingham TZ, et al.  Gamma frequency entrainment attenuates amyloid load and modifies microglia.  Nature. 2016;540(7632):230–+.  doi:10.1038/nature20587.
  * [43] Martorell AJ, Paulson AL, Suk HJ, Abdurrob F, Drummond GT, Guan W, et al.  Multi-sensory Gamma Stimulation Ameliorates Alzheimer’s-Associated Pathology and Improves Cognition.  Cell. 2019;177:256–271.e22.  doi:10.1016/j.cell.2019.02.014.
  * [44] Meng JH, Riecke H.  Synchronization by uncorrelated noise: interacting rhythms in interconnected neuronal networks.  BMC Neuroscience. 2018;19 (Suppl 2):116.
  * [45] Brunel N, Wang XJ.  What Determines the Frequency of Fast Network Oscillations With Irregular Neural Discharges? I. Synaptic Dynamics and Excitation-Inhibition Balance.  J Neurophysiol. 2003;90(1):415–430.

## Supplementary Information

Figure S1: Dependence of the phase dispersion on the heterogeneity of the bias
current. The phase dispersion was determined as the time difference between
the first and the last spike in the same spike volley normalized by the
period. Blue: fixed natural frequency ($f_{network}=40$Hz) for different
neuronal heterogeneity. Red: fixed mean input strength ($I^{(I)}=$15.8 pA) for
different neuronal heterogeneity. For $CV\geq$ 0.075 (dashed line), some
neurons spike more than once in a cycle.

Figure S2: fmPRC of heterogeneous ING networks for fixed steady current
($I^{(I)}=$15.8 pA) instead of fixed frequency (cf. Fig.2A). The paradoxical
phase advance increased with neuronal heterogeneity.

Figure S3: Subharmonic response of the ING rhythm with a longer synaptic delay
within the network ($\tau_{d}=5$ ms) receiving periodic inhibitory input. For
each value of the input heterogeneity, the natural frequency $f_{network}$ was
kept constant ($f_{network}=44$ Hz) by adjusting the mean input strength
$I^{(I)}$. The range of detuning where increasing heterogeneity induced a type
3 synchronization became wider compared to Fig.6B, where $\tau_{d}=3$ ms.
$W^{(ext)}=1.2\times 10^{-3}$.


* Whitfield et al. (2016) Margot Whitfield, Raza Mir Ali, and Deborah Fels. 2016.  LiveDescribe web redefining what and how entertainment content can be accessible to blind and low vision audiences. In _Computers Helping People with Special Needs: 15th International Conference, ICCHP 2016, Linz, Austria, July 13-15, 2016, Proceedings, Part I 15_. Springer, 224–230.  https://doi.org/10.1177/0145482X1210600304
  * Whittaker et al. (2019) Meredith Whittaker, Meryl Alper, Cynthia L Bennett, Sara Hendren, Liz Kaziunas, Mara Mills, Meredith Ringel Morris, Joy Rankin, Emily Rogers, Marcel Salas, et al. 2019\.  Disability, bias, and AI.  _AI Now Institute_ 8 (2019).  https://ainowinstitute.org/wp-content/uploads/2023/04/disabilitybiasai-2019.pdf
  * Wu et al. (2017) Shaomei Wu, Jeffrey Wieland, Omid Farivar, and Julie Schiller. 2017.  Automatic alt-text: Computer-generated image descriptions for blind users on a social network service. In _Proceedings of the 2017 ACM Conference on Computer Supported Cooperative Work and Social Computing_. 1180–1192.  https://doi.org/10.1145/2998181.2998364
  * Yao et al. (2015) Li Yao, Atousa Torabi, Kyunghyun Cho, Nicolas Ballas, Christopher Pal, Hugo Larochelle, and Aaron Courville. 2015.  Describing videos by exploiting temporal structure. In _Proceedings of the IEEE international conference on computer vision_. 4507–4515.  https://doi.org/10.1109/ICCV.2015.512
  * YouDescribe (2023) YouDescribe. 2023\.  _YouDescribe - Audio Description for YouTube Videos_.  https://youdescribe.org/
  * YouTube (2016) YouTube. 2016.  _The latest YouTube stats on when, where, and what people watch_.  https://www.thinkwithgoogle.com/data-collections/youtube-stats-video-consumption-trends/
  * Yuksel et al. (2020) Beste F Yuksel, Pooyan Fazli, Umang Mathur, Vaishali Bisht, Soo Jung Kim, Joshua Junhee Lee, Seung Jung Jin, Yue-Ting Siu, Joshua A Miele, and Ilmi Yoon. 2020\.  Human-in-the-Loop Machine Learning to Increase Video Accessibility for Visually Impaired and Blind Users. In _Proceedings of the 2020 ACM Designing Interactive Systems Conference_. 47–60.  http://doi.org/10.1145/3357236.3395433
  * Zhang et al. (2023b) Chaoyi Zhang, Kevin Lin, Zhengyuan Yang, Jianfeng Wang, Linjie Li, Chung-Ching Lin, Zicheng Liu, and Lijuan Wang. 2023b.  MM-Narrator: Narrating Long-form Videos with Multimodal In-Context Learning.  _arXiv preprint arXiv:2311.17435_ (2023).  https://arxiv.org/pdf/2311.17435.pdf
  * Zhang et al. (2022) Lotus Zhang, Jingyao Shao, Augustina Ao Liu, Lucy Jiang, Abigale Stangl, Adam Fourney, Meredith Ringel Morris, and Leah Findlater. 2022.  Exploring Interactive Sound Design for Auditory Websites. In _Proceedings of the 2022 CHI Conference on Human Factors in Computing Systems_. 1–16.  https://doi.org/10.1145/3491102.3517695
  * Zhang et al. (2023a) Zhuohao Jerry Zhang, Smirity Kaushik, JooYoung Seo, Haolin Yuan, Sauvik Das, Leah Findlater, Danna Gurari, Abigale Stangl, and Yang Wang. 2023a.  $\\{$ImageAlly$\\}$: A $\\{$Human-AI$\\}$ Hybrid Approach to Support Blind People in Detecting and Redacting Private Image Content. In _Nineteenth Symposium on Usable Privacy and Security (SOUPS 2023)_. 417–436.  https://www.usenix.org/conference/soups2023/presentation/zhang
  * Zhou et al. (2018) Luowei Zhou, Chenliang Xu, and Jason Corso. 2018.  Towards automatic learning of procedures from web instructional videos. In _Proceedings of the AAAI Conference on Artificial Intelligence_ , Vol. 32.  https://doi.org/10.1609/aaai.v32i1.12342

##  Appendix A All Participant-Specific Scenarios

Table 6. All 17 of our participant-specific video probes. Titles are
hyperlinked to the video on YouTube. Some videos are linked at specific
timestamps to capture the video start time used during our studies.

17 unique user-grounded video probes. Columns from left to right: Pseudonym,
Type, Platform, Goal, Description, and Title / Link.  Pseudonym Type Platform
Goal Description Title / Link Alice Vlog TVS Engage with friends /
entertainment A day in the life of a software engineer day in my life as a
software engineer in NYC * in-office edition * Blair Sports TVS Entertainment
A highlight reel of a professional women’s tennis match Beatriz Haddad Maia
vs. Leylah Fernandez — 2023 Montreal Round 2 — WTA Match Highlights Colin
Music video TVS Entertainment Music video for Take On Me by a-ha a-ha - Take
On Me (Official Video) [Remastered in 4K] Diana Video game / how-to TVS Learn
how to do something A playthrough of the video game, The Last of Us What Have
I Gotten Myself Into… * The Last of Us First Playthrough * Part 1 Comedic SNS
Entertainment A short video of a dog Dogs funny reaction to entering optical
illusion rug! #shorts Emily How-to / DIY TVS Learn how to do something
Tutorial for DIY desk upgrades DIY Desk Upgrades Felix Action / foreign
language SS Entertainment A fight scene from a Korean drama, Vincenzo Vincenzo
Cassano – Tailor Fight Scene Grace Cooking / how-to TVS Learn how to do
something A video explaining how to cook four different meals 4 Meals Anyone
Can Make Informational TVS Learn about a person, event, or idea A video about
Manhattan’s grid plan Where Manhattan’s grid plan came from Haley Music video
TVS Entertainment A montage of animated characters engaging in adventure Theme
Song — Elena of Avalor — @disneyjunior Isaac Music video TVS Entertainment
Music video for Bad Blood by Taylor Swift ft. Kendrick Lamar Taylor Swift -
Bad Blood ft. Kendrick Lamar Julia Live theater TVS Entertainment A musical
number from the Broadway show, Moulin Rouge Moulin Rouge! The Musical on Good
Morning America Karla Exercise / how-to TVS Learn how to do something A
fitness instructor goes through a warm-up routine Pumped Up Cardio Warmup!
(Easy, fun, at home workout) Layne Video game / how-to TVS Learn how to do
something A playthrough of two minigames in Mario Party Mario Party Superstars
ALL MINIGAMES!! Mason Cooking / how-to TVS Learn how to do something An
instructional video of five minute meals 7 Recipes You Can Make In 5 Minutes
Nicki Music video TVS Entertainment Music video for You Belong With Me by
Taylor Swift Taylor Swift - You Belong With Me Oscar Action SS Entertainment A
fight scene from The Avengers Thor vs Hulk - Fight Scene - The Avengers (2012)
Movie Clip HD

##  Appendix B Summary of Findings

Table 7. Summary of desired details and output modalities for nine different
scenarios, represented here by their video types.
Fantastical = Science Fiction, Fantasy, and Animation
Drama = Historical, Romance, Reality, and Drama

Summary of desired details and output modalities for nine different scenarios,
represented here by their video types. Desired details include actions, text
on screen, subjects, settings, clothing, visual aids, scene changes, facial
expressions, visual effects, and equipment. Output modalities include ambient
sound, additional resources, audio cues, tactile graphics, text to Braille, 3D
models, background music, and vibration.  How-To Info / Edu Short-Form Music
Live Personal Fantastical Comedy Drama Actions ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ Text on
Screen ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ Subjects ✓ ✓ ✓ ✓ ✓ ✓ ✓ Settings ✓ ✓ ✓ ✓ Clothing ✓ ✓
✓ ✓ ✓ ✓ Visual Aids ✓ ✓ ✓ Scene Changes ✓ ✓ ✓ Facial Expressions ✓ ✓ ✓ Visual
Effects ✓ ✓ ✓ Desired Details Equipment ✓ Ambient Sound ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓
Additional Resources ✓ ✓ ✓ ✓ ✓ ✓ Audio Cues ✓ ✓ ✓ ✓ ✓ Tactile Graphics ✓ ✓ ✓
Text to Braille ✓ ✓ 3D Models ✓ ✓ Background Music ✓ ✓ Output Modalities
Vibration ✓


# Modeling Generalized Specialist Approach To Train Quality Resilient Snapshot
Ensemble

Ghalib Ahmed Tahir , Chu Kiong Loo and Zongying Liu

###### Index Terms:

Deep learning, Ensemble learning, Data Augmentation, Malaysian Food database

##  I Abstract

Convolutional neural networks (CNNs) apply well with food image recognition
due to the ability to learn discriminative visual features. Nevertheless,
recognizing distorted images is challenging for existing CNNs. Hence, the
study modelled a generalized specialist approach to train a quality resilient
ensemble. The approach aids the models in the ensemble framework retain
general skills of recognizing clean images and shallow skills of classifying
noisy images with one deep expertise area on a particular distortion.
Subsequently, a novel data augmentation random quality mixup (RQMixUp) is
combined with snapshot ensembling to train G-Specialist. During each training
cycle of G-Specialist, a model is fine-tuned on the synthetic images generated
by RQMixup, intermixing clean and distorted images of a particular distortion
at a randomly chosen level. Resultantly, each snapshot in the ensemble gained
expertise on several distortion levels, with shallow skills on other quality
distortions. Next, the filter outputs from diverse experts were fused for
higher accuracy. The learning process has no additional cost due to a single
training process to train experts, compatible with a wide range of supervised
CNNs for transfer learning. Finally, the experimental analysis on three real-
world food and a Malaysian food database showed significant improvement for
distorted images with competitive classification performance on pristine food
images.

##  II INTRODUCTION AND CURRENT METHODOLOGIES

A vital aspect overlooked in existing methodology framework is the visual
quality of food images. In real-world computer vision applications, images
undergo various distortions, such as blur or additive noise during capturing
or transmission. Generally, current CNNs shows extreme sensitivity to these
image distortions [1][2][3][4]. Nevertheless, Dodge and Karam [6] stated that
the distorted images do not represent adversarial samples for CNNs but the
images cause a significant reduction in classification performance. Figure 1
shows the effect of image quality on the prediction performance of CNNs
trained on pristine food images. When the predicted label is correct, its
confidence decreases significantly with the increasing high distortion
severity, indicating that features learnt from the dataset of pristine images
are not invariant to image distortion nor useable for applications with
varying image quality.
Past studies addressed this challenge by fine-tuning with corrupted images at
all distortion levels [3][7], making the CNN model resilient to quality
distortions. Nonetheless, the method requires large distorted datasets at all
quality levels, which are not available and is time-consuming. Besides, fine-
tuning with low-quality images decreases the accuracy of the clean food photos
and increases the epochs to converge the model. Although other methods feature
quantization [9]], stability learning [8], and Deep correct [5] reduce the
impact of image quality on CNNs, requiring architectural level changes in the
pre-trained models. Similarly, stability learning and feature quantization
decrease the Top 1 accuracy, whereas Deep correct increases the floating-point
operations (FLOPS) during training and testing.

Figure 1: Effects of image quality on CNNs

The study extended the CNN-based approach by combining multiple networks with
a robust ensemble framework G-Specialist resilient against large variations in
quality distortions. Unlike image denoising methods, such as Block Matching 3D
denoising (BM3D) [16] and Non-locally Centralized Sparse Representation (NCSR)
[17] (which predict the known distortion intensity during training and
testing), the G-Specialist approach trains an ensemble for all distortion
levels simultaneously.

Accordingly, the study introduced a generalized specialist approach in deep
learning with each CNN in the ensemble carrying the general skills of
recognizing pristine images, specialist skills on a particular distortion
(Gaussian noise or blur) at all distortion levels, and shallow skills on other
distortions. Thus, a suitable solution to train generalized specialists during
the ensemble learning process is data augmentation.

The proposed RQMixup increases the quantifiable data value without increasing
the size of training data or change in the model architecture. Moreover, the
approach generates a synthetic sample coupled with snapshot ensembling [10],
creating a mixture of experts resilient to distortions and improves
classification performance on pristine images without extra training time or
increasing the training dataset.

The study has done experimentation on Malaysian Food database[ref], the
largest database yet to the authors’ best knowledge and experimented on other
three publicly available large food databases that will be publicly available
online [25].

##  III METHOD

###  III-A General Architecture

The proposed framework is based on the recent emergence of efficient neural
networks [18][19][20]. Furthermore, the study selected MobileNetV3 [18] as it
is resource-efficient, a crucial requirement of on-device computation for edge
devices. Next, RQMixup was combined with snapshot ensembling to train a
mixture of robust experts G-Specialist that are resilient to extreme image
distortions. Figure 2 shows the general architecture of the framework.

Figure 2: The proposed quality resilient ensemble for image recognition by
modeling generalized specialist approach during training

###  III-B Random MixUp Data Augmentation

Data involves driving and powering artificial intelligence (equivalent to fuel
for learning algorithms), a fundamental ingredient in supervised machine
learning problems with a quantified value in algorithmic classification tasks
measurable by the Data Shapley method [11] for supervised classification
tasks. Hence, the study developed a generic method that increases the
quantified value of data sets in algorithmic decisions to train robust models
resilient to different quality distortions without significant changes at the
architecture level.

Accordingly, the study chose RQMixup that generates a synthetic image to fine-
tune quality resilient deep learning models and used two copies of the same
training dataset, applying randomly chosen distortion levels to one copy of
the dataset during training. Next, raw vectors of cleaned and noisy images and
their corresponding labels were mixed using a strategy similar to mixup. For
example, raw vectors of a clean set of images in each batch are mixed by
applying blurring on the copy of the same set of the training dataset for the
Gaussian blur, randomly selecting the blur intensity (1, 3, 5, 7, 9, 11, 13,
and 15) for each batch. Consequently, fine-tuning the model with synthetic
images turns the model into an expert on all intensity levels of Gaussian blur
in the test samples. Equations 1 and 2 compute the mixup of feature-target
vectors of the pristine and corrupted images.

$\displaystyle\begin{array}[]{@{}l}x_{f}=\;\lambda
x_{i}+(1-\lambda)x_{n}\end{array}$ (1)
$\displaystyle\begin{array}[]{@{}l}y_{f}=\;\lambda
y_{i}+(1-\lambda)y_{n}\end{array}$ (2)

Specifically, $x_{n}$ is the raw input vector and $y_{n}$ is the one-hot label
encodings of the noisy image at the random distortion level, $x_{i}$ is the
raw input vector, and $y_{i}$ is the one-hot label encodings of the pristine
image. Meanwhile, Figure 3 shows samples of synthetic images generated to
fine-tune a model in a single cycle. Notably, the proposed method does not
increase the training images required to fine-tune the model

Figure 3: Samples of synthetic images to fine-tune a model in a single cycle.

###  III-C Snapshot Ensembling to Train G-Specialist by Modeling Generalized
Specialist Approach

Although the single model fine-tuned on synthetic images generated by the
RQMixup delivers all the distortion levels of particular noise, multiple types
of noises exist, such as Gaussian additive noise and blur. Therefore, the
study employed snapshot ensembling to train generalized specialists on various
distortions while inheriting the core benefit Train 1 Get M experts.

The process follows Loshchilov et al.’s [12] cyclic annealing schedule to
converge with multiple local minima and train M generalized specialists in a
single training process using synthetic images generated by RQMixup. The
learning rate was lowered at a very fast pace in each cycle, converging the
model towards local minima after a minimum of 32 epochs. The process was
repeated multiple times to obtain M generalized experts through multiple
convergences. The learning rate alpha is as follows:

$\displaystyle\alpha(t)\;=\;f(mod(t-1,\left[T/M\right]))$

Specifically, T is the total number of iterations and f is a monotonically
decreasing function. Thus, the training process of the whole framework was
split into M cycles whereby each cycle starts with a higher learning rate
annealed towards a smaller learning rate to converge towards local minima. A
larger learning rate at the start of the cycle provides the model energy to
escape a critical point while a smaller learning rate converges the model to
local minima. In the experiment, f is shifted cosine function following
Loshchilov et al. [12]

$\displaystyle\alpha(t)\;=\frac{\alpha_{0}}{2}\left(\cos\;\left({\textstyle\frac{\pi\;mod(t-1,\;\left[T/M\right])}{\left[T/M\right]}}\right)+1\right)$

The M-value is 2 for the Gaussian noise and blur. The total cycles were kept
equivalent to the number of distortions because during each process cycle, the
generalized specialist was trained on a particular distortion by applying
snapshot ensembling. The method enabled each snapshot to classify the clean
images with a particular distortion, such as Gaussian additive noise at all
intensity levels with confidence. Although the method is not an expert on
another distortion, such as blurring, it possesses shallow skills of
recognizing blur images. Hence, the data augmentation in the second cycle
mixes the clean images with a randomly chosen intensity level of blurring,
turning the second snapshot of the model to an expert on all distortion levels
in the images with shallow skills on other distortions.
As the model usually converges to local minima at the end of every cycle, a
snapshot of the model weights was taken before training the next generalized
experts on another distortion type. Furthermore, an M snapshot was used in the
final ensemble after training the model for M cycles. Besides making each
snapshot expert to the particular type of noise, the method also inherits the
benefits of ensembling for clean photos as each snapshot converges to
different local minima. Contrary to a single model, the ensemble of these
snapshots has lower test errors, as shown in Table II and III.
The ensemble prediction during testing is the average of M experts equivalent
to the total snapshots. Let I (w, h, and c) represent a pre-processed input
image of size w × h pixels, and c is the number of channels of the input
image. After feeding the image to the input conv2d of a snapshot in the
ensemble, g(x) will be the softmax score from a snapshot. Additionally, the
final output of the ensemble is an average of all models.

$\displaystyle g_{ensemble}=\;\frac{1}{M}\sum_{0}^{M-1}g(x)$

##  IV Experiment And Results

The experiments were performed on a server with 25 GB of random access memory
(RAM) equipped with 16 GB of the graphics processing unit (GPU). Besides, all
the experiments were performed using Python on three publically available
large food databases (FOOD101 [13], UECFOOD100 [14], UECFOOD256 [15]) and a
database of Malaysian foods. Table I summarizes the crucial attributes of the
datasets.

TABLE I: FOOD RECOGNITION DATASETS  Dataset |  Classes |  Instances |  Training/Testing
---|---|---|---
FOOD101 |  101 |  101330 |  75,750/25,250
UECFOOD100 |  100 |  14361 |  12864/1497
UECFOOD256 |  256 |  31148 |  28033/3115
Malaysian Food |  775 |  37198 |  33724/3474
Figure 4: The top row are the same class of Malaysian food while the bottom
row is random images of various classes of Malaysian food.

###  IV-A G-Specialist Analysis on Distorted Test Images

The G-Specialist was evaluated against the ensemble model on test images
distorted with the Gaussian noise and blurring, with a standard deviation of
10 to 100. Meanwhile, the Gaussian kernel during blurring of test images was
varied from 1 to 15. Figure 5 shows the test image results distorted with
Gaussian noise and blur for all datasets, denoting that the G-Specialist
significantly outperformed the ensemble fine-tuned on pristine images with a
significant difference for a higher noise level. Additionally, both models had
similar training costs.

The study also identified whether synthetic data adds value to the model by
employing a gradient explainer [21] to compute average Shapley values on test
images distorted with Gaussian noise and blur. Figure 6 shows that the noise
level increases while the test images average Shapley value decreases .
Nonetheless, the Shapley values for G-Specialists were significantly higher
than the ensemble fine-tuned with pristine images, suggesting improved
robustness by increasing the data value

Figure 5: The accuracy (%) on distorted data for Food101, UECFOOD256,
UECFOOD101, Malaysian Food and comparison with an ensemble of Mo- bileNetV3
after fine-tuning  Figure 6: The average Shapley value shows that as the noise
level increases, the average Shapley value on noisy images decreases

###  IV-B Results on Pristine Test Images

The results in Table II of Top 1 and Top 3 accuracy show that G-Specialist
outperformed all other variants of MobileNet including the single model of
MobileNet-V3 on the pristine test data. The study only reported the highest
performers, with more results seen in [24], [26], [27], [25], and [28].

###  IV-C Results on Malaysian Food Database

Malaysian Food database is comprised of 775 classes with 37,198 instances and
these classes were selected based on recommendations from an expert dietician.
The images were collected from Google search, Bing search, social media and
authors’ data using smartphones for capturing food images. Besides, the study
used the same experimental settings as the other three datasets. Figure 5
depicts that G-Specialist outperformed other methods on test images distorted
with Gaussian noise and blur. Meanwhile, Table III suggests that the
G-Specialist was better at recognizing food images than a single model and
other variants of MobileNet. Additionally, Figure 7 illustrates that the study
approach in considering relevant pixels for recognizing food images increases
the model transparency.

Figure 7: Red pixels show the SHAP values of the image positively contribute towards the class, whereas blue pixels show the pixels that reduce the probability of class. The top 3 outputs are explained with ranked_outputs = 3 TABLE II: ACCURACY (%) ON PRISTINE DATA FOR FOOD101, UECFOOD256, UECFOOD101 COMPARED WITH OTHER EFFICIENT NEURAL NETWORKS AFTER FINE-TUNING |  Top1 |  Top3
---|---|---
FOOD101
MobileNetV2 |  74.63 |  88.19
Ensemble Net [26] |  72.12 |  -
SSGAN [28] |  75.34 |  -
MobileNetV3 |  80.29 |  91.01
G-Specialist |  80.77 |  91.93
UECFOOD100
MobileNetV2 |  73.72 |  87.39
MobileNetV3 |  79.50 |  93.68
G-Specialist |  80.20 |  94.31
UECFOOD256
MobileNetV2 |  58.38 |  77.77
MobileNetV3 |  66.49 |  84.65
G-Specialist |  68.00 |  85.00
TABLE III: ACCURACY (%) ON PRISTINE DATA FOR MALAYSIAN FOOD COMPARED WITH OTHER EFFICIENT NEURAL NETWORKS AFTER FINE-TUNING  |  Top1 |  Top5
---|---|---
Malaysian Food
MobileNetV2 |  60.40 |  79.51
MobileNetV3 |  68.60 |  86.10
G-Specialist |  70.70 |  87.61

##  V Conclusion

Observably, efficient neural networks trained on pristine images had poor
classification performance for the distorted images affected by blur or
additive noise. The G-Specialist in the study addressed the issue by modeling
a generalized specialist approach to increase data value by generating
synthetic images on the fly during fine-tuning, helping each model in an
ensemble retain the general skill of recognizing clean images and deep
expertise on all intensity levels of a particular distortion. Significantly,
the study proved that each model in the ensemble individually contributed to a
better prediction for noisy test images by becoming an expert on a particular
noise, such as blur or additive noise.

Secondly, the study presented that the models in an ensemble converge to
different local minima improved classification performance on pristine data
compared to a single model. Most importantly, the study has experimented on
Malaysian Food database publicly available for further enrichment and
evaluation besided performing experiments on three publicly available real-
world image datasets. The results showed that the proposed G-Specialist
consistently outperformed other methods on pristine images and recognized
noisy test images with great precision than the other methodologies.

## References

  * [1] Lina J. Karam and Tong Zhu "Quality labeled faces in the wild (QLFW): a database for studying face recognition in real-world environments", Proc. SPIE 9394, Human Vision and Electronic Imaging XX, 93940B (17 March 2015); https://doi.org/10.1117/12.2080393
  * [2] S. Karahan, M. Kilinc Yildirum, K. Kirtac, F. S. Rende, G. Butun and H. K. Ekenel, "How Image Degradations Affect Deep CNN-Based Face Recognition?," 2016 International Conference of the Biometrics Special Interest Group (BIOSIG), 2016, pp. 1-5, doi: 10.1109/BIOSIG.2016.7736924.
  * [3] Rodner, Erik & Simon, Marcel & Fisher, Robert & Denzler, Joachim. (2016). Fine-grained Recognition in the Noisy Wild: Sensitivity Analysis of Convolutional Neural Networks Approaches.
  * [4] Deep object classification in low resolution LWIR imagery via transfer learning Abbott, R., Del Rincon, J. M., Connor, B., & Robertson, N. (2017). In Proceedings of the 5th IMA Conference on Mathematics in Defence
  * [5] Borkar, Tejas & Karam, Lina. (2017). DeepCorrect: Correcting DNN Models Against Image Distortions. IEEE Transactions on Image Processing. PP. 10.1109/TIP.2019.2924172.
  * [6] Dodge, Samuel & Karam, Lina. (2016). Understanding How Image Quality Affects Deep Neural Networks.
  * [7] Gazioğlu, Bilge & Kamasak, Mustafa. (2021). Effects of objects and image quality on melanoma classification using deep neural networks. Biomedical Signal Processing and Control. 67. 102530. 10.1016/j.bspc.2021.102530.
  * [8] Zheng, Stephan & Song, Yang & Leung, Thomas & Goodfellow, Ian. (2016). Improving the Robustness of Deep Neural Networks via Stability Training. 10.1109/CVPR.2016.485.
  * [9] Z. Sun, M. Ozay, Y. Zhang, X. Liu and T. Okatani, "Feature Quantization for Defending Against Distortion of Images," 2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2018, pp. 7957-7966, doi: 10.1109/CVPR.2018.00830.
  * [10] Huang, Gao & Li, Yixuan & Pleiss, Geoff & Liu, Zhuang & Hopcroft, John & Weinberger, Kilian. (2017). Snapshot Ensembles: Train 1, get M for free.
  * [11] Ghorbani, Amirata & Zou, James. (2019). Data Shapley: Equitable Valuation of Data for Machine Learning.
  * [12] Ilya Loshchilov and Frank Hutter. Sgdr: Stochastic gradient descent with restarts. arXiv preprint arXiv:1608.03983, 2016.
  * [13] Bossard, Lukas & Guillaumin, Matthieu & Van Gool, Luc., "Food-101 – Mining Discriminative Components with Random Forests," pp. 446-461, 2014.
  * [14] Y. Matsuda and K. Yanai, "Multiple-food recognition considering co-occurrence employing manifold ranking," in Proceedings of the 21st International Conference on Pattern Recognition (ICPR2012), Nov 2012.
  * [15] Automatic Expansion of a Food Image Dataset Leveraging Existing Categories with Domain Adaptation, "Proc. of ECCV Workshop on Transferring and Adapting Source Knowledge in Computer Vision (TASK-CV)","2014"
  * [16] Dabov, Kostadin & Foi, Alessandro & Katkovnik, Vladimir & Egiazarian, Karen. (2009). BM3D Image Denoising With Shape-Adaptive Principal Component Analysis. Proc. Workshop on Signal Processing with Adaptive Sparse Structured Representations (SPARS’09).
  * [17] Dong, Weisheng & Zhang, Lei & Shi, Guangming & li, Xin. (2012). Nonlocally Centralized Sparse Representation for Image Restoration. IEEE transactions on image processing : a publication of the IEEE Signal Processing Society. 22. 10.1109/TIP.2012.2235847.
  * [18] Howard, Andrew & Pang, Ruoming & Adam, Hartwig & Le, Quoc & Sandler, Mark & Chen, Bo & Wang, Weijun & Chen, Liang-Chieh & Tan, Mingxing & Chu, Grace & Vasudevan, Vijay & Zhu, Yukun. (2019). Searching for MobileNetV3. 1314-1324. 10.1109/ICCV.2019.00140.
  * [19] Howard, Andrew & Zhu, Menglong & Chen, Bo & Kalenichenko, Dmitry & Wang, Weijun & Weyand, Tobias & Andreetto, Marco & Adam, Hartwig. (2017). MobileNets: Efficient Convolutional Neural Networks for Mobile Vision Applications.
  * [20] Sandler, Mark & Howard, Andrew & Zhu, Menglong & Zhmoginov, Andrey & Chen, Liang-Chieh. (2018). MobileNetV2: Inverted Residuals and Linear Bottlenecks. 4510-4520. 10.1109/CVPR.2018.00474.
  * [21] Lundberg, Scott & Lee, Su-In. (2017). A Unified Approach to Interpreting Model Predictions.
  * [22] Ghalib Ahmed Tahir, https://github.com/ghalib2021/-G-Specialist
  * [23] Patrick McAllister, Huiru Zheng, Raymond Bond, Anne Moorhead,Combining deep residual neural network features with supervised machine learning algorithms to classify diverse food image datasets,Computers in Biology and Medicine,2018
  * [24] M. Niki & P. Claudio & M. Christian., "A supervised extreme learning committee for food recognition",Computer Vision and Image Understanding 2016\.
  * [25] K. Yanai and Y. Kawano, "Food image recognition using deep convolutional network with pre-training and fine-tuning,” inProc. IEEE Int. Conf.Multimedia Expo. Workshops, 2015, pp. 1–6
  * [26] Pandey,A.Deepthi,B.Mandal,and N.B.Puhan,“Foodnet:Recognizing foods using ensemble of deep networks,”IEEE Signal Process. Lett.,vol. 24, no. 12, pp. 1758–1762, Dec. 2017.
  * [27] Liu, Chang et al. "A New Deep Learning-Based Food Recognition System for Dietary Assessment on An Edge Computing Service Infrastructure," IEEE Transactions on Services Computing, pp. 1-1, 2017.
  * [28] Mandal, Bappaditya & Puhan, Niladri & Verma, Avijit. (2018). Deep Convolutional Generative Adversarial Network Based Food Recognition Using Partially Labeled Data. IEEE Sensors Letters. PP. 1-1. 10.1109/LSENS.2018.2886427.