# Evolution of electronic band reconstruction in thickness-controlled
perovskite SrRuO3 thin films

Byungmin Sohn  Center for Correlated Electron Systems, Institute for Basic
Science, Seoul 08826, Korea  Department of Physics and Astronomy, Seoul
National University, Seoul 08826, Korea  Changyoung Kim<EMAIL_ADDRESS>Center for Correlated Electron Systems, Institute for Basic Science, Seoul
08826, Korea  Department of Physics and Astronomy, Seoul National University,
Seoul 08826, Korea

(August 27, 2024)

###### Abstract

Transition metal perovskite oxides display a variety of emergent phenomena
which are tunable by tailoring the oxygen octahedral rotation. SrRuO3, a
ferromagnetic perovskite oxide, is well-known to have various atomic
structures and octahedral rotations when grown as thin films. However, how the
electronic structure changes with the film thickness has been hardly studied.
Here, by using angle-resolved photoemission spectroscopy and electron
diffraction techniques, we study the electronic structure of SrRuO3 thin films
as a function of the film thickness. Different reconstructed electronic
structures and spectral weights are observed for films with various
thicknesses. We suggest that octahedral rotations on the surface can be
qualitatively estimated via comparison of intensities of different bands. Our
observation and methodology shed light on how structural variation and
transition may be understood in terms of photoemission spectroscopy data.

##  I introduction

Transition metal oxides (TMO) have been intensively investigated over the past
decades owing to their novel phenomena such as high-temperature
superconductivity, multiferroicity, Mott transition, and magnetism [1, 2, 3,
4]. These phenomena are often understood to emerge from the strongly
correlated $d$ electrons in which the fundamental degrees of freedom – spin,
orbital, charge, and lattice – are entangled. Especially, in Ruddlesden-popper
(RP) oxide perovskites, which have a chemical formula of An+1MnO3n+1 ($n$ = 1,
2, 3 $\cdots$), transition metal $d$ electrons are combined with oxygen $p$
electrons, forming MO6 oxygen octahedron electronic structure. Controlling the
structure of oxygen octahedra determines the metallicity and magnetism in many
RP oxides such as Sr2RuO4 [5], Sr2IrO4 [6], and SrRuO3 (SRO) [7, 8].

In particular, MO6 octahedron rotation significantly alters the electronic
structure and may induce phase transitions, often accompanying changes in
physical properties. A prime example is a well-known ferromagnetic metallic
TMO perovskite SRO for which various characteristics appear with RuO6
octahedral rotations; a large octahedral rotation angle favors an insulating
antiferromagnetic phase whereas a small octahedral rotation prefers a metallic
ground state with ferromagnetism [9, 10, 11]. As such, numerous attempts to
control the MO6 octahedron angle and characteristics of TMO systems have been
reported [10, 12, 13, 14]. All these show that MO6 octahedral rotation can
play an important role in modifying the intrinsic properties. Accordingly,
finding suitable experimental methods to control oxygen octahedra in TMOs has
been of great interest.

One of the well-adopted experimental methods for tuning the MO6 octahedra is
growing TMOs in thin films. Due to the strain and interfacial effects induced
by the substrate, thin films may have octahedral rotation different from that
of the bulk [15]. Since thin films tend to have similar characteristics with
bulk when the thin films are thick enough, variation in the octahedral
rotation (and, crystal symmetry) as a function of the thickness of thin films
has been studied [16, 17, 18]. However, most studies were focused on atomic
structures through X-ray diffraction analysis, whereas only a few studies of
thickness-dependent electronic structures have been reported so far.

Figure 1: Reflection high-energy electron diffraction (RHEED) and low-energy
electron diffraction (LEED) patterns of various atomic structures. Schematics
of (a) cubic, (d) tetragonal, and (g) orthorhombic structures of perovskites.
RHEED patterns of (c) SrTiO3 (STO) (001) substrate, (f) 4 unit-cell (uc), and
(i) 50 uc SrRuO3 (SRO) thin films on grown on STO substrates are given. RHEED
patterns are measured along the (100)pc crystalline direction. [0
$\frac{1}{2}$] peaks marked with orange arrows appear due to the 2 $\times$ 2
reconstruction. LEED schematics and patterns of (b) STO (001) substrate, (e) 4
uc, and (h) 50 uc SRO thin films on STO substrates are plotted. LEED
measurements were done with an electron energy of 150 eV at 300 K. The 4 uc
and 50 uc SROs show $\sqrt{2}$ $\times$ $\sqrt{2}$ and 2 $\times$2
reconstruction, respectively (marked with white arrows). When the electron
($e^{-}$) beam comes along the (100)pc direction, [0 $\frac{1}{2}$] peaks are
observed in the RHEED pattern of orthorhombic perovskite structure (i) due to
the 2 $\times$ 2 reconstruction. f, i, and h are adapted from ref. [19].

Here, we investigate atomic and electronic structures of SRO thin films grown
on SrTiO3 (STO) (001) substrates as a function of the thickness with in-situ
reflection high-energy electron diffraction (RHEED), low-energy electron
diffraction (LEED) and angle-resolved photoemission spectroscopy (ARPES).
Atomic structural transitions appear with thickness variation, and
reconstructed electronic structures are observed in ARPES data. Based on our
simple theoretical model, we claim that the intensity of the reconstructed
band structure is intimately related to oxygen octahedral rotation.

##  II Methods

SRO thin films with thicknesses of 2.5, 4, 5, 10, 30, and 50 unit-cell (uc)
were grown on TiO2-terminated STO (001) substrates with RHEED-monitored pulsed
laser deposition (PLD). Here, a 2.5 uc film consists of 50 $\%$ of 2-uc-thick
SRO and 50 $\%$ of 3-uc-thick SRO. Since SRO is grown by step-flow mode, the
fractional thickness can be intentionally obtained by stopping the deposition
in the middle of the growth process [19, 20]. Before the growth, STO
substrates were prepared with deionized water etching and then in-situ pre-
annealing at 1,070 ∘C for 30 minutes under the oxygen partial pressure of 5
$\times$ 10-6 torr. SRO thin films were grown under the oxygen partial
pressure of 100 mTorr at 700 ∘C. A KrF excimer laser (248 nm) irradiated an
SRO target with a fluence of 1-2 J/cm2 and a repetition rate of 2 Hz. Samples
were in-situ transferred for ARPES and LEED measurements. ARPES measurements
were performed at 10 K equipped with a hemispherical electron analyzer and a
He discharge lamp. He-I$\alpha$ light with the photon energy of 21.2 eV was
used. Before ARPES measurements, samples were annealed at 550 ∘C to remove
surface contaminants.

##  III Results

We measured RHEED and LEED of a bare STO (001) substrate, 4 uc and 50 uc SRO
on STO (001) substrates at room temperature. Here, electrons come through the
(001)pc and (100)pc direction for LEED and RHEED, respectively. Since STO has
a cubic structure at room temperature (with a symmetry group of Pm-3m and a
glazer notation of a0a0a0) [21], oxygen octahedra are not expected to be
rotated in an STO substrate at all (Fig. 1a). Thus, the LEED pattern of STO
substrates shows a square-shaped Brillouin zone without any reconstruction
(Fig. 1b). Bright three spots from 0$\overline{1}$, 00, and 01 peaks, crossing
an Ewald sphere, are observed in the RHEED pattern (Fig. 1c).

On the other hand, the 4 uc SRO on an STO substrate is expected to have a
tetragonal (I4/mcm, a0a0c$-$) structure with in-plane octahedral rotation
(Fig. 1d) [22], resulting in in-plane unit-cell doubling with the $\sqrt{2}$
$\times$$\sqrt{2}$ periodicity. $\sqrt{2}$ $\times$$\sqrt{2}$ reconstruction
spots appear in LEED patterns (white arrows in Fig. 1e). Since an Ewald sphere
is not crossing the $\sqrt{2}$ $\times$$\sqrt{2}$ spots with our RHEED
geometry, the same RHEED patterns are observed for the bare STO and 4 uc SRO
(Fig. 1f).

It has been reported that a thickness-dependent structural phase transition
(from tetragonal to orthorhombic phases) occurs in SRO thin films on an STO
substrate when SRO films become thicker than 19 uc [17]. Thus, the 50 uc SRO
film is expected to be orthorhombic (Pbnm, a-b+a-) with in-plane and out-of-
plane octahedral rotation (Fig. 1g), accompanied by in-plane unit-cell
doubling with 2 $\times$ 2 reconstruction. Our LEED pattern clearly shows 2
$\times$ 2 spots (orange and white arrows in Fig. 1h). Contrary to the
tetragonal case, an Ewald sphere is crossing the 2 $\times$ 2 reconstruction
spots, which are observed in the RHEED pattern (marked with orange arrows in
Fig. 1i).

Figure 2: Observation of $\sqrt{2}$ $\times$$\sqrt{2}$ and 2 $\times$ 2
reconstructed bands in SRO thin films. a A schematic of $n$ uc SRO thin films
on an STO (001) substrate. b,c RHEED patterns of 10 and 30 uc SRO thin films
on STO substrates. [0 $\frac{1}{2}$] peaks are marked with orange arrows. d-f
Fermi surface (FS) schematics of cubic (1 $\times$ 1, d), tetragonal
($\sqrt{2}$ $\times$$\sqrt{2}$, e), and orthorhombic (2 $\times$ 2, f) SRO
thin films. In d, $\beta_{\rm min}$ is represented with black solid lines,
whereas other bands are plotted with gray solid lines. In e and f,
$\alpha^{{}^{\prime}}_{\rm maj}$ and $\beta^{{}^{\prime\prime}}_{\rm min}$
bands are represented with red and blue solid lines, respectively. To make it
easy to understand, we only plot specific bands in e and f. g,h FSs of 10 and
30 uc SRO thin films. i,j Energy-momentum cuts of 10 and 30 uc SRO thin films
along the $\Gamma$-X direction ($k_{x}$ = 0). Momentum distribution curves
(MDCs) at $E$ = $E_{F}$ \- 30 meV, where $E_{F}$ is the Fermi level, are also
plotted (red curves). The open (closed) inverted triangle represents a peak
from $\beta_{\rm min}$ ($\beta^{{}^{\prime\prime}}_{\rm min}$) band. Red
arrows in h and j indicate $\beta^{{}^{\prime\prime}}_{\rm min}$ band.

We next compare electronic structures of tetragonal and orthorhombic SRO thin
films by performing ARPES on 10 uc and 30 uc SRO thin films (Fig. 2a). To
investigate the structural phases of each film, RHEED patterns were observed
after the growth process (Fig. 2b and 2c). The 30 uc SRO thin film shows 2
$\times$ 2 reconstruction represented by orange arrows in Fig. 2c, which
supports that 10 uc (30 uc) SRO is tetragonal (orthorhombic) consistent with a
previous report [17].

Because Ru 4$d$ states cross the Fermi level ($E_{F}$), the electronic
structure of SRO is often described within the $t_{2g}$ manifold with four
electrons [23]. The $t_{2g}$ manifold consists of six states ($d_{xz}$,
$d_{yz}$, $d_{xy}$ with majority (maj) and minority (min) spins) due to the
ferromagnetism and six bands are expected to cross $E_{F}$. Based on previous
reports [24, 19, 25], we present FS schematics as depicted in Figs. 2d-f.
Without octahedral rotation (Fig. 2d), two hole-like spin-resolved bands
($\alpha_{\rm maj}$, $\alpha_{\rm min}$) and four electron-like spin-resolved
bands ($\beta_{\rm maj}$, $\beta_{\rm min}$, $\gamma_{\rm maj}$, $\gamma_{\rm
min}$) are expected to cross $E_{F}$. In tetragonal SRO films, in-plane
octahedral rotation induces $\sqrt{2}$ $\times$$\sqrt{2}$ reconstruction and
the Brillouin zone is reduced to the one indicated by a red-dotted square
(Fig. 2e). When SRO films become orthorhombic and have both in-plane and out-
of-plane octahedral rotations, the Brillouin zone is further reduced to the
one shown by a blue-dotted square (Fig. 2f) due to the 2 $\times$ 2
reconstruction.

Figures 2g and 2h show measured Fermi surfaces (FSs) of 10 uc and 30 uc SRO
films. Note that due to matrix element effects only part of the band structure
is often visible in ARPES measurements [26]. We noticed that the $\beta_{\rm
min}$ band is clearly observed in both FS maps (corresponds to the black
pocket in Fig. 2d). Note that the Fermi momentum ($k_{F}$) of the $\beta_{\rm
min}$ band is almost the same ($k_{F}$ = 0.51 ${\rm\AA^{-1}}$ along $k_{y}$ =
0) in the two maps. Strong spectral weight is observed at the $\Gamma$ point
which originates from the $\alpha^{\prime}_{\rm maj}$ band observed due to
$\sqrt{2}$ $\times$$\sqrt{2}$ reconstruction (the red pocket in Fig. 2e) [19].
In Fig. 2h, interestingly, additional spectral feature (marked with a red
arrow) appears, which is not observed for the 10 uc SRO FS.

To examine the additional spectra observed in the 30 uc SRO FS, energy-
momentum high-symmetry cuts along $k_{x}$ = 0 (red solid lines in Figs. 2g and
2h) are presented in Figs. 2i and 2j. Strong spectra from
$\alpha^{\prime}_{\rm maj}$ and the $\beta_{\rm min}$ band are clearly
observed in both high-symmetry cuts, whereas an additional dispersive band
(marked with a red arrow) is observed in the 30 uc SRO. Momentum distribution
curve (MDC) at $E$ = $E_{F}$ \- 30 meV is plotted for each figure, showing an
additional peak (marked with a closed inverted triangle) appears in 30 uc SRO.
As the additional dispersive band and $\beta_{\rm min}$ band are symmetric
with respect to $k_{y}$ = 0.4 ${\rm\AA^{-1}}$, we conclude the additional
dispersive band is $\beta^{\prime\prime}_{\rm min}$ band, the blue pocket in
Fig. 2f, induced by 2 $\times$ 2 reconstruction.

Figure 3: Observation of thickness-dependent intensity of $\sqrt{2}$
$\times$$\sqrt{2}$ reconstructed bands. a-d FSs of $n$ uc ($n$ = 2.5, 4, 5,
and 10) SRO thin films grown on STO (001) substrates. (Inset) RHEED patterns
measured along the (100)pc direction. e,f Energy momentum cuts of 2.5 and 10
uc SRO thin films along the $\Gamma$-X direction ($k_{y}$ = 0). g MDCs along
$k_{y}$ = 0 at $E$ = $E_{F}$ as a function of the thickness. Peaks from
$\beta_{\rm min}$ band are marked with a blue dotted line. All MDCs are
normalized with the intensity at the $\Gamma$ point. b is adapted from ref.
[19].

We also investigate how the electronic structures vary as a function of the
thickness for SRO films within the tetragonal structure. Figures 3a-d show FSs
of 2.5, 4, 5, and 10 uc SRO films. Based on the RHEED patterns in inset
figures and results in previous reports [27, 17, 22], we determine that all of
these SRO films are in tetragonal structure. We note that, regardless of the
thickness, the fermiology of the films remains the same. Figures 3e and 3f
show energy-momentum band dispersion of 2.5 and 10 uc SRO films along $k_{y}$
= 0. A spectral weight near the $\Gamma$ point is strongly observed in the 10
uc SRO film, whereas the spectral weight is suppressed in the 2.5 uc SRO film.

We plot MDCs of these tetragonal SRO films along $k_{y}$ = 0 at $E_{F}$ (Fig.
3g). Three peaks are seen in each MDC. As we discussed above for Fig. 2, the
spectral weight at the $\Gamma$ point ($k_{x}$ = 0.51 ${\rm\AA^{-1}}$) stems
from $\alpha^{\prime}_{\rm maj}$ ($\beta_{\rm min}$). Notably, the peak ratio
between $\alpha^{\prime}_{\rm maj}$ and $\beta_{\rm min}$ becomes larger as
the thickness increases. Meanwhile, the position of the $\beta_{\rm min}$ peak
remains almost the same at $k_{x}$ = 0.51 ${\rm\AA^{-1}}$, regardless of the
thickness.

##  IV Discussion

Now, we wish to discuss our observations. ARPES results from 2.5, 4, 5, and 10
uc tetragonal SRO films do not show significant variation in the band
structure but the spectral weight intensity differs. Since electronic
structures are highly sensitive to the structure of MO6 octahedra in SRO
systems [28, 29], FS should appear different if MO6 octahedral rotations are
different. In fact, ARPES is highly surface-sensitive; less than 1 nm ($\sim$
2.5 uc of SRO layers) is expected to be probed with 21.2 eV photon [30]. Thus,
our observations indicate that the atomic structure near the surface may
remian almost the same regardless of the thickness, resulting in the same
fermiology and $k_{F}$.

However, if the atomic structure is exactly the same, the different intensity
ratios in MDCs (Fig. 3g) cannot be explained. Since SRO films we studied are
ultrathin ones and the lattice constant mismatch between bulk STO and SRO is
0.06 % [7], SRO films are expected to be fully strained, i.e. the in-plane
lattice constant of the SRO films is the same as that of STO substrate [22,
17, 12]. Instead, we speculate that the MO6 octahedral rotation angle near the
surface layers is different. Namely, we claim that the octahedral rotations
are different among the SRO films we studied. The difference in the octahedral
rotation angle must be responsible for the difference in the intensity ratio.

Figure 4: Intensity variation of reconstructed bands with MO6 octahedral
rotation. a Schematic of a single tetragonal SRO layer with an in-plane
octahedral rotation angle, $\theta$. b Schematic of $d_{xz}$ orbitals with and
without in-plane oxygen octahedral rotation. c,d Schematics of $\pi$\- and
$\delta$-bondings between two $d_{xz}$ orbitals. e Intensity ratio between the
$\sqrt{2}$ $\times$$\sqrt{2}$ reconstructed and original band as a function of
the thickness. $I_{\alpha^{\prime},{\rm maj}}$ and $I_{\beta,{\rm min}}$ are
defined as $E_{F}$ MDC intensity at the $\Gamma$ point and the peak position
of $\beta_{\rm min}$ band, respectively. f $r$ vs $\theta$ plot where $r$ is
defined as $\theta^{2}$($t_{\delta}$/$t_{\pi}$)2. $t_{\delta}$ and $t_{\pi}$
are hopping integrals of $\delta$\- and $\pi$-bonding electrons between two
$d_{xz}$ orbitals.

We attempt to explain the thickness-dependent suppression of the spectral
weight near the $\Gamma$ point with a simple toy model. Especially, we discuss
the intensity suppression along $k_{x}$ = 0 given in Fig. 3g. The spectral
weight near the $\Gamma$ point stems from $\alpha^{\prime}_{\rm maj}$; thus,
$d_{xz}$ and $d_{yz}$ orbitals mainly contribute to the spectral weight [19].
Since the electron hopping integral between $d_{yz}$ orbitals along the x-axis
is approximately ten times smaller than that of between $d_{xz}$ orbitals
[31], we only consider $d_{xz}$ orbitals in our model. Figure 4a is a
schematic of a single in-plane tetragonal SRO layer. We define $\theta$ as the
in-plane octahedral rotation angle. Chains of $d_{xz}$ orbitals in an SRO
layer are shown in Fig. 4b for cubic ($\theta$ = 0) and tetragonal ($\theta$ >
0) structure. When $\theta$ is zero, $d_{xz}$ orbitals are not rotated, so
electrons in $d_{xz}$ hop only through a $\pi$ bonding (Fig. 4c). On the other
hand, $d_{xz}$ orbitals rotate when $\theta$ > 0 and electrons hop through
both $\pi$ and $\delta$ bondings (Fig. 4d). When the $d_{xz}$ chain is equally
spaced but $d_{xz}$ orbitals are distinct and rotated, the intensity ratio
between the folded and original band (or, the low and high intensity band)
along the $x$-axis is expressed as

$\frac{I^{\prime}}{I}=\frac{(\braket{k_{f}}{d_{xz,\theta}}-\braket{k_{f}}{d_{xz,-\theta}})^{2}}{(\braket{k_{f}}{d_{xz,\theta}}+\braket{k_{f}}{d_{xz,-\theta}})^{2}},$
(1)

where $I^{\prime}$, $I$, and $\ket{k_{f}}$ denote the intensity of a folded
band, the intensity of the original band, and a single plane wave component
with a momentum $k_{f}$, respectively [26]. $\ket{d_{xz,\theta}}$ and
$\ket{d_{xz,-\theta}}$ are a $d_{xz}$ orbital rotated by $\theta$ in the
clockwise and counter-clockwise direction, respectively. With an assumption
that $\theta$ is small,
$\braket{k_{f}}{d_{xz,\pm\theta}}\simeq\braket{k_{f}}{d_{xz}\pm\theta
d_{yz}}$, where $d_{xz}$ ($d_{yz}$) represents $d_{xz,\pm\theta}$
($d_{yz,\pm\theta}$) with $\theta$ = 0. Thus, the intensity ratio between the
bands with the $d_{xz}$ orbital component is expressed as

$\frac{I^{\prime}}{I}\simeq(\frac{\braket{k_{f}}{\theta
d_{yz}}}{\braket{k_{f}}{d_{xz}}})^{2}\propto\theta^{2}(\frac{t_{\delta}}{t_{\pi}})^{2},$
(2)

where $t_{\pi}$ ($t_{\delta}$) denotes an electron hopping integral between
the same $d_{xz}$ orbitals through $\pi$ ($\delta$) bonding. Note that along
the $y$-axis, the electron hopping integral between $d_{xz}$ orbitals is
approximately ten times smaller than that of between $d_{yz}$ orbitals [31].
Thus, the intensity suppression can be well explained with chains of $d_{yz}$
orbitals instead of those of $d_{xz}$ orbitals. For the y-axis case, the
$d_{yz}$ ($d_{xz}$) orbital is substituted by $d_{xz}$ ($d_{yz}$) in
equations.

We compare our experimental observation with the theoretical equation. The
intensity ratio between the $\alpha^{\prime}_{\rm maj}$ and $\beta_{\rm min}$
band spectral weights, $I_{\alpha^{\prime},{\rm maj}}$/$I_{\beta,{\rm min}}$,
is plotted in Fig. 4e, and $r$ = $(\theta t_{\delta}/t_{\pi})^{2}$ is shown as
a function of $\theta$ in Fig. 4f. We assume a value of 0.1 for
$t_{\delta}/t_{\pi}$ based on the reported Slater-Koster parameters in bulk
SRO [31]. As the SRO thin film becomes thicker, $I_{\alpha^{\prime},{\rm
maj}}$/$I_{\beta,{\rm min}}$ increases. Similarly, $r$ also increases with
$\theta$. Since $I_{\alpha^{\prime},{\rm maj}}$/$I_{\beta,{\rm min}}$ is
proportional to $r$, we claim that the large (small) intensity ratio stems
from the large (small) octahedral rotation angle.

It has been reported that octahedral rotations in ultrathin films highly
depend on the thickness due to the interfacial effects; the MO6 octahedra of
films tend to adopt the same octahedral rotations of the substrate [32, 33],
which relaxes as it propagates through the film. Since STO has a smaller
octahedral rotation angle ($\theta$ $\sim$ 2 deg at 10 K [34, 35, 36]) than
that of bulk SRO ($\theta$ $\sim$ 6 deg at room temperature [37]), we expect
that thicker SRO films have a larger rotation angle on the surface than
thinner films. We suggest that our observation of the thickness-dependent
intensity ratio variation can be an indication of different octahedral
rotation on the surface of SRO ultrathin films. Further studies may be desired
to address the issue in a more quantitative manner.

##  V Conclusion

In summary, we investigated how the electronic structure on the surface of SRO
thin film changes with the atomic reconstruction through thin-film thickness
dependent studies. Differently reconstructed band structures were observed in
orthorhombic (thick) and tetragonal (thin) SRO films. As for SRO ultrathin
films (n $\leqq$ 10 uc), they show almost thickness independent Fermi surface
topology but exhibit different spectral weight intensities. Based on our
theoretical model, we conclude that the rotation angle of the surface
octahedra becomes smaller as the thickness is reduced in the ultrathin limit.
Our observation and methodology may shed light on how structural change and
transition can be understood in the view of photoemission spectroscopy and
electronic structures.

## Acknowledgements

We gratefully acknowledge insightful discussions with Younsik Kim. This work
is supported by the Institute for Basic science in Korea (Grant No.
IBS-R009-G2). CK acknowledges the supported by the National Research
Foundation of Korea(NRF) grant funded by the Korea government(MSIT). (No.
2022R1A3B1077234)

## References

  * Tokura and Nagaosa [2000] Y. Tokura and N. Nagaosa, Science 288, 462 (2000).
  * Damascelli  _et al._ [2003] A. Damascelli, Z. Hussain, and Z.-X. Shen, Rev. Mod. Phys. 75, 473 (2003).
  * Spaldin and Ramesh [2019] N. A. Spaldin and R. Ramesh, Nat. Mater. 18, 203 (2019).
  * Ahn  _et al._ [2021] C. Ahn, A. Cavalleri, A. Georges, S. Ismail-Beigi, A. J. Millis, and J.-M. Triscone, Nat. Mater. 20, 1462 (2021).
  * Anisimov  _et al._ [2002] V. I. Anisimov, I. A. Nekrasov, D. E. Kondakov, T. M. Rice, and M. Sigrist, Eur. Phys. J. B 25, 191 (2002).
  * Kim  _et al._ [2009] B. J. Kim, H. Ohsumi, T. Komesu, S. Sakai, T. Morita, H. Takagi, and T.-h. Arima, Science 323, 1329 (2009).
  * Koster  _et al._ [2012] G. Koster, L. Klein, W. Siemons, G. Rijnders, J. S. Dodge, C.-B. Eom, D. H. Blank, and M. R. Beasley, Rev. Mod. Phys. 84, 253 (2012).
  * [8] J. R. Kim, B. Sohn, H. J. Lee, S. Lee, E. K. Ko, S. Hahn, S. Lee, Y. Kim, D. Kim, H. J. Kim, _et al._ , arXiv:2203.04244 [cond-mat.str-el] .
  * Ahn  _et al._ [1999] J. S. Ahn, J. Bak, H. S. Choi, T. W. Noh, J. E. Han, Y. Bang, J. H. Cho, and Q. X. Jia, Phys. Rev. Lett. 82, 5321 (1999).
  * Jin  _et al._ [2008] C.-Q. Jin, J.-S. Zhou, J. B. Goodenough, Q. Q. Liu, J. G. Zhao, L. X. Yang, Y. Yu, R. C. Yu, T. Katsura, A. Shatskiy, _et al._ , PNAS 105, 7115 (2008).
  * Ali  _et al._ [2022] Z. Ali, Z. Wang, A. O’Hara, M. Saghayezhian, D. Shin, Y. Zhu, S. T. Pantelides, and J. Zhang, Phys. Rev. B 105, 054429 (2022).
  * Vailionis  _et al._ [2008] A. Vailionis, W. Siemons, and G. Koster, Appl. Phys. Lett. 93, 051909 (2008).
  * Lu  _et al._ [2015] W. Lu, W. Song, P. Yang, J. Ding, G. M. Chow, and J. Chen, Sci. Rep. 5, 1 (2015).
  * Balachandran and Rondinelli [2013] P. V. Balachandran and J. M. Rondinelli, Phys. Rev. B 88, 054101 (2013).
  * Hwang  _et al._ [2012] H. Y. Hwang, Y. Iwasa, M. Kawasaki, B. Keimer, N. Nagaosa, and Y. Tokura, Nat. Mater. 11, 103 (2012).
  * Siwakoti  _et al._ [2021] P. Siwakoti, Z. Wang, M. Saghayezhian, D. Howe, Z. Ali, Y. Zhu, and J. Zhang, Phys. Rev. Mater. 5, 114409 (2021).
  * Chang  _et al._ [2011] S. H. Chang, Y. J. Chang, S. Jang, D. Jeong, C. Jung, Y.-J. Kim, J.-S. Chung, and T. Noh, Phys. Rev. B 84, 104101 (2011).
  * Roh  _et al._ [2021] C. J. Roh, J. R. Kim, S. Park, Y. J. Shin, B.-J. Yang, T. W. Noh, and J. S. Lee, Appl. Surf. Sci. 553, 149574 (2021).
  * Sohn  _et al._ [2021a] B. Sohn, E. Lee, S. Y. Park, W. Kyung, J. Hwang, J. D. Denlinger, M. Kim, D. Kim, B. Kim, H. Ryu, _et al._ , Nat. Mater. 20, 1643 (2021a).
  * Wang  _et al._ [2020] L. Wang, Q. Feng, H. G. Lee, E. K. Ko, Q. Lu, and T. W. Noh, Nano Lett. 20, 2468 (2020).
  * Rimai and DeMars [1962] L. Rimai and G. DeMars, Phys. Rev. 127, 702 (1962).
  * Sohn  _et al._ [2021b] B. Sohn, B. Kim, S. Y. Park, H. Y. Choi, J. Y. Moon, T. Choi, Y. J. Choi, H. Zhou, J. W. Choi, A. Bombardi, _et al._ , Phys. Rev. Res. 3, 023232 (2021b).
  * Chen  _et al._ [2013] Y. Chen, D. L. Bergman, and A. A. Burkov, Phys. Rev. B 88, 125110 (2013).
  * Shai  _et al._ [2013] D. E. Shai, C. Adamo, D. W. Shen, C. M. Brooks, J. W. Harter, E. J. Monkman, B. Burganov, D. G. Schlom, and K. M. Shen, Phys. Rev. Lett. 110, 087004 (2013).
  * Hahn  _et al._ [2021] S. Hahn, B. Sohn, M. Kim, J. R. Kim, S. Huh, Y. Kim, W. Kyung, M. Kim, D. Kim, Y. Kim, _et al._ , Phys. Rev. Lett. 127, 256401 (2021).
  * Moser [2017] S. Moser, J. Electron Spectrosc. Relat. Phenom. 214, 29 (2017).
  * Sohn  _et al._ [2021c] B. Sohn, J. R. Kim, C. H. Kim, S. Lee, S. Hahn, Y. Kim, S. Huh, D. Kim, Y. Kim, W. Kyung, _et al._ , Nat. Commun. 12, 6171 (2021c).
  * Kim  _et al._ [2022] M. Kim, J. Kwon, C. H. Kim, Y. Kim, D. Chung, H. Ryu, J. Jung, B. S. Kim, D. Song, J. D. Denlinger, _et al._ , npj Quantum Mater. 7, 1 (2022).
  * Ko  _et al._ [2007] E. Ko, B. J. Kim, C. Kim, and H. J. Choi, Phys. Rev. Lett. 98, 226401 (2007).
  * Seah and Dench [1979] M. P. Seah and W. A. Dench, Surf. Interface Anal. 1, 2 (1979).
  * Mazin  _et al._ [2000] I. I. Mazin, D. A. Papaconstantopoulos, and D. Singh, Phys. Rev. B 61, 5223 (2000).
  * Fowlie  _et al._ [2019] J. Fowlie, C. Lichtensteiger, M. Gibert, H. Meley, P. Willmott, and J.-M. Triscone, Nano Lett. 19, 4188 (2019).
  * Liao  _et al._ [2016] Z. Liao, M. Huijben, Z. Zhong, N. Gauquelin, S. Macke, R. Green, S. Van Aert, J. Verbeeck, G. Van Tendeloo, K. Held, _et al._ , Nat. Mater. 15, 425 (2016).
  * Unoki and Sakudo [1967] H. Unoki and T. Sakudo, J. Phys. Soc. Japan 23, 546 (1967).
  * Müller  _et al._ [1968] K. A. Müller, W. Berlinger, and F. Waldner, Phys. Rev. Lett. 21, 814 (1968).
  * Tsuda and Tanaka [1995] K. Tsuda and M. Tanaka, Acta Crystallogr. A: Found. Adv. 51, 7 (1995).
  * Gao  _et al._ [2016] R. Gao, Y. Dong, H. Xu, H. Zhou, Y. Yuan, V. Gopalan, C. Gao, D. D. Fong, Z. Chen, Z. Luo, _et al._ , ACS Appl. Mater. Interfaces 8, 14871 (2016).


# Rapid Changes in Synchronizability in Conductance-based Neuronal Networks
with Conductance-based Coupling

Wilten Nicola  Hotchkiss Brain Institute  University of Calgary, Faculty of
Science, Department of Physics and Astronomy

###### Abstract

Real neurons connect to each other non-randomly. These connectivity graphs can
potentially impact the ability of networks to synchronize, along with the
dynamics of neurons and the dynamics of their connections. How the
connectivity of networks of conductance-based neuron models like the classical
Hodgkin-Huxley model, or the Morris-Lecar model, impacts synchronizability
remains unknown. One powerful tool to resolve the synchronizability of these
networks is the Master Stability Function (MSF). Here, we apply and extend the
MSF approach to networks of Morris-Lecar neurons with conductance-based
coupling to determine under which parameters and graphs synchronous solutions
are stable. We consider connectivity graphs with a constant row-sum, where the
MSF approach can be readily extended to conductance-based synapses rather than
the more well studied diffusive connectivity case, which primarily applies to
gap junction connectivity. In this formulation, the synchronous solution is a
single, self-coupled or ’autaptic’ neuron. We find that the primary
determining parameter for the stability of the synchronous solution is,
unsurprisingly, the reversal potential, as it largely dictates the
excitatory/inhibitory potential of a synaptic connection. However, the change
between an “excitatory” and “inhibitory” synapses is rapid, with only a few
millivolts separating stability and instability of the synchronous state for
most graphs. We also find that for specific coupling strengths (as measured by
the global synaptic conductance), islands of synchronizability in the MSF can
emerge for inhibitory connectivity. We verified the stability of these islands
by direct simulation of pairs of neurons coupled with eigenvalues in the
matching spectrum. These results were robust for different transitions to
spiking (Hodgkin Class I vs Class II), which displayed very similar
synchronizability characteristics.

##  1 Introduction

Brain cells, like other complex interacting systems can readily synchronize
and fire their action potentials or spikes simultaneously, under the right
conditions. Sometimes, this synchronizability is a normal part of brain
function. For example, pyramidal neurons in the hippocampus collectively fire
during the 100-150 millisecond hippocampal sharp-wave ripples, a so-called
cognitive biomarker for memory consolidation and memory replay [1, 2, 3]. This
synchronization is strong enough to be observed even in the hippocampal local
field potential, a macroscopic observable of collective neuronal activity.
However, neurons can also pathologically synchronize [1]. In fact, the same
hippocampal neurons that synchronize in sharp-wave ripples in the hippocampus
often synchronize excessively leading to epileptic seizures [1]. The myriad of
conditions that can lead to synchronization or de-synchronization remains an
active area of research.

One hypothesis is that the specific characteristics of the connectivity
between neurons can promote or obstruct the ability of neurons to otherwise
synchronize. Thus, with a network that is only constrained by a fixed number
of connections, or with a fixed global connection strength, different
connectivity profiles can lead to synchronization or any number of
asynchronous states.

Fortunately, there is a tool to analyze how different networks of neuronal
models can synchronize: the master stability function (MSF) [4, 5, 6, 7, 8].
The MSF allows one to determine the stability of a synchronous solution across
any connectivity graph through a three-step process. First, the master
stability function is computed in the complex plane. Next, the eigenvalues are
computed for any particular connectivity graph. Finally, the sign of the MSF
evaluated at the eigenvalues of the graph on the complex plane dictates the
local asymptotic stability of the synchronous solution. If all eigenvalues
fall in regions where the MSF is negative, the synchronous solution is locally
asymptotically stable while a single eigenvalue falling into a region where
the MSF is positive destabilizes the synchronous solution [4, 5].

However, the application of the MSF in networks of spiking neurons with
conductance-based coupling or chemical synapses has remained under studied for
a few reasons. The first reason is that many neuron models or synaptic models
are non-smooth. For example, all integrate-and-fire neurons utilize a
discontinuity to reset the membrane potential of a neuron after a “spike”, and
possibly change other state variables [9, 10, 11, 12, 13]. While there have
been considerable advances in extending the MSF approach to non-smooth spiking
networks and non-smooth differential equations in general [6, 7, 8, 14], real
neurons do not have the types of membrane discontinuities exhibited by
integrate-and-fire neurons. The second reason is that the connections that
real neurons form are not entirely “excitatory” or “inhibitory” [15]. The net
effect of a chemical synapse depends on the driving force, which is influenced
by ionic reversal potentials and the voltage and spike-shape of a neuron
itself. Thus, a connection can switch between excitatory to inhibitory during
a single spike depending on the voltage of the neuron at any moment. Finally,
the master stability function is primarily limited to cases where the row-sum
of the connectivity matrix is 0, which forces both positive and negative
connection weights [4, 5]. This constraint, which is termed “diffusive
connectivity” is mathematically convenient as it implies that the synchronous
solution in the network is also simultaneously a solution to a single
uncoupled neuron. Unfortunately, this constraint is incompatible with chemical
synapses where the unitary synaptic conductances, being physical quantities,
are non-negative. This is however compatible with gap-junction based
connectivity or electrical synapses. Indeed, this latter case has been
extensively studied [16, 17, 18, 19, 20, 21, 22].

Here, we apply the MSF approach to networks of smooth conductance-based
neurons with smooth conductance-based or chemical synaptic coupling. By
utilizing a constant row-sum, rather than 0 row-sum constraint we successfully
applied the MSF approach to analyzing the synchronizability of networks of
Morris-Lecar neurons with chemical synapses [23, 24]. This constant row-sum
corresponds to a global conductance strength, but not necessarily a measure of
inhibition or excitation strength. We found that independent of the uncoupled
neurons’ bifurcation to spiking (Hodgkin Class I or Hodgkin Class II), the
reversal potential of the synapse would most strongly dictate the stability of
synchronous solutions. We also found that the MSF would rapidly change sign as
a function of the reversal potential over the complex plane. Only a few
millivolts separated large-scale stability or instability of the synchronous
solution. The global portrait of synchronizability looked similar for both
Hodgkin Class I and Class II parameter regimes. However, the MSF deviated from
this result in two ways. First, the synchronous solution could lose stability
depending on the value of the global unitary synaptic conductance and the
reversal potential for excitatory coupling parameter ranges. This loss of
stability was readily observed with sufficiently large ring structures.
Second, we found that for both classes of firing, islands of stability would
emerge for the synchronous solution for inhibitory synapses. We tested this
finding directly with simulations of networks with constrained connectivity
matrices that forced the eigenvalues of the matrices into the
synchronizability islands. The synchronizability of any particular graph
exhibits regimes that are highly parameter dependent in conductance-based
neurons with chemical synapses. Stable configurations can emerge with
increasing inhibition and unstable configurations can emerge with increasing
excitation.

##  2 Results

###  2.1 The Master Stability Function for Conductance Models with
Conductance-Based Coupling

To investigate the synchronization of spiking networks with conductance-based
synapses, we used the general form for a conductance-based neuron with
conductance-based coupling. The network equations are given by:

$\displaystyle C\frac{dV_{i}}{dt}$ $\displaystyle=$ $\displaystyle
F(V,\bm{x}_{i})-\sum_{j=1}^{N}g_{ij}r_{j}(t)(V_{i}-E)$ (1)
$\displaystyle\frac{d\bm{x}_{i}}{dt}$ $\displaystyle=$ $\displaystyle
G(V_{i},\bm{x}_{i})$ (2) $\displaystyle\frac{dr_{i}}{dt}$ $\displaystyle=$
$\displaystyle a_{r}T(V_{i})(1-r_{i})-a_{d}r_{i},$ (3)

where $V$ corresponds to the voltage of the neuron, and $\bm{x}$ consists of a
vector of gating variables. Although our derivation was for a general
conductance-based neuron, we restricted our numerical analysis primarily to
the Morris-Lecar neuron model [23, 24] which we describe below. These neurons
are coupled with a smooth synaptic gating variable equation (3) which was
first introduced in [15]. The synaptic connection between neuron $j$ to neuron
$i$ is given by $\bar{g}_{ij}r_{j}(t)(V_{i}-E)$, where the
inhibitory/excitatory nature of the synapse is determined by the reversal
potential $E$. The function $T(V_{j})$ models the amount of neurotransmitter
released in the synaptic cleft by neuron $j$, and is given by:

$\displaystyle T(V)=\frac{\bar{T}_{max}}{1+\exp\left(-(V-V_{T})/K_{p}\right)}$

Finally, the conductance matrices $\bm{g}$ is an $N\times N$ matrix with the
following constraints:

$\displaystyle g_{ij}\geq 0,\forall i,j=1,2,\ldots
N,\quad\sum_{j=1}^{N}\bar{g}_{ij}=G,\forall i=1,2,\ldots N$ (4)

Thus, all the unitary conductances must be positive. The constraints 4 are
critical for the application of a modified Master Stability Function (MSF)
analysis of the synchronized solutions.

With the constraint (4) in hand, then the synchronous solution corresponds to
the dynamics of a self-coupled or autaptic neuron:

$\displaystyle C\frac{dV_{S}}{dt}$ $\displaystyle=$ $\displaystyle
F(V_{S},\bm{x}_{S})-\bar{g}r_{S}(t)(V_{S}-E_{E})$
$\displaystyle\frac{d\bm{x}_{S}}{dt}$ $\displaystyle=$ $\displaystyle
G(V_{S},\bm{x}_{S})$ $\displaystyle\frac{dr_{S}}{dt}$ $\displaystyle=$
$\displaystyle a_{r}T(V_{S})(1-r_{S})-a_{d}r_{S}$

where $(V_{S}(t),\bm{x}_{S}(t),r_{S}(t))$ is the solution to the synchronous
differential equations (2.1)-(2.1).

Next, we perturbed around the synchronous solution with:

$V_{i}=\epsilon^{V}_{i}+V_{S}(t),\bm{x}_{i}=\epsilon^{\bm{x}}_{i}+\bm{x}_{S}(t),r^{E}_{i}=\epsilon^{r^{E}}_{i}+r_{S}(t)$

which yielded the following linearization

$\displaystyle C\frac{d\bm{\epsilon}^{V}}{dt}$ $\displaystyle=$
$\displaystyle\left(\frac{\partial F}{\partial
V}-\bar{g}r_{S}\right)\bm{\epsilon}^{V}+\sum_{j=1}^{m}\frac{\partial
F}{\partial x_{j}}\bm{\epsilon}^{x_{j}}-(V_{S}-E)\bm{g}\bm{\epsilon}^{r}$ (5)
$\displaystyle\frac{d\bm{\epsilon}^{\bm{x}_{i}}}{dt}$ $\displaystyle=$
$\displaystyle\frac{\partial G}{\partial
V}\bm{\epsilon}^{V}+\sum_{j=1}^{m}\frac{\partial G_{i}}{\partial
x_{j}}\bm{\epsilon}^{x_{j}}$ (6)
$\displaystyle\frac{d\bm{\epsilon}^{r^{E}}}{dt}$ $\displaystyle=$
$\displaystyle
a_{r}T^{\prime}(V_{s})(1-r_{s})\bm{\epsilon}^{V}-(a_{r}T(V_{S})+a_{d})\bm{\epsilon}^{r}$
(7)

Note that all of the partial derivative terms (e.g. $\frac{\partial
F}{\partial V},\frac{\partial G}{\partial x_{j}}$) are evaluated with respect
to the synchronized solution $(V_{S}(t),\bm{x}_{S}(t),r_{S}(t))$.

Next, we will make the standard assumption in MSF-applications that the matrix
$\bm{g}$ is diagonalizable

$\displaystyle\bm{g}=\bm{P}\bm{D}\bm{P}^{-1}$

Then consider the substitution:

$\displaystyle\bm{\eta}^{V}=\bm{P}^{-1}\bm{\epsilon}^{V},\quad\bm{\eta}^{x_{i}}=\bm{P}^{-1}\bm{\epsilon}^{x_{i}},\quad\bm{\eta}^{r}=\bm{P}^{-1}\bm{\epsilon}^{r}.$
(8)

Substituting (8) into equations (5)-(7) yields:

$\displaystyle C\frac{d\bm{\eta}^{V}}{dt}$ $\displaystyle=$
$\displaystyle\left(\frac{\partial F}{\partial
V}-\bar{g}r_{S}\right)\bm{\eta}^{V}+\sum_{j=1}^{m}\frac{\partial F}{\partial
x_{j}}\bm{\eta}^{x_{j}}-(V_{S}-E)\bm{D}\bm{\eta}^{r}$
$\displaystyle\frac{d\bm{\eta}^{\bm{x}_{i}}}{dt}$ $\displaystyle=$
$\displaystyle\frac{\partial G}{\partial
V}\bm{\eta}^{V}+\sum_{j=1}^{m}\frac{\partial G_{i}}{\partial
x_{j}}\bm{\eta}^{x_{j}}$ $\displaystyle\frac{d\bm{\eta}^{r^{E}}}{dt}$
$\displaystyle=$ $\displaystyle
T^{\prime}(V_{S})(1-r_{S})\bm{\eta}^{V}-(a_{r}T(V_{S})+a_{d})\bm{\eta}^{r}.$

The key insight drawn from the MSF approach is that equation (9) is now
effectively uncoupled as the system has been block diagonalized (see Appendix
1 for further details). This implies that to determine the stability of the
synchronized solution, we can determine Lyapunov exponents of the system

$\displaystyle C\frac{d\eta^{V}}{dt}$ $\displaystyle=$
$\displaystyle\left(\frac{\partial F}{\partial
V}-\bar{g}r_{S}\right)\eta^{V}+\sum_{j=1}^{m}\frac{\partial F}{\partial
x_{j}}\eta^{x_{j}}-\lambda_{i}(V_{S}-E)\eta^{r}$ (9)
$\displaystyle\frac{d\eta^{x_{i}}}{dt}$ $\displaystyle=$
$\displaystyle\frac{\partial G}{\partial
V}\eta^{V}+\sum_{j=1}^{m}\frac{\partial G_{i}}{\partial x_{j}}\eta^{x_{j}}$
(10) $\displaystyle\frac{d\eta^{r}}{dt}$ $\displaystyle=$ $\displaystyle
a_{r}T^{\prime}(V_{S})(1-r_{S})\eta^{V}-(a_{r}T(V_{S})+a_{d})\eta^{r}.$ (11)

over a mesh in $\lambda$. The system (9)-(11) is numerically integrated along
with the synchronous solution (2.1)-(2.1), in conjunction with a numerical
estimation of the Lyapunov exponents (see Supplementary Information for
Details). Then, the stability the synchronized solution can be determined
readily by first computing eigenvalues of the connectivity matrix, and then
“looking up” the values on the mesh produced by the MSF .

###  2.2 The Morris-Lecar Neuron Model

To test the predictions of the MSF function, we primarily considered the
Morris-Lecar neuron model [23, 24] :

$\displaystyle C\frac{dV}{dt}$ $\displaystyle=$ $\displaystyle
I-g_{L}(V-E_{L})-g_{K}n(V-E_{K})-g_{Ca}m_{\infty}(V)(V-E_{Ca})$
$\displaystyle\frac{dn}{dt}$ $\displaystyle=$
$\displaystyle\phi\left(\frac{n_{\infty}(V)-n}{\tau_{n}(V)}\right)$
$\displaystyle m_{\infty}(V)$ $\displaystyle=$
$\displaystyle\frac{1}{2}(1+\tanh(V-V_{1})/V_{2})$ $\displaystyle
n_{\infty}(V)$ $\displaystyle=$
$\displaystyle\frac{1}{2}(1+\tanh(V-V_{3})/V_{4})$ $\displaystyle\tau_{n}(V)$
$\displaystyle=$ $\displaystyle\frac{1}{\cosh(V-V_{3})/(2V_{4})}$

The parameters for this model are given in Table 1 and correspond to two
classical regimes, the Hodgkin Class I regime which corresponds to a Saddle
Node on an Invariance Cycle (SNIC) bifurcation from quiescence to spiking, and
the Hodgkin Class II regime which corresponds to a subcritical Hopf
bifurcation followed by a saddle-node of limit cycles from quiescence to
spiking [24]. The network equations are given by

$\displaystyle C\frac{dV_{i}}{dt}$ $\displaystyle=$ $\displaystyle
I-g_{L}(V_{i}-E_{L})-g_{K}n_{i}(V_{i}-E_{K})-g_{Ca}m_{\infty}(V_{i})(V_{i}-E_{Ca})-\sum_{j=1}^{N}g_{ij}r_{j}(t)(V_{i}-E)$
$\displaystyle\frac{dn_{i}}{dt}$ $\displaystyle=$
$\displaystyle\phi\left(\frac{n_{\infty}(V_{i})-n_{i}}{\tau_{n}(V_{i})}\right)$
$\displaystyle\frac{dr_{i}}{dt}$ $\displaystyle=$ $\displaystyle
a_{r}T(V_{i})(1-r_{i})-r_{i}a_{d}$ $\displaystyle\sum_{j=1}^{N}g_{ij}$
$\displaystyle=$ $\displaystyle\bar{g},\forall i,\quad g_{ij}\geq 0$

while the synchronized solution corresponds to the autaptic Morris-Lecar
neuron:

$\displaystyle C\frac{dV_{S}}{dt}$ $\displaystyle=$ $\displaystyle
I-g_{L}(V_{S}-E_{L})-g_{K}n_{S}(V_{S}-E_{K})-g_{Ca}m_{\infty}(V_{S})(V_{S}-E_{Ca})-\bar{g}r_{S}(t)(V_{S}-E)$
(12) $\displaystyle\frac{dn_{S}}{dt}$ $\displaystyle=$
$\displaystyle\phi\left(\frac{n_{\infty}(V_{S})-n_{S}}{\tau_{n}(V_{S})}\right)$
(13) $\displaystyle\frac{dr_{S}}{dt}$ $\displaystyle=$ $\displaystyle
a_{r}T(V_{S})(1-r_{S})-r_{S}a_{d}$ (14)

Prior to computing the MSF, we investigated which regions of $\bar{g}$ lead to
stable spiking solutions as a function of the driving current $I$, and the
reversal of the synapses $E$ for the synchronous solution in equations
(12)-(14). First, we found that excitatory self-coupling ($E=0$ mV) did not
change the overall bifurcation types for Hodgkin Class I or Hodgkin Class II
parameter regimes (Figure 1). However, we did find changes to the overall
bifurcation structure with inhibitory self-coupling would either shift the
spiking regimes to non-physical negative $\bar{g}$, or reduce these regimes to
narrow parameter ranges in $\bar{g}$ depending on how close the driving
current was to the bifurcation point of the non-autaptic neuron. Thus, we
primarily focused on parameter regimes with higher driving currents for both
excitatory ($E=0$ mV) and inhibitory ($E=-70$ mV) reversal potentials.
Furthermore, we found that for these regimes, we could systematically vary
$\bar{g}$ over the interval $(0,3]$ across both Class I and Class II parameter
regimes, albeit with different applied currents ($I=50$ pA for Class I,
$I=115$ pA for Class II). For the currents we considered, the local
bifurcation structure appears largely identical for both classes of firing
(Figure 1), with both parameter regimes stopping spiking via a supercritical
Hopf ($E=-70$ mV) bifurcation for $\bar{g}>3$ or subcritical Hopf ($E=0$ mV)
regimes.

Figure 1: The bifurcation diagrams for self-coupled Morris-Lecar neuron
models. The left column corresponds to the Hodgkin Class I parameter regime,
while the right column corresponds to the Hodgkin Class II parameter regime.
The parameters correspond to excitatory self-coupling near the onset to
spiking (low $I$), excitatory self-coupling far from the onset to spiking
(high $I$), and inhibitory self-coupling far from the onset to spiking (high
$I$). The inhibitory self-coupling regime near the onset to spiking was not
considered as it typically leads to narrower parameter regions in $\bar{g}$,
or can sometimes lead to spiking in non-physical regimes (e.g. $\bar{g}<0$).
The master stability functions were computed for Hodgkin Class I and Hodgkin
Class II parameters in the large $I$ regimes for similar ranges in $\bar{g}$.

###  2.3 Computing and Validating the Master Stability Function for the
Morris-Lecar Model

The block linearization for the MSF function of the Morris-Lecar network is
given by

$\displaystyle C\frac{d\eta^{V}}{dt}$ $\displaystyle=$
$\displaystyle\left(-g_{L}-g_{Ca}m^{\infty}(V_{S})-g_{Ca}(V_{S}-E_{Ca})\frac{dm^{\infty}}{dV}g_{K}n_{S}-\bar{g}r_{S}\right)\eta_{V}$
(15) $\displaystyle-$ $\displaystyle
g_{K}(V_{s}-E_{k})\eta^{n}-\bar{g}\lambda(V_{S}-E)\eta^{r}$
$\displaystyle\frac{d\eta^{n}}{dt}$ $\displaystyle=$
$\displaystyle\phi\left(\frac{1}{\tau_{n}(V_{S})}\frac{dn^{\infty}(V_{S})}{dV}-\frac{d\tau_{n}}{dV_{S}}\frac{n_{\infty}(V_{S})-n_{S}}{\tau_{n}(V_{S})^{2}}\right)\eta_{V}-\frac{\phi}{\tau_{n}(V_{S})}\eta^{n}$
(16) $\displaystyle\frac{d\eta^{r}}{dt}$ $\displaystyle=$ $\displaystyle
a_{r}\frac{dT(V_{S})}{dV}(1-r_{s})\eta^{V}-(a_{r}(T(V_{S}))+a_{d})\eta^{r}$
(17)

where $(V_{S},\eta_{S},r_{S})$ correspond to the synchronized solution for the
self-coupled node:

To compute the master stability function, we simulated the system of equations
(12)-(14) while computing the Lyapunov exponents in parallel for each of the
blocks in (15)-(17). The eigenvalues $\lambda$ were selected over a $101\times
101$ mesh over the unit cube $[a,b]\in[0,1]^{2}$ with $\lambda=a+bi$.

A single computation for the MSF for $\bar{g}=2.1$ nS and $E=30$ mV is shown
in Figure 2. For these particular parameter regimes, the reversal potential
indicates a predominantly excitatory synaptic coupling. This excitatory
connectivity leads to large regions in the eigenvalue space where the MSF is
negative, (Figure 2A-B). For eigenvalues with larger magnitudes, and positive
real components, the MSF exhibits a sign change indicating a loss of stability
in the synchronous solution. To test this loss of stability, we used ring
networks of different sizes as the eigenvalues of a ring of $N$ neurons lie on
the unit circle as the $N$th roots of unity. A ring with 5-neurons has all
eigenvalues laying in the negative MSF area indicating all negative Lyapunov
exponents while a ring with 7 neurons has a pair of complex conjugate
eigenvalues in the positive MSF area. This indicates that the synchronous
solution is stable for a ring of $N=5$ neurons but unstable for $N=7$ neurons.
The ring of $N=6$ neurons has eigenvalues that are extremely close to the sign
change of the MSF, and thus was not considered. The $N=5$ Morris-Lecar ring
and the $N=7$ Morris-Lecar ring were both simulated with initial conditions
near the synchronous solution to test its local asymptotic stability (Figure
2C-D). After the initial transient, the $N=5$ ring converges to the
synchronous solution (Figure 2E) while the $N=7$ ring diverges (Figure 2F).

Figure 2: Testing the master stability function in simulated ring networks (A)
The MSF for the Morris-Lecar neuron under Hodgkin Class I Excitability. The
MSF was computed over a $101\times 101$ evenly distributed mesh over the
complex plane for $\bar{g}=2.1$ nS and $E=30$ mV (excitatory synapses) (B) The
sign of the MSF function. Black denotes negative values, indicating a stable
Lyapunov exponent associated with the eigenvalue, while white denotes positive
Lyapunov exponents associated with the eigenvalue. The MSF was tested with two
rings, a 5-neuron ring and a 7 neuron ring. The 6-neuron ring lay close to the
sign change transition point. (C) The 5-neurons in the ring are simulated near
the basin of attraction of the synchronous solution. The voltage of the
neurons is initialized with a normally distributed random variable with mean
-60 mV, and a standard deviation of 5 mV. The $n$ and $r$ variables are all
set to 0.5. (D) Identical as in (A), only with the 7 neuron ring. (E) The
asymptotic behaviour of the 5-neuron ring is a synchronous solution. (F) The 7
neuron ring desynchronizes after a suitably long period of time.

###  2.4 The Master Stability Function over the ($E,\bar{g}$) parameter space
for Class I and Class II Neurons

Next, we determined how the inhibitory/excitatory valence of the conductance-
based synapse would impact synchronizability. For a sufficiently large
positive reversal potential, the current induced by presynaptic spikes
primarily serves to depolarize the cell, and therefore initiates subsequent
spikes. For a sufficiently large negative reversal potential, the cell becomes
hyperpolarized by presynaptic spikes. The conventional wisdom would be that
increasing the reversal potential would lead to more synchronization as all
the synapses transitioned from inhibitory to excitatory. To test this
hypothesis, we computed the MSF over a discrete mesh in $(\bar{g},E)$ space
for the Morris-Lecar network under both Class I (Figure 3) and Class II
parameters (Figure 4).

We found that the conventional wisdom largely prevails here, with increasing
reversal potentials leading to synchronizability (Figure 3-4). However, there
are some important caveats and unexpected deviations that occur. First, we
noticed that predominantly inhibitory connectivity ($E\leq-50mV$) can lead to
islands of synchronizability when the global conductance strength $\bar{g}$ is
sufficiently high. This is a common feature of MSFs where one can find
localized and compact region(s) of stability [25]. These islands are “born”
and die through a process we describe below. Second, the transition from a
predominantly inhibitory synapse to an excitatory synapse appears quite
suddenly, somewhere in between $E\in[-10,10]$ mV. For the most part, there
does not appear to be a gradual change in synchronizability, but an abrupt
one. Finally, even for “excitatory” synapses, increasing the global
conductance strength can decrease the area where the MSF is negative. This
occurs for example in the $E=10$ mV reversal potential, where a higher
$\bar{g}$ progressively erodes the synchronizability regime. Interestingly,
the loss of stability with increasing $\bar{g}$ can also be reversed for a
sufficiently high reversal potential. Finally, we note that in the weak
coupling regime ($\bar{g}\ll O(1)$), there is a vertical slice of
synchronizability where neurons can synchronize provided that they connect
primarily with strong autaptic connections and weak cross connections to each
other

Figure 3: The sign of the master stability function (MSF) as the network
parameters $\bar{g},E$ are varied The Morris-Lecar parameters were taken to be
in the Hodgkin Class I parameter regime. The sign of the MSF function was
computed over a mesh in the $(\bar{g},E)$ parameter space. The mesh points
correspond to $E=-70+20j$ mV, for $j=0,1,2,\ldots 5$ and$g=0.1+0.5k$ nS for
$k=0,1,2\ldots 7$. In between $E=-30$ and $E=-10$ mV, there is a large-scale
transition from predominantly unstable synchronized solutions and
predominantly stable synchronized solutions.  Figure 4: The sign of the master
stability function (MSF) as the network parameters $\bar{g},E$ are varied. The
Morris-Lecar parameters were taken to be in the Hodgkin Class II parameter
regime. The sign of the MSF function was computed over a mesh in the
$(\bar{g},E)$ parameter space. The mesh points correspond to $E=-70+20j$ mV,
for $j=0,1,2,\ldots 5$ and$g=0.1+0.5k$ nS for $k=0,1,2\ldots 7$. In between
$E=-10$ and $E=10$ mV, there is a large-scale transition from predominantly
unstable synchronized solutions and predominantly stable synchronized
solutions. Note the similarities to Figure 3 Figure 5: The sign of the master
stability function (MSF) as the network parameters $\bar{g},E$ are varied. The
Morris-Lecar parameters were taken to be in the Hodgkin Class II parameter
regime. The sign of the MSF function was computed over a mesh in the
$(\bar{g},E)$ parameter space. The mesh points correspond to $E=-4+j$ mV, for
$j=0,1,2,\ldots 9$ and$g=0.1+0.5k$ nS for $k=0,1,2\ldots 7$. In between $E=-4$
and $E=5$ mV, there is a large-scale transition from predominantly unstable
synchronize solutions and predominantly stable synchronized solutions.

###  2.5 Abrupt Changes to the MSF in the ($E,\bar{g}$) parameter space.

For both Hodgkin Class I and Hodgkin Class II parameter sets, there was a near
global change in the sign of the MSF. For Class I excitability, the MSF
evolved from predominantly positive Lyapunov exponents for $E=-10$ mV, to
predominantly negative Lyapunov exponents for $E=10$ mV. We investigated how
abrupt this transition was by considering the Hodgkin Class I parameters on a
finer mesh. We computed the MSF and plotted its sign from the interval $E=-4$
mV to $E=5$ mV with increments of 1 mV, with $\bar{g}$ fixed as in Figures 3-4
(Figure 5). We found that the MSF was predominantly positive, aside from the
small $\bar{g}$ region. As $E$ is increased however from $E=0$ to $E=5$ mV the
MSF rapidly changes sign over a large area of the admissible eigenvalue space
$|\lambda|\leq 1$. For every millivolt of increase in the reversal potential,
the region of stability for the synchronous solution grows rapidly from the
left side of the admissible eigenvalue domain. For non-weak coupling
$(\bar{g}\gg O(10^{-1})nS$, most admissible weight matrices go from unstable
synchronous solutions to stable synchronous solutions with the change of just
a few millivolt’s in the reversal potential of the synapse.

Next, we tested these abrupt changes in synchronizability by using randomly
coupled networks of Morris-Lecar neurons (Figure 6). The network was simulated
at two parameter values, with $\bar{g}=2.1$ nS and $E=0$ mV, and $\bar{g}=2.1$
nS and $E=4$ mV. The former was predicted to yield unstable synchronized
solutions while the latter had a large region of stability. The network
consisted of 50 neurons, with 85% sparse random coupling, and all gating
variables initialized to 0.5, while the voltage variable was drawn from a
normal random variable with a mean of -60 mV and a standard deviation of 3 mV
(Figure 6A-B). The eigenvalues of the matrix were verified to lie in the
positive sign MSF ($E=0$ mV) and the negative sign MSF ($E=4$ mV),
respectively (Figure 6A-B). The network simulated with the lower reversal
potential desynchronized after 1.5 seconds while the network at the slightly
higher reversal potential synchronized (Figure 6C-D). These results
demonstrate how a change of only a few millivolts in the reversal potential
changes the stability of the synchronous solution for many types of network
coupling.

Figure 6: A change of a few mV in the reversal potential changes the stability
of the synchronous solution. (A) The sign of the MSF function for
$\bar{g}=2.1$, $E=0$ mV. The eigenvalues (red dots) of a randomly generated
$N=50$ node Erdős–Rényi model network that is 85% sparse is generated are also
plotted. (B) The sign of the MSF function for $\bar{g}=2.1$, $E=4$ mV. The
eigenvalues (red dots) are identical as in (A). (C) A simulation of a network
Morris-Lecar neurons with $\bar{g}$ and $E$ identical as in (A), coupled with
the weight matrix from (A). Every neuron is generated with $n_{j}=0.5$,
$r_{j}=0.5$, and a random initial voltage drawn from a normal distribution
with a mean of -60 mV, and standard deviation of 3 mV. (D) An identical
simulation as in (C), only with $\bar{g}=2.1$ and $E=4$ mV. The initial
conditions are identical as in (C). The synchronous solution is now locally
asymptotically stable. The network parameters were in the Hodgkin Class I
regime.

###  2.6 Islands of Synchronizability in the ($E,\bar{g}$) parameter space.

To witness the evolution of an island of synchronizability, we first computed
the MSF function over a finer mesh in the $\bar{g}$ space with $E$ fixed to
$-70$ mV (Figure 7). In these islands, the MSF has a local minimum which can
emerge from the left-hand side of the eigenvalue domain (Figure 7). The island
transitions to the right-hand side of the eigenvalue domain while expanding in
size. Eventually, the island is “absorbed” by the neutral stability eigenvalue
at $\lambda=1$ (Figure 7). Multiple islands (negative minima of the MSF
function) can co-exist (Figure 7).

Next, we investigated if these islands of synchronizability could be directly
tested with a network simulation (Figure 8). As the islands bound small areas
in the complex plane, we tested if networks of pairs $N=2$ of oscillators
would synchronize (Figure 82A, 2D). We constrained the connection strengths of
these matrices to have a constant row sum of 1, which forces the maximum
eigenvalue of $\lambda=1$. It is then trivial to constrain the connections
such that the second eigenvalue lies in an arbitrary position on the real line
(Figure 8 2B, E). For each parameter regime considered, two networks were
tested. One network was selected with an eigenvalue inside the island, and
another with an eigenvalue outside of the island. As predicted by the MSF
function, we found that the connectivity matrices with eigenvalues within the
islands of synchronizability lead to local asymptotic stability of the
synchronous solution. Given the small size of these islands, our results
suggest that for conductance-based synapses, the connectivity graphs in which
inhibition can induce synchronization may be highly constrained and strongly
parameter dependent.

Figure 7: The “birth” and “death” of an island of synchronizability as
$\bar{g}$ is increased. The MSF was computed for the Morris-Lecar neuron model
under Hodgkin Class I excitability. The island of stability emerges from the
negative real part of the eigenvalue mesh (1st row). Then, the island
undergoes a period of expansion and lateral movement to larger real components
(2nd to 4th rows). The island subsequently collides with neutral stability
eigenvalue $\lambda=1$ and begins to be “absorbed” by it (5th row). A second,
very small island of synchronizability has emerged in the last column. The
reversal potential was -70 mV.  Figure 8: Testing the islands of stability
with simulated pairs of coupled oscillators. (A) Two connectivity matrices for
a network of two coupled Morris-Lecar neurons with Hodgkin Class I
excitability. The two matrices yield eigenvalues that lie in the island of
synchronizability (top), or are outside and adjacent to it (bottom). Note that
these matrices satisfy both the unity row-sum constraint, and all positive
elements constraint of a conductance matrix. (B) The computed master stability
function for $\bar{g}=2.25$ nS (left), $E=-70$ mV, and the sign of the MSF
(right). Both matrices have eigenvalues of 1 (neutral stability) and a second
eigenvalue less than 1. The matrix on the top of (A) has a second eigenvalue
within the island of stability (red), while the matrix on the bottom has an
eigenvalue outside of the island of stability (blue). (C) Simulation of the
Morris-Lecar neurons coupled as in (A), with only the voltage plotted. The
matrix with an eigenvalue inside of the island has a locally asymptotically
stable synchronous solution. The matrix on the bottom of (A) has an eigenvalue
outside of the island, and has an unstable synchronous solution. (D) Identical
as in (A), only with a larger conductance value $\bar{g}=2.75$ nAS. (E)
Identical as in (B), with the connectivity matrices determined by (D) (F)
Identical as in (C), with the connectivity matrices determined by (D).

##  3 Discussion

We have found that the synchronizability of conductance-based neurons is a
complex affair. By applying a modified MSF approach, we analyzed how networks
of chemically coupled Morris-Lecar neurons synchronize in Class I and Class II
regimes. We found remarkable consistency in the stability regimes across these
two parameter regimes. As a general rule of thumb, higher reversal potentials
tend to lead to the (local) stability of synchronous autaptic solutions, and
lower reversal potentials tend to lead to the (local) instability of
synchronous autaptic solutions. However, this rule of thumb is often deviated
from as islands of stability under inhibition, and wedge-shaped regions of
instability under excitation were both observed. The actual change between
predominantly stable or predominantly unstable synchronized solutions occurs
very rapidly as a function of the reversal potential (a few mVs). Evidently,
the stability of synchronous solutions is strongly parameter dependent with
small changes to the global conductance or small changes to the reversal
potential stabilizing or destabilizing the synchronizability of different
connectivity graphs.

We note that we are not the only authors to consider applying the MSF function
to neurons with chemical synapses [26, 27]. In [26] the authors consider
networks of Hindmarsh-Rose neuron models with both electrical and chemical
coupling, thereby making a direct comparison to the work here difficult. In
[27], the authors consider networks of Izhikevich neuron models with chemical
synapses and utilize non-smooth analysis to determine the MSF. This was
performed for electrical only, chemical only, and simultaneous
electrical/chemical synapses. For chemical synapses, the estimated MSF appears
positive (Figure 5 in [27]). The modification to the MSF to allow for non-
diffusive coupling appears in multiple sources in the literature [27, 28, 29,
30, 31].

The islands of stability observed here occurred for primarily inhibitory (low
$E$) connection strengths. Interestingly, work with classical integrate-and-
fire neurons with current based synapses has also demonstrated that inhibition
can induce synchronization [32]. In particular, the authors consider leaky
integrate-and-fire neurons with alpha-function like synaptic connectivity,
where every spike at time $t^{*}$ increases the current arriving to a neuron
by $E_{s}(t)=E_{s}(t-t^{*})$ where $E_{s}=g\alpha^{2}t\exp(-\alpha t)$. The
authors find that for inhibitory synapses ($g<0$), the synchronous solution is
always stable, although at a critical value of $\alpha$ (faster synapses), a
pitchfork bifurcation occurs that stabilizes for the asynchronous solution. As
the synapses become faster, the basin of attraction for the synchronous
solution shrinks. We primarily considered synapses with a rise time of $~{}1$
($a_{r}=1.1$) ms and a decay time of $~{}5$ ms ($a_{d}=0.19$). For
conductance-based synapses, it appears that the synchronous solution is stable
under inhibition dominated regimes only for very specific connectivity graphs,
in contrast with the behaviour of an inhibitory coupled integrate-and-fire
network. We leave the general analysis of the impacts of synaptic timing for
future work, however we do remark that the findings in [32] were supported and
extended by subsequent modelling work [33] and even experimental findings [34,
35].

By thoroughly exploring parameter space, we were able to shed insight into the
impacts of having conductance-based synaptic coupling on synchronizability.
Chemical synapses are not strictly excitatory or inhibitory as the flow of
current is dictated by the driving force, which leads to complex behaviours
and synchronizability regimes, where excitation can lead to desynchronization
or inhibition can lead to synchronization.

### Acknowledgement

WN is supported by an NSERC Discovery Grant (DGECR/00334-2020) and a Canada
Research Chair (CRC-2019-00416)

## Appendix A

### The Master Stability Function

The Master Stability Function (MSF) approach to resolving the stability of
synchronous solutions was originally proposed in [4], with subsequent advances
to the analysis of clustered systems in [5]. Briefly, the approach considers
coupled non-linear systems of the following form:

$\displaystyle\frac{d\bm{x}_{i}}{dt}=Q(\bm{x}_{i})+\sum_{j=1}^{N}g_{ij}H(\bm{x}_{j}),\quad
i=1,2,\ldots N,\quad\bm{x}_{i}\in\mathbb{R}^{k},$

with the constraint that

$\displaystyle\sum_{j=1}^{N}g_{ij}=0,\forall i$

and that the coupling matrix $\bm{g}$ be diagonalizable. This is sometimes
referred to as “diffusive connectivity”

One can express the system in (The Master Stability Function) with the direct
product $\otimes$ or tensor operator as

$\displaystyle\frac{d\bm{x}}{dt}=\bm{Q}(\bm{x})+\bm{G}\otimes\bm{H}(\bm{x})$

where

$\bm{Q}(\bm{x})=[Q(\bm{x}_{1}),Q(\bm{x}_{2}),\ldots
Q(\bm{x}_{N})],\quad\bm{H}(\bm{x})=[H(\bm{x}_{1}),H(\bm{x}_{2})\ldots
H(\bm{x}_{N})].$

Under these constraints, the synchronous solution $\bm{x}_{s}(t)$ is a
solution to the differential equation

$\displaystyle\frac{d\bm{x}_{s}}{dt}=Q(\bm{x}_{s})$ (18)

Then, perturbations off of the synchronous solution are considered

$\displaystyle\bm{x}_{i}(t)=\bm{x}_{s}(t)+\bm{\epsilon}_{i}(t).$

which yields the following block diagonal system

$\displaystyle\frac{d\bm{\epsilon}}{dt}=\left[\bm{I}_{N}\otimes\frac{\partial
Q}{\partial\bm{x}}+\bm{G}\otimes\frac{\partial
H}{\partial\bm{x}}\right]\bm{\epsilon}$ (19)

The coupled non-autonomous system in (19) determines the linear stability of
the system (18) via computation of the Lyapunov exponents of the system. In a
master-stability function approach, the problem of computing these Lyapunov is
simplified greatly by assuming that $\bm{G}$ is diagonalizable:

$\displaystyle\bm{G}=\bm{P}\bm{D}\bm{P}^{-1}$

By applying the substitution $\bm{\eta}=\bm{P}^{-1}\bm{\eta}$, the system
simplifies into a diagonalized block-system:

$\displaystyle\frac{d\bm{\eta}_{i}}{dt}=\left[\frac{\partial
Q}{\partial\bm{x}}+\lambda_{i}\frac{\partial
H}{\partial\bm{x}}\right]\bm{\eta}_{i},\quad i=1,2,\ldots N$ (20)

where $\lambda_{i}$ is an eigenvalue of the matrix $\bm{G}$. The next step in
a MSF function approach is to compute the Lyapunov exponents of equation 20
over a mesh in the complex eigenvalue space of $\lambda=a+bi$. Thus, one
computes the mesh of lyapunov exponents, $\mu$ as a function of the
eigenvalues $\lambda$:

$\displaystyle\frac{d\bm{\eta}}{dt}=\left[\frac{\partial
Q}{\partial\bm{x}}+\lambda\frac{\partial
H}{\partial\bm{x}}\right]\bm{\eta},\rightarrow\mu(\lambda)$ (21)

The Lyapunov exponents $\mu(\lambda)$ can be computed with established methods
for their approximation (e.g. [37]).

The maximum Lyapunov exponent of the block 20 as function of $\lambda$ is the
Master Stability Function (MSF). With the mesh computed, one can use any
diagonalizable connection matrix and simply ”look up” the value of the maximum
Lyapunov exponent as a function of $\lambda$ with $\mu(\lambda)$, the MSF. In
this work, we refer to the final diagonalized block structure in (21) as the
MSF equations, as they are necessary for the numerical approximation of
$\mu(\lambda)$.

### Supplementary Material

The parameters for the Morris-Lecar neuron model under Hodgkin Class I and
Hodgkin Class 2 regimes are shown in Table 1, along with the synaptic
parameters.

Parameter | Value | Units
---|---|---
$C$ | 20 | pF
$g_{L}$ | 2 | nS
$E_{L}$ | 60 | mV
$g_{K}$ | 8 | nS
$E_{K}$ | -84 | mV
$g_{Ca}$ | 4 | nS
$E_{Ca}$ | 120 | mV
$V_{1}$ | -1.2 | mV
$V_{2}$ | 18 | mV
$V_{3}$ | 12 (Class I), 2 (Class II) | mV
$V_{4}$ | 17.4 (Class I), 30 (Class II) | mV
$\phi$ | 0.067 (Class I), 0.04 (Class II) | mV
$I$ | 50 (Class I), 115 (Class II) | pA
$a_{r}$ | 1.1 | ms-1
$a_{d}$ | 0.19 | ms-1
$T_{max}$ | 1 | unitless
$K_{p}$ | 5 | mV
$V_{t}$ | 2 | mV
Table 1: Parameter values used for the Morris-Lecar neuron model, unless
otherwise specified in a figure

### Numerical Integration and Computation of the Lyapunov Exponents

All ODEs for the Morris-Lecar system(s) under consideration were integrated
with the MATLAB2023a function _ode45_. The ’RelTol’ (relative error tolerance)
and ’AbsTol’ absolute error tolerance parametes were set to $10^{-14}$ for all
direct integration of the Morris-Lecar network equations. The Lypaunov
exponents were computed with the algorithm in [37] with code modified from
[30, 31].

## References

  * [1] Buzsáki, G.  Hippocampal sharp wave-ripple: A cognitive biomarker for episodic memory and planning.  _Hippocampus_ 25, 1073–1188 (2015).
  * [2] Lee, A. K. & Wilson, M. A.  Memory of sequential experience in the hippocampus during slow wave sleep.  _Neuron_ 36, 1183–1194 (2002).
  * [3] Wilson, M. A. & McNaughton, B. L.  Reactivation of hippocampal ensemble memories during sleep.  _Science_ 265, 676–679 (1994).
  * [4] Pecora, L. M. & Carroll, T. L.  Master stability functions for synchronized coupled systems.  _Physical review letters_ 80, 2109 (1998).
  * [5] Pecora, L. M., Sorrentino, F., Hagerstrom, A. M., Murphy, T. E. & Roy, R.  Cluster synchronization and isolated desynchronization in complex networks with symmetries.  _Nature communications_ 5, 4079 (2014).
  * [6] Coombes, S. & Thul, R.  Synchrony in networks of coupled non-smooth dynamical systems: Extending the master stability function .
  * [7] Lai, Y. M., Thul, R. & Coombes, S.  Analysis of networks where discontinuities and nonsmooth dynamics collide: understanding synchrony.  _The European Physical Journal Special Topics_ 227, 1251–1265 (2018).
  * [8] Coombes, S. & Wedgwood, K. C.  _Neurodynamics: An Applied Mathematics Perspective_ , vol. 75 (Springer Nature, 2023).
  * [9] Izhikevich, E. M.  Simple model of spiking neurons.  _IEEE Transactions on neural networks_ 14, 1569–1572 (2003).
  * [10] Izhikevich, E. M.  _Dynamical systems in neuroscience_ (MIT press, 2007).
  * [11] Touboul, J.  Bifurcation analysis of a general class of nonlinear integrate-and-fire neurons.  _SIAM Journal on Applied Mathematics_ 68, 1045–1079 (2008).
  * [12] Brette, R. & Gerstner, W.  Adaptive exponential integrate-and-fire model as an effective description of neuronal activity.  _Journal of neurophysiology_ 94, 3637–3642 (2005).
  * [13] Naud, R., Marcille, N., Clopath, C. & Gerstner, W.  Firing patterns in the adaptive exponential integrate-and-fire model.  _Biological cybernetics_ 99, 335–347 (2008).
  * [14] Ashwin, P., Coombes, S. & Nicks, R.  Mathematical frameworks for oscillatory network dynamics in neuroscience.  _The Journal of Mathematical Neuroscience_ 6, 1–92 (2016).
  * [15] Destexhe, A., Mainen, Z. F., Sejnowski, T. J. _et al._ Kinetic models of synaptic transmission.  _Methods in neuronal modeling_ 2, 1–25 (1998).
  * [16] Lehnert, J., Dahms, T., Hövel, P. & Schöll, E.  Loss of synchronization in complex neuronal networks with delay.  _Europhysics Letters_ 96, 60013 (2011).
  * [17] Bonacini, E. _et al._ How single node dynamics enhances synchronization in neural networks with electrical coupling.  _Chaos, Solitons & Fractals_ 85, 32–43 (2016).
  * [18] Jalili, M.  Enhancing synchronizability of diffusively coupled dynamical networks: A survey.  _IEEE transactions on neural networks and learning systems_ 24, 1009–1022 (2013).
  * [19] Belykh, I., De Lange, E. & Hasler, M.  Synchronization of bursting neurons: What matters in the network topology.  _Physical review letters_ 94, 188101 (2005).
  * [20] Neefs, P., Steur, E. & Nijmeijer, H.  Network complexity and synchronous behavior—an experimental approach.  _International Journal of Neural Systems_ 20, 233–247 (2010).
  * [21] Keplinger, K. & Wackerbauer, R.  Transient spatiotemporal chaos in the morris-lecar neuronal ring network.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 24 (2014).
  * [22] Rossoni, E., Chen, Y., Ding, M. & Feng, J.  Stability of synchronous oscillations in a system of hodgkin-huxley neurons with delayed diffusive and pulsed coupling.  _Physical Review E_ 71, 061904 (2005).
  * [23] Morris, C. & Lecar, H.  Voltage oscillations in the barnacle giant muscle fiber.  _Biophysical journal_ 35, 193–213 (1981).
  * [24] Ermentrout, B. & Terman, D. H.  _Mathematical foundations of neuroscience_ , vol. 35 (Springer, 2010).
  * [25] Huang, L., Chen, Q., Lai, Y.-C. & Pecora, L. M.  Generic behavior of master-stability functions in coupled nonlinear dynamical systems.  _Physical Review E_ 80, 036204 (2009).
  * [26] Jhou, F.-J., Juang, J. & Liang, Y.-H.  Multistate and multistage synchronization of hindmarsh-rose neurons with excitatory chemical and electrical synapses.  _IEEE Transactions on Circuits and Systems I: Regular Papers_ 59, 1335–1347 (2012).
  * [27] Aristides, R. P. & Cerdeira, H. A.  Master stability functions of networks of izhikevich neurons.  _arXiv preprint arXiv:2303.13921_ (2023).
  * [28] Checco, P., Righero, M., Biey, M. & Kocarev, L.  Synchronization in networks of hindmarsh–rose neurons.  _IEEE Transactions on Circuits and Systems II: Express Briefs_ 55, 1274–1278 (2008).
  * [29] Nicola, W., Hellyer, P. J., Campbell, S. A. & Clopath, C.  Chaos in homeostatically regulated neural systems.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 28 (2018).
  * [30] Nicola, W. & Campbell, S. A.  Normalized connectomes show increased synchronizability with age through their second largest eigenvalue.  _SIAM Journal on Applied Dynamical Systems_ 20, 1158–1176 (2021).
  * [31] Al-Darabsah, I., Chen, L., Nicola, W. & Campbell, S. A.  The impact of small time delays on the onset of oscillations and synchrony in brain networks.  _Frontiers in Systems Neuroscience_ 15, 688517 (2021).
  * [32] Van Vreeswijk, C., Abbott, L. & Bard Ermentrout, G.  When inhibition not excitation synchronizes neural firing.  _Journal of computational neuroscience_ 1, 313–321 (1994).
  * [33] White, J. A., Chow, C. C., Rit, J., Soto-Treviño, C. & Kopell, N.  Synchronization and oscillatory dynamics in heterogeneous, mutually inhibited neurons.  _Journal of computational neuroscience_ 5, 5–16 (1998).
  * [34] Elson, R. C., Selverston, A. I., Abarbanel, H. D. & Rabinovich, M. I.  Inhibitory synchronization of bursting in biological neurons: dependence on synaptic time constant.  _Journal of Neurophysiology_ 88, 1166–1176 (2002).
  * [35] Khazipov, R.  Gabaergic synchronization in epilepsy.  _Cold Spring Harbor perspectives in medicine_ 6 (2016).
  * [36] McDougal, R. A. _et al._ Twenty years of modeldb and beyond: building essential modeling tools for the future of neuroscience.  _Journal of computational neuroscience_ 42, 1–10 (2017).
  * [37] Wolf, A., Swift, J. B., Swinney, H. L. & Vastano, J. A.  Determining lyapunov exponents from a time series.  _Physica D: nonlinear phenomena_ 16, 285–317 (1985).