# A multiple timescale network model of intracellular calcium concentrations in coupled neurons: Insights from ROM simulations

Alejandro Bandera Moreno, Soledad Fernández-García, María Macarena Gómez-Mármol, Alexandre Vidal

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Alejandro Bandera Moreno, Soledad Fernández-García, María Macarena Gómez-Mármol, Alexandre Vidal. A multiple timescale network model of intracellular calcium concentrations in coupled neurons: Insights from ROM simulations. Mathematical Modelling of Natural Phenomena, 2022, 17, pp.11. 10.1051/mmnp/2022016. hal-04462433

HAL Id: hal-04462433

https://hal.science/hal-04462433v1

Submitted on 16 Feb 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# A MULTIPLE TIMESCALE NETWORK MODEL OF INTRACELLULAR CALCIUM CONCENTRATIONS IN COUPLED NEURONS: INSIGHTS FROM ROM SIMULATIONS

A. BANDERA<fnref n="1" />,<fnref n="2" />,* , S. FERNÁNDEZ-GARCÍA<fnref n="1" />,<fnref n="2" /> , M. GÓMEZ-MÁRMOL<fnref n="1" />
AND A. VIDAL<fnref n="3" />

**Abstract.** In Fernández-García and Vidal [*Physica D* **401** (2020) 132129], the authors have analyzed the synchronization features between two identical 3D slow-fast oscillators, symmetrically coupled, and built as an extension of the FitzHugh–Nagumo dynamics generating Mixed-Mode Oscillations. The third variable in each oscillator aims at representing the time-varying intracellular calcium concentration in neurons. The global model is therefore six-dimensional with two fast variables and four slow variables with strong symmetry properties. In the present article, we consider an extension of this model in two different directions. First, we consider heterogeneity among cells and analyze the coupling of two oscillators with different values for one parameter which tunes the intrinsic frequency of the output. We therefore identify new patterns of antiphasic synchronization, with non trivial signatures and that exhibit a Devil’s Staircase phenomenon in signature transitions when varying the coupling gain parameter value. Second, we introduce a network of *N* cells divided into two clusters: the coupling between neurons in each cluster is excitatory, while the coupling between the two clusters is inhibitory. Such system aims at modelling the interactions between neurons tending to synchronization in each of two subpopulations that inhibit each other, like ipsi- and contra-lateral motoneurons assemblies. To perform the numerical simulations in this case when *N* is large, as an initial step towards the network analysis, we consider Reduced Order Models in order to save computational costs. We present the numerical reduction results in a network of 100 cells. For the sake of validation of the numerical reduction method, we both compare the outputs and CPU times obtained with the original and the reduced models in different cases of network coupling structures.

**Mathematics Subject Classification.** 34C15, 34C25, 34C27, 34C28, 92B20, 34K28, 37M05, 70K70.

Received November 11, 2021. Accepted April 11, 2022.

<footnote>*Keywords and phrases:* Slow-fast dynamics, coupled oscillators, mixed-mode oscillations, synchronization, network neuron model, reduced order models.</footnote>

<footnote n="1">¹ Ecuaciones Diferenciales y Análisis Numérico, Universidad de Sevilla, Calle Tarfia s/n, Seville 41012, Spain.</footnote>

<footnote n="2">² IMUS, Universidad de Sevilla, Calle Tarfia s/n, Seville 41012, Spain.</footnote>

<footnote n="3">³ Laboratoire de Mathématiques et Modélisation d'Évry (LAMME), Univ Evry, CNRS, Université Paris-Saclay, IBGBI, 23 Bld de France, Evry 91037, France.</footnote>

<footnote>* Corresponding author: abandera@us.es</footnote>
# 1. INTRODUCTION

Synchronization features between neuron activities play a key role in the treatment of information by the nervous system. The rich panel of synchronization types in neuron network results both from the complexity of the neural ionic dynamics and from the network structure and scale, from a few neuronal cells (microscopic scale), through neuron populations (mesoscopic scale), to large areas (macroscopic scale). The complexity of the ionic dynamics is often reflected in the models of neural activity by the presence of different timescales for the different variables and by the non linearity of the dynamics, as the initial approach of Hodgkin-Huxley [20] and the numerous neural activity models designed afterwards for different scales. Yet, the synchronization features between coupled systems with multiple timescales – relaxation oscillators, bursters, Mixed-Mode Oscillations (MMOs) – differs strongly from those between coupled harmonic oscillators [5, 8, 21, 28, 29, 36, 37, 40] and the role of the coupling strength differs between the two cases. In particular, canard phenomena arising in multiple timescale systems plays a prominent role in organizing the synchronization of weakly coupled slow-fast systems. The studies linked with this problematic focus on different questions, for instance synchronization/desynchronization [10, 12–14, 33], complete synchronization [1, 3, 4, 32, 35], local oscillations and clustering [12, 33, 45]. Yet, the problem of synchronization of coupled multiple timescale systems, especially when MMOs are involved, remains mainly open.

In many context, synchronization of neural activities can be experimentally observed using fluorescence techniques to obtain recordings of Intracellular Calcium Concentrations (ICC) in single neurons connected as a population. In [24], a compact 3D model, extending the classical FitzHugh-Nagumo dynamics [17, 30], has been built for reproducing the ICC variations along time in a single neuron and the method for tuning the parameters to reproduce given intrinsic frequencies, amplitudes and quiescence phase durations in the generated ICC time series. A network model has also been built for reproducing episodic synchronization identified in GnRH neuron populations. Afterwards, authors in [16] have proposed a 6D system obtained by symmetrically coupling two such 3D oscillators generating MMOs and performed a theoretical study of the synchronization features in oscillatory ICC patterns of two coupled cells. The coupling type chosen for building this model relies on the underlying dynamics of ICC, on which the synaptic transmission as well as gap-junctions between neurons has little effect (see [22] for instance). Increase of ICC mainly results from calcium-induced calcium-release (CICR) mechanism from the endoplasmic reticulum of the neural cell. Glial cells such as astrocytes are involved in the CICR channel through the increase of Inositol trisphosphate (IP3). Since this mechanism is slower than neuronal transmission, the coupling in [16] was chosen on the recovery variables that evolves on a consistent timescale for modeling this mechanism [42]. The choice of a linear coupling, although it consists in a strong simplification of the biological interactions, allows us to obtain a mathematically-tractable model while preserving the qualitative properties of the dynamical interactions between neuron ICC levels. Considering both excitatory and inhibitory coupling, the aim of the study was two-fold : (i) identify the conditions for obtaining in-phase and perfectly in-phase synchronization between two coupled cells (for the excitatory case), (ii) analyze the synchronization behaviors of two identical populations of neurons inhibiting each others, considering that one 3D oscillator is also a reduced model for a population of perfectly synchronized neurons, which occurs (following (i)) if the excitatory coupling is strong enough. The analysis performed in [16] is based on geometric singular perturbation theory and takes advantage of the symmetry of the system. We summarize the identified behaviors and some of the results obtained in this case in Section 2, since they are mandatory for understanding the present study.

Although the results obtained in [16] represent a first step in the theoretical understanding of the panel and repartition of synchronization types according to the coupling strength, two assumptions are too restrictive for dealing with realistic problems: first, the two cells are assumed to be identical while heterogeneity in intrinsic quantitative features of single cell ICC oscillations is a well-known fact [38, 44]; second, if we consider each of the 3D systems as a cluster, as mentioned above, we are assuming that the single cells dynamics constituting the cluster are perfectly synchronized. In order to advance in this line towards more realistic problems, in the present article, we consider an extension of the system analyzed in [16] in two different directions:
1. We look upon a heterogeneous case of two cells coupling, *i.e.* two oscillators with different parameter values. We study the synchronization patterns of the system with respect to two parameters: the ratio between the intrinsic frequencies of the oscillators and the coupling strength. Therefore, the system is no longer symmetric and several arguments used for the theoretical analysis in [16] are no longer valid. Thus, we proceed through direct numerical simulations through a fourth order Runge-Kutta solver.

2. We consider a network of $N > 2$ cells divided into two clusters. Then, we have to cope with a stiff system of $3N$ equations. In order to solve the system, in principle, we can use the same solver as in the two-cells system. However, although this technique approximates well the solution, in some cases it is computationally excessive to solve. Therefore, as a first step towards the network analysis, we adapt and apply a Reduced Order Method (ROM) approximation [9, 26, 27, 31, 41]. Such method consists in defining an equivalent dynamical system of lower dimension that approximates the full order model, namely, the reduced problem.

The choice of the connectivity matrix structure in the network system is strongly motivated by its potential application to a precise neuronal system. It has been emphasized in several experiments using fluorescence techniques for recording ICC time series that populations of motoneurons in different species are organized as clusters formed by neurons tending to synchronize their activities (ICC pulses) in-phase or almost in-phase in each cluster, while the clusters seems to mutually inhibit the other ones. In simple animal models, such as the zebrafish, for which several studies have been made, two clusters are identified corresponding to the ipsi- and contra-lateral sides of the spinal cord respectively, with strong correlations between neuron activities in each of them. Yet, the synchronization features, the underlying dynamical mechanisms, as well as the transient situation where both clusters can activate at the same time are misunderstood. The theoretical and numerical analysis of a clustered network model of ICC dynamics mimicking such connectivity will be helpful in understanding the dynamical interaction of motoneurons.

The article outline is the following. In Section 2, we first summarize the essential properties of the 3D model of ICC changes in a single neuron introduced in [24]. We recall the synchronization patterns of the 6D coupled homogeneous system identified in [16] and we describe the new features associated with the heterogeneous case. Then, we introduce the network model ($N > 2$) and explain the choice of the coupling type and network structure. Section 3 focus on the method used to define ROMs. In Section 4 we provide numerical evidences of the accuracy of the reduced model obtained with ROM by showing its ability to approximate closely the time series generated by the original system for both the main behaviors and the most singular ones. We also show the reduction of CPU time needed for the numerical integration of the reduced system compared to the original one. Finally, Section 5 is devoted to present conclusions and perspectives of the present work. The more technical questions of this work have been relegated to the Appendix A.

## 2. MODELS OF NEURONAL ICC DYNAMICS

We recall quickly the dynamics of the 3D one-cell model and the formerly identified behaviors of the 6D system in the homogeneous symmetric coupling case. Then, we focus on the heterogeneous case and present the new specific synchronization features according to the ratio between the values of the parameter that are different in each of the two cells. We show the structure of the changes in frequencies while the value of the coupling gain parameter varies. Finally, Section 2.3 is devoted to set the clustered network model involving heterogeneity between cells as well as in the coupling gain value between cells.

### 2.1. Single-cell model of neuronal ICC dynamics

The single-cell system used in [24] is a 3D slow-fast oscillator designed for modeling the ICC changes in a single neuron. This simple model has been proved, by comparison with experimental data, to be efficient in capturing the essential qualitative and quantitative properties of ICC oscillations in GnRH neurons as well as motoneurons. Moreover, in [24], the authors have shown how to tune the parameters in a simple way for reproducing given frequencies, amplitudes and quiescence phase durations in ICC oscillations.
The system reads

$$ 
\left\{ 
\begin{aligned}
\dot{x} &= \tau(-y + f(x) - \phi_f(z)), \\
\dot{y} &= \tau\varepsilon(x + a_1y + a_2), \\
\dot{z} &= \tau\varepsilon\left(\phi_r(x) - \frac{z - z_b}{\tau_z}\right),
\end{aligned}
\right. 
\qquad (2.1) $$

where,

$$ \phi_f(z) = \frac{\mu z}{z + z_0}, \qquad (2.2) $$

$$ \phi_r(x) = \frac{\lambda}{1 + \exp(-\rho(x - x_{on}))}, \qquad (2.3) $$

with $(x, y, z) \in \mathbb{R}^3$. The overdot stands for the derivative with respect to $t$ (the time variable). Parameter $\varepsilon$ fits the timescale separation and therefore $0 < \varepsilon \ll 1$. As in the FitzHugh-Nagumo system (see [17, 30]) and 3D extensions of such dynamics (see for instance [23] for a model of dopaminergic neurons), fast variable $x$ stands for the cell electrical activity, while slow variable $y$ acts as a recovery variable. Slow variable $z$ stands for the ICC, and its dynamics is essentially driven by sigmoid $\phi_r$. Variable $z$ feedbacks onto the electrical activity through the Hill function $\phi_f$ with superior bound $\mu$. For $\phi_r(x)$ close to 0, $z$ decreases down to a value close to $z_b$ that is a quasi-steady state for $z$ in case of electric inactivation and stands for ICC baseline. The decrease of variable $z$ is an exponential decay with $\tau\varepsilon/\tau_z$ rate. The coupling term $\phi_f(z)$ onto $x$ aims at reproducing the reduction neuronal electrical activity when ICC increases. Such coupling has already been introduced in single neuron models to mimick the neuron hyperpolarization induced by calcium (see [11] for instance). The essential properties of ICC dynamics depend on the qualitative features of the coupling functions. These properties and the results of our study are preserved for other increasing functions $\phi_f$ and $\phi_r$ that are bounded on $\mathbb{R}_+$. For reproducing the excitable feature of calcium release from the endoplasmic reticulum, function $\phi_r$ is chosen stiff and a sigmoid is well-designed for this aim.

In this article, for parameters that are fixed among the simulations, we have kept the values used in [16]. Therefore, we consider $f(x) = -x^3 + 4x$ and

$$ 
\begin{array}{l}
a_1 = -0.1, \quad a_2 = 0.8, \quad \varepsilon = 0.06, \quad \tau = 37, \quad \mu = 2.4, \quad z_0 = 5, \\
z_b = 1, \quad \tau_z = 2, \quad x_{on} = -0.45, \quad \lambda = 1.75, \quad \rho = 4.5.
\end{array} 
\qquad (2.4) $$

Note that parameter $\tau$ was originally introduced in [24] for fitting the timescale to the real time and mimick experimental data. This parameter does not impact the qualitative behaviors presented in this study, yet we keep it for the sake of continuity with former studies.

For describing the main system behaviors, we define the S-shaped critical surface

$$ S \equiv \{y = f(x) - \phi_f(z)\}. $$

Its two fold lines are

$$ \mathcal{F}^{-} \equiv \{x = -x_{f}, y = f(-x_{f}) - \phi_{f}(z)\}, $$

and

$$ \mathcal{F}^{+} \equiv \{x = x_{f}, y = f(x_{f}) - \phi_{f}(z)\}. $$
<figure><img src="image_3.png" /><figcaption>(a) Attractive relaxation limit cycle of the single-cell ICC model.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(c) MMOs generated by the single-cell ICC model. Left panel: Periodic orbit. Right panel: magnified view of the small oscillations near the lower fold of the critical manifold.</figcaption></figure>

<caption>FIGURE 1. Asymptotic behaviors of system (2.1) depending on the position of its singular point together with the S-shaped critical manifold. Single and double arrows materialize the slow and fast parts of the limit cycle, respectively.</caption>

The critical surface consists of three sheets and the fold lines that separate them. The left sheet $S_l$ contained in $x < -x_f$ as well as the right sheet $S_r$ contained in $x > x_f$ are attracting manifolds for the boundary-layer system (obtained for the limit case $\varepsilon = 0$), i.e. the fast $x$-dynamics with $y, z$ considered as parameters. The middle sheet $S_m$, contained in $|x| < x_f$, is repulsive for the fast dynamics.

According to the position of system (2.1) singular point, three main asymptotic behaviors can be identified using slow-fast analysis, i.e. considering $\varepsilon > 0$ small enough and using geometric singular perturbation theory. Rigorously, the phase portraits described below are obtained for any $\varepsilon \in (0, \varepsilon_0]$ for a fixed $\varepsilon_0 > 0$. In particular, in [16] it is proved that for $\mu = 2.51$, system (2.1) admits a folded saddle-node type II singularity in the singular limit (see [25]). This implies that a singular Hopf bifurcation takes place within a $O(\varepsilon)$-distance in the parameter space. We refer the reader to [16, 24] for more detailed explanations of the following phase portraits but display Figure 1a and b for visual support.
1. The singular point lies on $S_m$ (the middle sheet) and far from the folds: then, it is unstable and there exists a globally attractive limit cycle of system (2.1) generating relaxation oscillations in the time series (Fig. 1a).

2. The singular point lies on $S_m$ but close to a fold: then, it is unstable and system (2.1) generates Mixed-Mode Oscillations (Fig. 1c).

3. The singular point lies on $S_l$: then, it is asymptotically stable and any orbit of system (2.1) admits this point as $\omega$-limit (Fig. 1b). The phase portrait has the same properties if the singular point lies on $S_r$.

## 2.2. 6D model of ICC dynamics in two coupled cells

The two-cells model is built as a linear coupling of two copies of system (2.1). As explained in the Introduction, the timescale of CICR mechanisms is consistent with a coupling involving the fast variables onto the slow recovery variables dynamics between the two oscillators:

$$O_1 \left\{ \begin{aligned} \dot{x}_1 &= \tau(-y_1 + f(x_1) - \phi_f(z_1)), \\ \dot{y}_1 &= \tau\varepsilon k_1(x_1 + a_1y_1 + a_2 + c(x_1 - x_2)), \\ \dot{z}_1 &= \tau\varepsilon \left( \phi_r(x_1) - \frac{z_1 - z_b}{\tau_z} \right), \\ \dot{x}_2 &= \tau(-y_2 + f(x_2) - \phi_f(z_2)), \\ \dot{y}_2 &= \tau\varepsilon k_2(x_2 + a_1y_2 + a_2 + c(x_2 - x_1)), \\ \dot{z}_2 &= \tau\varepsilon \left( \phi_r(x_2) - \frac{z_2 - z_b}{\tau_z} \right), \end{aligned} \right\} \quad (2.5)$$

with positive or negative values for coupling parameter $c$ and parameters $k_i, i = 1, 2$, adding possibly heterogeneity between both oscillators. In the following, as in [16], we will denote $(x_1, y_1, z_1)$-subsystem as $O_1$, which depends on $x_2$, and $(x_2, y_2, z_2)$-subsystem as $O_2$, which depends on $x_1$.

Let us describe the main behaviors patterns obtained in system (2.5). We distinguish two cases, depending if $k_i = 1, i = 1, 2$ that is, the homogeneous case, or $k_i \neq 1$ for some $i = 1, 2$, that is, the heterogeneous case.

### 2.2.1. Homogeneous case of the two-cells model

The case with $k_i = 1, i = 1, 2$, that is, the homogeneous case, was already studied in [16]. The authors have outlined the different synchronization types between $O_1$ and $O_2$ within the attractive limit cycle of system (2.5) depending on the coupling strength $c \neq 0$.

Here, we summarize the main asymptotic behaviors obtained in [16]: Figures 2a-4b illustrate each case. We set in every case the initial conditions of each oscillator as $(x_0, y_0, z_0) = [r, 4r - r^3, 1]$, with $r = 1.75$ for $O_1$ and $r = 1.25$ for $O_2$.

1. Total oscillation death: $O_1$ and $O_2$ reach asymptotically a stable equilibrium. Figure 2a shows the behavior of the system in the time interval $[0, 30]$ for $c = -0.7$.

2. Relaxation loss: one oscillator produces relaxation oscillations, i.e. the orbit approaches both $S_l$ and $S_r$, while the other one produces small oscillations and the orbit stays entirely near $S_l$. Figure 2b shows the behavior of the system in the time interval $[0, 30]$ for $c = -0.502$.

3. Antiphase synchronization: $O_1$ and $O_2$ reach asymptotically the same orbit, and the phase-shift between them is exactly one half-period. Figure 3a shows the behavior of the system in the time interval $[0, 30]$ for $c = -0.25$.

4. Uncoupled case: when $c = 0$ the oscillators $O_1$ and $O_2$ reach asymptotically the exact same orbit with MMOs. There is a phase shift between the oscillators, which depends on the initial conditions, see Figure 3b. The phase-shift will vanish if the initial conditions are identical.
<figure><img src="image_5.png" /><figcaption>(a) Total oscillation death, with $c = -0.7$. The orbit tends to a stable equilibrium.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(b) Relaxation loss, with $c = -0.502$. Only $O_1$ generates spikes in the corresponding ICC time series, $O_2$ does not undergo relaxation.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(a) Antiphase synchronization, with $c = -0.25$. Each spike of each oscillator occurs precisely at the middle time between two successive spikes of the other oscillator.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) No correlation between the cells, with $c = 0$. Both oscillators reach asymptotically the exact same orbit with MMOs.</figcaption></figure>

<caption>FIGURE 3. The *x* and *z* time series (top and bottom panel respectively) of each oscillator $O_1$ and $O_2$ are obtained from the two-cells model (2.5) for $k_1 = k_2 = 1$ and different coupling parameters $c \le 0$.</caption>

5. *Almost-in-phase synchronization:* system (2.5) admits an attractive limit cycle along which $O_1$ and $O_2$ generate MMOs. The two 3D oscillators are synchronized with a small phase-shift. Figure 4a shows the behavior of the system in the time interval [15, 45] for $c = 0.1$, as well as a zoom of one of the peaks where we can see that the oscillators are not perfectly in phase.
<figure><img src="image_9.png" /><figcaption>(a) Almost-in-phase synchronization, with $c = 0.1$. The two oscillators are synchronized with a slight phase-shift between their spikes.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) In-phase synchronization, with $c = 1$. The two oscillators spikes are perfectly synchronized. More precisely, the time difference between corresponding spikes of $O_1$ and $O_2$ tends asymptotically to 0.</figcaption></figure>

FIGURE 4. The *x* and *z* time series (top and bottom panel respectively) of each oscillator $O_1$ and $O_2$ are obtained from the two-cells model (2.5) for $k_1 = k_2 = 1$ and $c > 0$.

6. In-phase locking synchronization: $O_1$ and $O_2$ reach asymptotically the exact same orbit that is the precisely the cycle existing for the uncoupled case. Therefore, the attractive limit cycle of system (2.5) lies on $\{x_1 = x_2, y_1 = y_2, z_1 = z_2\}$ and is a global attractor of the system. Figure 4b shows the behavior of the system in the time interval [15, 45] for $c = 1$, as well as a zoom of one of the peaks.

In [16], the distribution of these behaviors according to the coupling parameter $c$ has been theoretically studied and numerically shown for symmetric coupled system (2.5) (with $k_1 = k_2$). Additionally, from a numerical investigation and heuristic numerical arguments, we have identified singular behaviors arising for very narrow intervals of $c$ values and corresponding to transitions between the main behaviors described above. For better understanding the organization of the behaviors and the role of the most singular ones, it is helpful to consider a larger case and erase the constraint of symmetry between the two oscillators. Despite the theoretical arguments based in this symmetry are no longer valid, the evanescent behaviors identified in [16] become main behaviors if we consider $k_1 \neq k_2$ and the roles of these timescale parameters and of the coupling parameter in organizing the dynamical structure are clarified by the numerical investigation presented in the next section.
<figure><img src="image_11.png" /><figcaption>FIGURE 5. Signature of the two-cells model (2.5) limit cycle according to *k* values for *c* = −0.05 (*left panel*), *c* = −0.25 (*central panel*) and *c* = −0.40 (*right panel*). The grey bands materialize the intervals of *k* values for which the signature is not well-defined.</figcaption></figure>

2.2.2. Two-cells model with heterogeneity in recovery dynamics

We introduce heterogeneity between the two oscillators by considering different values for the recovery variable times-scale parameters $k_1$ and $k_2$. Hence, if $k_1 = 1$ and $k_2 = k \neq 1$, the uncoupled oscillators ($O_1$ and $O_2$ with $c = 0$) admit attractive limit cycles with different intrinsic frequencies. Moreover, the oscillatory time series generated by the ICC variables $z_1$ and $z_2$ along the corresponding limit cycles feature different amplitudes. The theoretical study performed in [16] for the homogeneous case is strongly based on symmetry, therefore it cannot be directly adapted to this case. In this article, the study of the heterogeneous case is carried out using numerical simulations.

We introduce some notations that will be used from now on.

**Remark 2.1.** In the following, we say that the signature of a periodic orbit of system (2.5) is *s*/*l* if, along one period of this orbit, one 3D oscillator spikes *s* times while the other one spikes a single time. Hence, with this notation, the antiphase behavior of the homogeneous case (see Fig. 3a) would be denoted as 1/1.

In the following, we fix $k_1 = 1$ and consider $k_2 = k > 1$, since the results for $k < 1$ can be deduced by exchanging the roles of $O_1$ and $O_2$. The different simulations are performed as follows: for a fixed (either positive or negative) value of *c*, we solve the system for $k \in (1, k_0]$, where the upper limit $k_0$ changes depending on the case, with a step $dk = 0.01$ in *k* values and if one obtains a change in the behavior, we reduce the step to $dk = 0.0001$ to specify better the transition.

We further explain the synchronization patterns obtained in the heterogeneous case.

- For $c > 0$, while $k \in (1, 2]$, we obtain the same qualitative behavior as in the homogeneous case ($k = 1$).

- For $c < 0$, we find the three behaviors obtained in the homogeneous case, but in the antiphase region $c \in [c^0 + \delta, -\delta]$ with $c^0 \approx -0.5$ and $0 < \delta \ll 1$, apart from the antiphase behavior we obtain new synchronization patterns for $k \in (1, 3]$ that we explain subsequently.

In Figure 5, we represent the signature of system (2.5) limit cycle for $c = -0.05$ (left panel), $c = -0.25$ (central panel) and $c = -0.40$ (right panel) according to the value of parameter *k*. The blue lines represent the signature value when system (2.5) admits a limit cycle, while grey bands materialize the intervals of *k* value where the signature is not well-defined, corresponding to signature transitions. We note that the representation of the different behaviors in Figure 5 exhibits a Devil Staircase-like behavior. Devil Staircases appear in different fields and, in the case of dynamical systems, are characterized by different stable behaviors for large interval of one parameter values, separated by transition intervals where the system can exhibit a mixture of these stable behaviors, and even chaotic behaviors, see [2]. In our case, we have different stable behaviors (see Rem. 2.1) that are separated by transition intervals where the system exhibits a mixture of the stable behaviors among others that are discussed further in the following.
<figure><img src="image_12.png" /><figcaption>(a) Signature 1/1, with $k_2 = 1.5$. Each oscillator spikes a single time between two consecutive spikes of the other one.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(c) Signature 3/1, with $k_2 = 3$. Oscillator $O_2$ spikes precisely three times between two consecutive spikes of oscillator $O_1$.</figcaption></figure>

FIGURE 6. Phase-shifted synchronization with different signatures. The time series are generated with the two-cells model (2.5) with $k_1 = 1$ and varying $k_2$.

Let us focus on case $c = -0.25$ and describe the main and transitory patterns.

- For $k \in (1, 1.5563]$, we obtain a behavior qualitatively similar to the one obtained in the homogeneous case, one oscillation for one oscillation (signature 1/1), see Figure 6a for a specific instance where $k = 1.5$. However, as we can observe, the behavior is not antiphase, as it is the case in the homogeneous case (Fig. 3a), because the phase shift is not half of the period. This is an indicator of how the heterogeneity is affecting the system. We refer to the behavior shown in Figure 6a as “phase-shifted synchronization”.

- For $k \in [1.6468, 2.7350]$, one oscillator generates two pulses while the other one generates a single one along a limit cycle period, which corresponds to signature 2/1. See Figure 6b for a specific instance where $k = 2$.
<figure><img src="image_14.png" /><figcaption>FIGURE 7. Transitory signatures of system (2.5) outputs while varying $k$. Parameter $k_1 = 1$. Left panel: for $k = 1.6$, the signature alternates between 1/1 and 2/1. Right panel: for $k = 2.74$, the signature alternates between 2/1 and 3/1.</figcaption></figure>

- For $k \in [2.7582, 3]$, the behavior of the system changes and one oscillator generates three pulses while the other one generates a single one along a limit cycle period, which corresponds to signature 3/1. See Figure 6c for a specific instance where $k = 3$.

Naturally, any signature s/1 can be obtained by increasing the value of ratio $k$, and the same given signature is obtained for a whole interval of $k$ values. In the transition between two main behaviors, the system exhibits different combinations of both cases. We summarize here some of them.

- For $k \in [1.5592, 1.6454]$, in the transition from signature 1/1 to 2/1, the system can exhibit a mixture of two signatures: between successive spikes of $O_1$, oscillator $O_2$ spikes one and two times, alternatively. Therefore, two signatures 1/1 and 2/1 alternate along the limit cycle. See left panel of Figure 7.

- For $k \in [2.7372, 2.7563]$, in the transition from signature 2/1 to 3/1, a similar case occurs: between successive spikes of $O_1$, oscillator $O_2$ spikes two and three times, alternatively. Therefore, two signatures 2/1 and 3/1 alternate along the limit cycle. See Figure 7 right.

Note that there is some ranges of the parameter values inside the transitions where the system does not exhibit one of the behaviors described above, but a more complex combination of the main behaviors. In fact, we also have a canard explosion, that is, a very fast transition from small oscillations to relaxation oscillations, with a small variation of the parameters. This transition happens *via* a family of cycles whose characteristic property is that they contain *canard segments*, i.e. segments of the orbit that follow the repulsive slow manifold for a time $O(1)$ [6]. We observe canard orbits in Figure 8, for $k = 1.55698903$, and Figure 9, for $k = 1.55698902$. In the left panels is represented the $x_1y_1$-projection on the plane $z = z_b$ of the orbit and shows, in black dotted line the approximated $x_1$-nullcline, in dashed lines the limits of the $y_1$-nullcline and in blue a canard orbit without head (Fig. 8) or with head (Fig. 9) for the first oscillator. Right panels show the $x$ and $z$ components of both oscillators, $O_1$ in blue and $O_2$ in green.

The canard explosion takes place in the slow oscillator (in blue orbit in Figs. 8 and 9). In Figure 8, we observe a canard “without head” orbit and in Figure 9, a canard with head. Note that, between the two instances shown in these figures, only the 8th decimal digit of $k$ value changes from 3 to 2. One expects that this behavior is due to a folded-saddle singularity in the system [7, 18, 19, 25, 39, 43]. The detailed analysis of this fact, which can be tackled through desingularization techniques as those used in [13], is left for a future work.
<figure><img src="image_15.png" /><figcaption>FIGURE 8. Canard without head in $(x_1, z_1)$-plane (left panel) and resulting $x$ and $z$ time series for $O_1$ and $O_2$ (right panels). Parameter $k = 1.55698903$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIGURE 9. Canard without head in $(x_1, z_1)$-plane (left panel) and resulting $x$ and $z$ time series for $O_1$ and $O_2$ (right panels). Parameter $k = 1.55698902$.</figcaption></figure>

Now, that we have detailed the observed patterns in the case $c = -0.25$, let us come back to the general case. We study the behaviors for $k \in (1, 3]$ in the range of the coupling parameter $c \in [-0.4, 0)$, we obtain the subdomains in the $(c, k)$-plane where the system (2.5) exhibits each main behavior, as shown in Figure 10 right, where the region $c \in [-0.4, -0.05]$ and $k \in [1, 3]$ is displayed. We can observe that the transition intervals shrink while the coupling paramater values are lowered.

Let us focus on the first transition of the central panel of Figure 5. If we solve the system in the transition for different values of the heterogeneity parameter inside $[1.5, 1.7]$, we obtain Figure 11. We remark that although the signature is not well-defined in the transition we can associate an intuitive value for almost every behavior, as it is explained in the caption of Figure 11, so we can graphically see the Devil Staircase. In Figure 11, main behaviors are represented in blue, whereas the non well-defined signatures are represented in red. Although, we have reached the forth decimal of the heterogeneity parameter $k$ there are still intervals where we are not able to select an intuitive value for the signature, for instance, the canard behavior explained before.

Finally, we focus now on the range near the transition to the relaxation loss behavior, namely $c \in [-0.5, -0.4]$. All ranges of the main behaviors converge to the transition point for $k = 1$, see Figure 10 left showing
<figure><img src="image_17.png" /><figcaption>FIGURE 10. Regions of $(c, k)$-plane associated with fixed signature value of the attractive limit cycle of the two-cells model (2.5). The left panel shows that the intervals of $k$ values corresponding to a given signature converge to the singleton $\{1\}$ while $c$ decreases to $-0.5$. The right panel emphasizes the small strips in $(c, k)$-plane for which the signature is not well-defined.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>FIGURE 11. Devil Staircase in the signature with respect to $k$ value for $c = -0.25$. Outside the intervals of $k$-values for which the signature is well-defined, i.e. for $k \in [1.56, 1.65]$, a mean number of spikes generated by $O_2$ between two successive spikes of $O_1$ is calculated over a time interval $[0, 60]$ (red points and line). For a large subinterval of $k$ value, the signature alternates along time between 1/1 and 2/1, leading to a mean value 1.5 for $s$.</figcaption></figure>

$c \in [-0.5, -0.4]$ and $k \in [1, 3.5]$. We can also observe that the transition to the relaxation loss regime depends on the heterogeneity parameter $k$. The detailed analysis of this observed behavior together with the Devil's Staircase, which could be done by using similar techniques to those considered in [34], that is, the analysis of an appropriate transition map, is left to a future work.
## 2.3. Network model

In more realistic problems, we have more than two cells. We consider here the case where we have $N$ different cells and symmetric coupling, so we can model the system with $N$ coupled oscillators:

$$ O_i \left\{ \begin{array}{l} \dot{x}_i = \tau(-y_i + f(x_i) - \phi_f(z_i)), \\ \dot{y}_i = \tau\varepsilon k_i(x_i + a_1 y_i + a_2 + \frac{1}{N/2} \sum_{j=1}^{N} c_{ij}(x_i - x_j)), \\ \dot{z}_i = \tau\varepsilon \left( \phi_r(x_i) - \frac{z_i - z_b}{\tau_z} \right), \end{array} \right. \qquad (2.6) $$

where $i = 1, 2, \dots, N$ and $c_{ij} \in [-1, 1]$ with $c_{ij} = c_{ji}$. For $N = 2$, we recover system (2.5).

In a first approach, we consider the cells separated in two groups or clusters, namely first cluster, are the cells with indexes $i = 1, \dots, N_-$, and second cluster, the cells with indexes $i = N_- + 1, \dots, N$. A similar approach was first introduced for the zebrafish embryo in [15]. We set $N_- = \lfloor N/2 \rfloor$, that is, the greatest integer lower than or equal to $N/2$, and $N_+ = \lceil N/2 \rceil$, that is, the lowest integer greater than or equal to $N/2$.

We set that the cells of the same cluster are coupled with $c_{ij} = c_\alpha = 1$, i.e. in-phase synchronization and that every two cells of different clusters are coupled with the same parameter $c_\beta$, that is $c_{ij} = c_\beta$, $\forall i = 1, \dots, N_-, \forall j = N_- + 1, \dots, N$. We consider two different cases:

- Homogeneous case ($k_i = 1$). All the cells have the same intrinsic amplitude and frequency.

- Heterogeneous case ($k_i = 1$ for $i = 1, \dots, N_-$ and $k_i > 1$ for $i = N_- + 1, \dots, N$). All the cells in one cluster have the same intrinsic amplitude and frequency but the cells of the other cluster have different intrinsic amplitudes and frequencies.

## 3. NUMERICAL METHOD DESCRIPTION

As said in the introduction, although high-fidelity solvers approximate the solution, in some cases it is computationally excessive to solve the system, so we introduce a ROM approximation [9, 26, 27, 31, 41]. This seeks to define an equivalent dynamical system of lower dimension that approximates the full order model, namely, the reduced problem.

Although methods presented in the literature [9, 26, 27, 31, 41], reduce the number of equations of the system, in some cases there is no translation in a reasonable reduction in the computational time. In this section, we develop the reduced problem in a way we obtain a reasonable speed up in the computational times.

In order to build the reduced problem, there are two main steps:

1. Construction of the reduced basis. First, we need to determine the reduced basis dimension. To do so, we apply the Proper Orthogonal Decomposition (POD), the most popular approach for reduced basis formulation [9, 27, 31, 41]. The POD technique gives us the basis of the lower dimensional space.

2. Construction of the reduced problem. Once we have a basis, we can formulate a reduced problem equivalent to the full order problem.

We remark here that methods in the literature [9, 27, 31, 41] only seek for a projector onto the lower dimensional space. Our approach is slightly different as we try to define the solution of the original problem as the linear combination of some (and not many) solutions of a reduced problem. Before explaining how we construct the ROM once we have obtained a basis, we need to represent the original problem in an algebraic form.
### 3.1. Matrix representation of the original problem

In this section, we build the algebraic form of the original network problem (2.6). We name $\boldsymbol{w} = (x_1, y_1, z_1, \dots, x_N, y_N, z_N)$, and we can write the original problem (2.6) as

$$ \dot{\boldsymbol{w}} = f(\boldsymbol{w}, t), \text{ with } \boldsymbol{w}(0) = \boldsymbol{w}_0. $$

We can rewrite the system as

$$ \dot{\boldsymbol{w}} = (A + B(c_{\alpha}, c_{\beta}))\boldsymbol{w} + \boldsymbol{b} + \boldsymbol{b}_{nl}(\boldsymbol{w}), $$

where

$$ A = \boldsymbol{k} \begin{pmatrix} A_l & & \\ & \ddots & \\ & & A_l \end{pmatrix}, \text{ with } A_l = \tau \begin{pmatrix} 4 & -1 & 0 \\ \varepsilon & \varepsilon a_1 & 0 \\ 0 & 0 & -\varepsilon/\tau_z \end{pmatrix}, $$

and

$$ \frac{\tau \varepsilon}{N/2} \boldsymbol{k} \cdot \begin{pmatrix} -\Sigma_1 M & c_\alpha M & \dots & c_\alpha M & c_\beta M & \dots & \dots & c_\beta M \\ c_\alpha M & \ddots & \ddots & \vdots & \vdots & & & & \vdots \\ \vdots & \ddots & \ddots & c_\alpha M & \vdots & & & & \vdots \\ c_\alpha M & \dots & c_\alpha M & -\Sigma_{N-1} M & c_\beta M & \dots & \dots & c_\beta M \\ c_\beta M & \dots & \dots & c_\beta M & -\Sigma_{N-1+1} M & c_\alpha M & \dots & c_\alpha M \\ \vdots & & \vdots & & c_\alpha M & \ddots & \ddots & \vdots \\ \vdots & & \vdots & & \vdots & & \ddots & c_\alpha M & \vdots \\ c_\beta M & \dots & \dots & c_\beta M & c_\alpha M & \dots & c_\alpha M & -\Sigma_N M \end{pmatrix}, $$

$$ \text{with } M = \begin{pmatrix} 0 & 0 & 0 \\ -1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \quad \boldsymbol{k} = \begin{pmatrix} 1 & 0 & 0 & \dots & 0 & 0 & 0 \\ 0 & k_1 & 0 & \dots & 0 & 0 & 0 \\ 0 & 0 & 1 & \dots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \dots & 1 & 0 & 0 \\ 0 & 0 & 0 & \dots & 0 & k_N & 0 \\ 0 & 0 & 0 & \dots & 0 & 0 & 1 \end{pmatrix} \quad \text{and} $$

$$ \Sigma_i = \begin{cases} (N_- - 1)c_\alpha + N_+c_\beta, & \text{for } i = 1, \dots, N_-, \\ (N_+ - 1)c_\alpha + N_-c_\beta, & \text{for } i = N_- + 1, \dots, N, \end{cases} $$

where $N_-$ and $N_+$ have been defined in Section 2.3. Thus, top left block is $3N_- \times 3N_-$, top right is $3N_- \times 3N_+$, bottom left is $3N_+ \times 3N_-$, and bottom right is $3N_+ \times 3N_-$. We remark here that matrix $B(c_\alpha, c_\beta)$ admits an affine decomposition depending on $c_\alpha, c_\beta$ and $\Sigma_i$.
$$b = k \begin{pmatrix} b_0 \\ \vdots \\ b_0 \end{pmatrix}, \text{ with } b_0 = \tau \varepsilon \begin{pmatrix} 0 \\ a_2 \\ z_b / \tau_z \end{pmatrix}.$$

$$b_{nl}(\mathbf{w}) = \tau \begin{pmatrix} d_1(\mathbf{w}) \\ \vdots \\ d_N(\mathbf{w}) \end{pmatrix}, \text{ with } d_i(\mathbf{w}) = \begin{pmatrix} -x_i^3 - \phi_f(z_i) \\ 0 \\ \varepsilon \phi_r(x_i) \end{pmatrix}, i = 1, \dots, N.$$

## 3.2. Reduced problem

In order to build the reduced problem, we aim at expressing the original variable $\mathbf{w}$ as a linear combination of certain functions. Following the procedure of the POD technique we obtain a basis of data set, called snapshots, so we can write the original variable as

$$\mathbf{w} = \sum_{i=1}^{N_s} a_i \phi_i,$$

where $N_s$ is the number of snapshots considered. The procedure is slightly different from the one presented in [31] since we do not use straightforwardly the eigenvectors of the correlation matrix as the vector of the basis of the reduced space. Instead, we compute the left-multiplication by the snapshot matrix and we perform a division by a constant depending on the corresponding eigenvalue and the number of snapshots to ensure the orthonormality. Therefore, we obtain the eigenvectors of the snapshot matrix and we can approximate the original solution as a combination of these. We remark here that the vector $\mathbf{a} = (a_i, \dots, a_{N_s})$ depends on the time, whereas the vector of the basis does not.

Each of these basis functions (or vectors) is associated with a certain weight. These weights are the singular values of the snapshots matrix, which are a measure of the information of the original system solution each basis function retain. If we sort the singular values in descending order, we can finish to apply the POD technique to select a reduced number of basis functions so that the information retained by them is maximized, and we can approximate the original variable by

$$\mathbf{w} \approx \sum_{j=1}^{N_{rb}} a_j \phi_j = \underline{\phi} \mathbf{a}, \qquad (3.1)$$

where $\underline{\phi} = (\phi_1 | \dots | \phi_{N_{rb}})$ and $\mathbf{a} \in \mathbb{R}^{N_{rb}}$. Introducing approximation (3.1) into the algebraic form and naming $D = \bar{A} + B(c_\alpha, c_\beta)$, we obtain,

$$\underline{\phi} \dot{\mathbf{a}} = D \underline{\phi} \mathbf{a} + \mathbf{b} + \mathbf{b}_{nl}(\mathbf{w}),$$

now we make use of the fact that $\underline{\phi}^T \underline{\phi} = 1_{N_{rb}}$, as the vector of the basis are orthonormal. We get to

$$\dot{\mathbf{a}} = \underline{\phi}^T D \underline{\phi} \mathbf{a} + \underline{\phi} \mathbf{b} + \underline{\phi} \mathbf{b}_{nl}(\mathbf{w}).$$
We name now $\hat{D} = \phi^T D \phi$, $\hat{b} = \phi \hat{b}$, and to deal with the nonlinearity involved in the last term, we apply the Empirical Interpolation Method (EIM) to approximate it. Finally, we obtain:

$$ \dot{a} = \hat{D}a + \hat{b} + EIM(b_{nl}). $$

To sum up, we have started from a system with $3N$ equations and we have built an equivalent reduced problem of $N_{rb}$ equations where the nonlinear terms are approximated by an EIM. Therefore, in problems that require a high number of evaluations for different instances of the parameters, we obtain $N_{rb} \ll 3N$ and the reduction is drastic.

## 4. NUMERICAL TESTS

In this section, we test the numerical method described in the previous section for different cases with various number of cells, focusing on the antiphase synchronization interval. The main objectives of this section are the following:

- Check that the reduced systems built via the reduction method can reproduce the behaviors of the original systems, although the original systems are stiff and features a slow-fast structure.

- Apply the reduction to some cases where the theoretical analysis is too difficult to handle directly.

The section is divided as follows: in a first part, we perform different tests to validate our method with the aim to achieve the first goal. Once our method is validated, we focus on the heterogeneous model (2.5) with two cells and $k_1 = 1$ and $k_2 = k > 1$. Finally, we consider the network model (2.6) with a high number of cells.

Following the information of Appendix A, for system (2.5) we select initial conditions $(x_0, y_0, z_0) = [r, 4r - r^3, 1]$, with $r = 1.75$ for $O_1$ and $r = 1.25$ for $O_2$, and we solve the system along time interval $[0, 30]$. Throughout the section, time is measured in seconds, so we omitted the notation.

### 4.1. Validation tests

For the first test, we compare the solutions of various instances of the homogeneous case, i.e. system (2.5) with $k_1 = k_2 = 1$, and check that the solution of our reduced problem truly reproduces the original solution with a reasonable reduction in needed CPU times. We are using the same adaptive solver for integrating both original and reduced systems.

In order to be able to compare the results, we impose that the reduced basis dimension is the dimension of the original system, i.e. $N_{rb} = 3N = 6$ and the nonlinear part is approximated by EIM with the highest possible number of vectors, in this case 2 for each. We study the CPU-times needed to solve a hundred instances of the same problem, represented in Table 1.

<table><caption>TABLE 1. Comparison of CPU times.</caption><thead><tr><th rowspan="2">c</th><th colspan="3">Original</th><th colspan="3">Reduced</th></tr><tr><th>T = 10</th><th>T = 20</th><th>T = 30</th><th>T = 10</th><th>T = 20</th><th>T = 30</th></tr></thead><tbody><tr><td>-0.05</td><td>21.32</td><td>50.77</td><td>120.11</td><td>1.41</td><td>3.62</td><td>10.96</td></tr><tr><td>-0.10</td><td>21.29</td><td>50.66</td><td>120.45</td><td>1.40</td><td>3.65</td><td>11.10</td></tr><tr><td>-0.15</td><td>21.29</td><td>50.78</td><td>120.38</td><td>1.38</td><td>3.64</td><td>11.06</td></tr><tr><td>-0.20</td><td>21.36</td><td>50.77</td><td>120.43</td><td>1.42</td><td>3.63</td><td>11.05</td></tr><tr><td>-0.25</td><td>21.32</td><td>50.66</td><td>120.32</td><td>1.41</td><td>3.64</td><td>11.06</td></tr><tr><td>-0.30</td><td>21.32</td><td>50.64</td><td>120.22</td><td>1.41</td><td>3.63</td><td>11.10</td></tr><tr><td>-0.35</td><td>21.33</td><td>50.76</td><td>120.17</td><td>1.42</td><td>3.63</td><td>11.03</td></tr><tr><td>-0.40</td><td>21.33</td><td>50.63</td><td>120.45</td><td>1.42</td><td>3.62</td><td>11.03</td></tr><tr><td>-0.45</td><td>21.30</td><td>50.79</td><td>120.33</td><td>1.42</td><td>3.64</td><td>11.02</td></tr><tr><td>Mean</td><td>21.32</td><td>50.72</td><td>120.32</td><td>1.41</td><td>3.64</td><td>11.05</td></tr></tbody></table>
<figure><img src="image_19.png" /><figcaption>FIGURE 12 Comparison between the *x* and *z* time series of $O_1$ and $O_2$ obtained from the reduced systems in case of antiphasic synchronization.</figcaption></figure>

We note that the CPU times does not depend on the value of the coupling parameter *c* inside the antiphase synchronization regime. However, the reduced problem requires less CPU time to be solved, with a speed-up around 15 for *T* = 10 and *T* = 20 and slightly over 10 for *T* = 30. We remark that although the number of equations is not reduced, the computational times is consistently lowered.

Figure 12 compares *x* and *z* time series generated by the oscillators between the original and reduced models for coupling parameter value *c* = -0.25. The reduced problem reproduces truthfully the original solution.

At this point, the tests are finished in a positive manner, the reduced problem reproduces the original solution while performing a reduction in the CPU times, so our reduction procedure is validated. Although this procedure only allows us to build a reduced problem for a fixed value of *c*, it is still useful when trying to catch a certain behaviour for a particular instance of the coupling parameter, for instance, the canard phenomena which we reproduce in Section 4.2, and it is also the first step in the construction of a ROM for a range of the parameters, which is left for a future work.

## 4.2. ROM in time for the heterogeneous model (2.5)

We now apply the reduction method to the heterogeneous system (2.5) with $k_1 = 1$ and $k_2 = k > 1$. The goal is to reproduce the canard phenomena of the original problem, but with the reduced model built by a reduction in time. We proceed as follows: we first study the reduced system behavior for values of *k* of the two first main behaviours of the central panel in Figure 5, we set *k* = 1.5 for signature 1/1 and *k* = 2 for signature 2/1. Then, we study the solution for a value of *k* that lies in the transition interval, in particular *k* = 1.6, that exhibits a mixture of both signatures and finally we try to reproduce the canard phenomena for *k* = 1.55698902 and *k* = 1.55698903.

We set the same initial conditions as in the homogeneous case (Sect. 4.1), and we solve the system in the time interval [0, 30] (see A) for the coupling parameter *c* = -0.25. We are following the procedure described in Section 3 and we impose that the reduced model has the same number of equations as the original system.

As shown in Figure 13a and b, the reduced model reproduces the new behaviors identified from the heterogeneous system (2.5) outputs: signature 1/1 for *k* = 1.5 in the left panel and signature 2/1 for *k* = 2 in the right panel.
<figure><img src="image_20.png" /><figcaption>FIGURE 13. Comparison between the *x* and *z* time series of *O*₁ and *O*₂ obtained from the reduced systems in case of phase-shifted synchronization with different signatures.</figcaption></figure>

We conclude that our reduced model still exhibits the new main behaviours obtained for the heterogeneous system (2.5), so we proceed with the study of the transition interval. We consider a value of *k* inside the transition interval, for which the system does not feature one of the main behaviours (with well defined signature), but a mixture of signatures or even canard phenomena. For *k* = 1.6, Figure 13c shows that the reduced system also generates the same mixed behaviors as the original system, with an alternation of two signatures 1/1 and 2/1.

Finally, we compute the reduced system for $k = 1.55698902$ and $k = 1.55698903$ with the aim of reproducing canard phenomena. Figure 14a and b displays such canard phenomena, where the orbit follow for a time $O(1)$ the repulsive slow manifold: canard without head in Figure 14a and canard with head in Figure 14b.

### 4.3. ROM in time for the two clusters network model (2.6)

The ultimate goal of this work is to build a reduced model that we can compute faster than the original one for a large number of cells while representing its behaviors, since theoretical analysis of this original model is still hardly tractable. In this section, we aim at applying a reduction in time to system (2.6).
<figure><img src="image_21.png" /><figcaption>FIGURE 14. Comparison between the *x* and *z* time series of *O*₁ and *O*₂ obtained from the reduced systems in case of the transition behavior involving a canard.</figcaption></figure>

In a first approach, we build a reduced model for an instance in which we can heuristically know how the system behaves. For this instance, we have considered homogeneous, $k_i = 1$, system (2.6) for $N = 100$ cells, so we have 300 equations separated in two clusters. We consider the same initial conditions $[r, 4r - r^3, 1]$ for every cell of each cluster, with $r = -1.25$ for the first cluster and $r = -1.75$ for the second cluster, and the coupling parameters to be $c_\alpha = 1$ for cells inside the same cluster and for cells in different clusters $c_\beta$ fixed in the range antiphase synchronization regime, more precisely we select $c_\beta = -0.25$. Under these assumptions, we have copies of the same cell in each cluster and the system behaves as two coupled cells.

We consider a second case dealing with the homogeneous system (2.6) but with slightly different initial conditions in each cluster: for every cell, we choose an initial conditions $[r, 4r - r^3, 1]$ with $r \in [-1.25, -1.2]$ for the cells within the first cluster and $r \in [-1.8, -1.75]$ for the cells within the second one.

For a third approach, we consider that cells within the second cluster are heterogeneous, so we select a random vector for the second cluster $\mathbf{k} = [k_{N-1} + 1, \dots, k_N]$ in a certain range where we assume the system behaves in a qualitatively similar way for a fixed interval of the coupling parameter. We select the range in the heterogeneity
<figure><img src="image_22.png" /><figcaption>FIGURE 15. Relative sum of neglected singular values for different cases.</figcaption></figure>

parameter to be [1.7, 1.8], assuming the system would behave as the two cells case, see Figure 10 right. We consider the same initial conditions as in the first case.

It is worth noting that our method has been developed for capturing efficiently the asymptotic behaviors. Indeed, to avoid long transients in the second and third cases, we have selected the initial conditions as explained in Appendix A.1. Therefore, the following simulations show the robustness of the method in case of weak heterogeneity in initial conditions. We note that considering a wider range of initial conditions would result in a less efficient reduction of the model.

We set the integration interval as [0, 30]. Following our procedure, based on the POD technique, we represent in Figure 15 the relative sum of neglected singular values versus the number of selected singular values for the first 30 singular values and for the cases described above. In every case we apply the following criterion for the selection of the reduced model dimension: we select the reduced basis dimension as the lowest integer for which the relative sum of the neglected singular values is below a certain tolerance, in our case, we have selected a tolerance equal to $10^{-6}$. In Figure 15 the tolerance is represented by the straight line.

For the first case, as we can see in Figure 15, the reduced model corresponding to this case has only 6 equations, the equivalent of two cells, as we expected. We note here that for this simple case the reduction carried out is from 300 equations to just 6, a 2% of the original equations, which is a very good reduction in our case. Speaking of the speed-up in the CPU time, we obtain a speed-up of $\times 35.30$ comparing the CPU times of solving the system in the time interval [0, 30]. Comparison of other time intervals are summarized in Table 2.

<table><caption>TABLE 2. Comparison of speed-ups for the cluster model (2.6).</caption><thead><tr><th>Case</th><th><i>T</i> = 10</th><th><i>T</i> = 20</th><th><i>T</i> = 30</th></tr></thead><tbody><tr><td>1</td><td>89.80</td><td>66.15</td><td>35.30</td></tr><tr><td>2</td><td>78.42</td><td>61.37</td><td>32.23</td></tr><tr><td>3</td><td>56.03</td><td>43.29</td><td>23.10</td></tr></tbody></table>

In Figure 15, we can see that in the second case we need a higher dimension as we need to select 9 singular values to be below the stopping tolerance. So, the reduced model would have 9 equations, this is a reduction to just 3% of the equations of the original problem. The speed-up in the CPU time for the time interval [0, 30]
<figure><img src="image_23.png" /><figcaption>FIGURE 16. Comparison between the *x*-components of the reduced system, top, and of the original network system (2.6), bottom, for an instance of case 2 where *N* = 100. Cells within the first cluster are set in blue and the second cluster are set in green. *Left panel:* Whole time interval. *Right panel:* Zoom at the beginning of the time interval.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>FIGURE 17. Comparison between the *x*-components of the reduced system, top, and of the original network system (2.6), bottom, for an instance of case 3 where *N* = 10. Cells within the first cluster (in blue) are homogeneous and the range of the heterogeneity parameter in the second cluster (in green) is set to [1, 1.5].</figcaption></figure>

is ×32.23. Comparison of other time intervals are summarized in Table 2. In Figure 16 we show a comparison between time-series obtained with the original homogeneous network system and the associated reduced system when considering different initial conditions for each cell. First, we observe that the network cells synchronize in each cluster very fast, and therefore the network systems approaches quickly its asymptotic behavior. Second, the ROM method is particularly efficient to approximate such limit cycle, since it is somehow dedicated to this aim. Third, we can observe that the reduced system also approximates accurately the transient behavior even if the error (at each time-step) during the short interval of time at the beginning of the simulations is slightly larger than it is once the orbits becomes close to the limit cycle.
In the third case, as can be seen in Figure 15, we need to consider 15 singular values to lay under the tolerance, that means that the reduced model has 15 equations, that is, 5% of the equations of the original system. The speed-up in this case for the time interval [0, 30] is ×23.10. Comparison of other time intervals are summarized in Table 2. For illustration purposes, we present in Figure 17 an instance of this case with N = 10 and the range of the parameter as [1, 1.5]. We can see that while in the first cluster every cell behaves in the same way, every cell in the second cluster has a different frequency from each other. We can observe that the reduced model reproduces truly these differences.

5. CONCLUSIONS AND FUTURE WORK

In this study, we have extended the numerical analysis of the model introduced in [16] in the case of two
non identical coupled ICC oscillators, precisely with different value of the timescale parameter for the recov-
ery variable k₁ and k₂. In this case, we have identified new types of synchronization patterns, in particular
asymmetric antiphasic synchronization and, more generally, the ability of one oscillator to spike several times
between two successive spikes of the other one. We have introduced the notion of signature of an asymptoti-
cally stable oscillatory pattern and that the signature follows a Devil's Staircase organization when varying the
value of k₂/k₁ ratio. The role of canard phenomenon in the transition between k₂/k₁ intervals corresponding to
well-defined and constant signature of the oscillations has been emphasized and numerically shown. Moreover,
we have shown the convergence of all the intervals k₂/k₁ associated with fixed signature to the singleton {1}
when the (negative) coupling strength is decreased down to the value corresponding to the transition towards
relaxation loss (of one oscillator). This value k₂/k₁ = 1 naturally splits the two cases (i) O₁ relaxation loss
(k₁ < k₂) and (ii) O₂ relaxation loss (k₁ > k₂).

We have also introduced a network model using the same coupling type considered for the two-cells model. We have especially considered the connectivity matrix between cells associated with the case of two populations, each of them formed by excitatory coupling of neurons, and inhibiting each other. We have adapted a ROM in time that, from a given instance of the network model, builds a system that approximates accurately the original one and requires much less time to be numerically integrated. We have tested the accuracy of the ROM by comparing simulations on the original and reduced system in various cases serving as benchmarks and inspired by the first part of the study: (i) two identical cells inhibitory coupling leading to symmetric antiphasic synchronization, (ii) two cells inhibitory coupling with different parameter values leading to phase-shifted synchronization (iii) canard phenomenon both with and without head (iv) network model (with 100 cells) of two populations inhibiting each others. We have shown both the ability of the system built by the ROM to reproduce accurately the outputs of the original system, even the slow-fast transitions of the orbits and the canard phenomena involved in the signature transitions. Let us emphasize that the numerical tracking of canard is a non trivial problem since they appear for very narrow interval of parameter. We have finally shown the drastic reduction of CPU needed to integrate the reduced system compared to the original one, especially for large cells number.

The ROM in time is therefore a useful tool for numerical investigation of such network model since it captures
efficiently the complex behaviors of the coupled model and well approximate the quantitative data of the time
series. The theoretical investigation of synchronization in the two-cells model with different parameter values and
in the network model is a hard but tremendously interesting problem. We may benefit from efficient numerical
tracking of precise phenomena – such as canard explosion, signature transitions, loss of relaxation – that is
consequently accelerated by application of the ROM.

In future studies, we intend to analyze the organization of the dynamics close to the folded singularity of the
6D systems and how it underlies the Devil’s Staircase of the signature changes with k₂/k₁ parameter value. We
also intend to analyze the repartition of the network behaviors – more precisely, the type of synchronization
intra- and inter-clusters – according to the structure of the connectivity matrix. Furthermore, we intend to
perform ROM in both time and coupling parameter for the network model in order to reduce the simulation
times.
<figure><img src="image_25.png" /><figcaption>FIGURE A.1. Transition time versus the absolute difference between the *x* initial components for various values of *c*. From left to right, *c* = -0.05, *c* = -0.25, *c* = -0.45.</figcaption></figure>

# APPENDIX A.

We need to fix a set of initial conditions and the integration interval, so before starting we perform a little study of the dependence of the solution on the initial conditions and the dependence of the period between peaks of the same cell on the coupling parameter *c*.

## A.1 Initial conditions selection

We compute the time spent before the first pulse, namely transition time, in the homogeneous case for different values of *c* and for different selections of initial conditions. We set: $x_0 = [r_1, r_2, 4r_1 - r_1^3, 4r_2 - r_2^3, 1, 1]$, where $r_1$ and $r_2$ are random numbers over [-2, -1].

In Figure A.1, we represent the transition time versus the absolute difference between $r_1$ and $r_2$ for different values of *c*.

We see that for those initial conditions with $|r_1 - r_2| \ge 0.5$, the transition time is less than half the period, see Figure A.1 and Table A.1.

<table><caption>TABLE A.1. Numerical periods for various values of *c*.</caption><thead><tr><th>c</th><th>T<sub>1</sub></th><th>T<sub>2</sub></th><th>Period</th></tr></thead><tbody><tr><td>-0.45</td><td>2.820</td><td>2.447</td><td>10.534</td></tr><tr><td>-0.4</td><td>2.268</td><td>2.080</td><td>8.696</td></tr><tr><td>-0.35</td><td>1.925</td><td>1.846</td><td>7.542</td></tr><tr><td>-0.3</td><td>1.682</td><td>1.687</td><td>6.738</td></tr><tr><td>-0.25</td><td>1.509</td><td>1.557</td><td>6.132</td></tr><tr><td>-0.2</td><td>1.368</td><td>1.461</td><td>5.658</td></tr><tr><td>-0.15</td><td>1.249</td><td>1.395</td><td>5.288</td></tr><tr><td>-0.1</td><td>1.157</td><td>1.343</td><td>5.000</td></tr><tr><td>-0.05</td><td>1.074</td><td>1.327</td><td>4.802</td></tr></tbody></table>

## A.2 Time interval determination

Taking into account that the maximum period obtained is near 10, see Table A.1 and that the transition time is less than the half of it for our selection of initial conditions, see Figure A.1, we select an integration interval at least $[t_0, t_f] = [0, 30]$. Therefore, in the worst case we catch 2 periods.

We can compute the phase durations numerically obtaining the time difference between the maximum peaks of each cell and then obtaining the mean value. The obtained results are resumed in Table A.1. We can clearly see that as *c* decreases the period increases, as proved in [16].
*Acknowledgements.* This work has been supported by the Spanish Government Project RTI2018-093521-B-C31.

## REFERENCES

[1] A. Arenas, A. Diaz-Guilera, J. Kurths, Y. Moreno and C. Zhou, Synchronization in complex networks. *Phys. Rep.* **469** (2008) 93–153.

[2] P. Bak, Devil's staircase. *Phys. Today* **39** (1986) 38–45.

[3] V.N. Belykh, I.V. Belykh and M. Hasler, Connection graph stability method for synchronized coupled chaotic systems. *Physica D* **195** (2004) 188–206.

[4] V.N. Belykh and E.V. Pankratova, Chaotic synchronization in ensembles of coupled neurons modeled by the FitzHugh-Rinzel system. *Radiophys. Quantum Electr.* **4** (2006) 910–921.

[5] J. Bélair and P. Holmes, On linearly coupled relaxation oscillations. *Quart. Appl. Math.* **42** (1984) 193–219.

[6] E. Benoit, J. Callot, F. Diener and M. Diener, Chasse au canard (première partie). *Collectanea Math.* **32** (1981) 37–76.

[7] M. Brøns, M. Krupa and M. Wechselberger, Mixed mode oscillations due to the generalized canard phenomenon. *Fields Inst. Commun.* **49** (2006) 39–63.

[8] S.A. Campbell and M. Waite, Multistability in coupled FitzHugh-Nagumo oscillators. Science Direct Working Paper (2001) (S1574-0358), 04.

[9] D. Chapelle, A. Gariah and J. Sainte-Marie, Galerkin approximation with proper orthogonal decomposition: new error estimates and illustrative examples. *ESAIM: M2AN* **46** (2012) 731–757.

[10] J. Drover, J. Rubin, J. Su and B. Ermentrout, Analysis of a canard mechanism by which excitatory synaptic coupling can synchronize neurons at low firing frequencies. *SIAM J. Appl. Math.* **65** (2004) 69–92.

[11] B. Ermentrout, M. Pascal and B. Gutkin, The effects of spike frequency adaptation and negative feedback on the synchronization of neural oscillators. *Neural Comput.* **13** (2001) 1285–1310.

[12] B. Ermentrout and M. Wechselberger, Canards, clusters, and synchronization in a weakly coupled interneuron model. *SIAM J. Appl. Dyn. Syst.* **8** (2009) 253–278.

[13] E.K. Ersöz, M. Desroches, M. Krupa and F. Clément, Canard-mediated (de) synchronization in coupled phantom bursters. *SIAM J. Appl. Dyn. Syst.* **15** (2016) 580–608.

[14] E.K. Ersöz, M. Desroches and M. Krupa, Synchronization of weakly coupled canard oscillators. *Physica D* **349** (2017) 46–61.

[15] F.D.V. Fallani, M. Corazzol, J.R. Sternberg, C. Wyart and M. Chavez, Hierarchy of neural organization in the embryonic spinal cord: Granger-causality graph analysis of in vivo calcium imaging data. *IEEE Trans. Neural Syst. Rehab. Eng.* **23** (2014) 333–341.

[16] S. Fernández-García and A. Vidal, Symmetric coupling of multiple timescale systems with mixed-mode oscillations and synchronization. *Physica D* **401** (2020) 132129.

[17] R. FitzHugh, Impulses and physiological states in theoretical models of nerve membrane. *Biophys. J.* **1** (1961) 445–466.

[18] J. Guckenheimer, Singular Hopf bifurcation in systems with two slow variables. *SIAM J. Appl. Dyn. Syst.* **7** (2008) 1355–1377.

[19] J. Guckenheimer, and P. Meerkamp, Unfoldings of singular Hopf bifurcation. *SIAM J. Appl. Dyn. Syst.* **11** (2012) 1325–1359.

[20] A.L. Hodgkin and A.F. Huxley, A quantitative description of membrane current and its application to conduction and excitation in nerve. *J. Physiol.* **117** (1952) 500–544.

[21] E.M. Izhikevich, Phase equations for relaxation oscillators. *SIAM J. Appl. Math.* **60** (2000) 1789–1804.

[22] K. Jahn, J. Grosskreutz, K. Haastert, E. Ziegler, F. Schlesinger, C. Grothe and J. Bufler, Temporospatial coupling of networked synaptic activation of AMPA-type glutamate receptor channels and calcium transients in cultured motoneurons. *Neuroscience* **142** (2006) 1019–1029.

[23] M. Krupa, N. Popović, N. Kopell and H. Rotstein, Mixed-mode oscillations in a three time-scale model for the dopaminergic neuron. *Chaos* **18** (2008) 015106.

[24] M. Krupa, A. Vidal and F. Clément, A network model of the periodic synchronization process in the dynamics of calcium concentration in GnRH neurons. *J. Math. Neurosci.* **3** (2013) 1–24.

[25] M. Krupa and V. Wechselberger, Local analysis near a folded saddle-node singularity. *J. Differ. Equ.* **248** (2010) 2841–2888.

[26] T. Kostova-Vassilevska and G.M. Oxberry, Model reduction of dynamical systems by proper orthogonal decomposition: error bounds and comparison of methods using snapshots from the solution and the time derivatives. *J. Comput. Appl. Math.* **330** (2018) 553–573.

[27] K. Kunisch and S. Volkwein, Galerkin proper orthogonal decomposition methods for parabolic problems. *em Numer. Math.* **90** (2001) 117–148.

[28] E. Lee and D. Terman, Stable antiphase oscillations in a network of electrically coupled model neurons. *SIAM J. Appl. Dyn. Syst.* **12** (2013) 1–27.

[29] E. Lee and D. Terman, Stability of antiphase oscillations in a network of inhibitory neurons. *SIAM J. Appl. Dyn. Syst.* **14** (2015) 448–480.

[30] J. Nagumo, S. Arimoto and S. Yoshizawa, An active pulse transmission line simulating nerve axon. *Proc. IRE* **50** (1962) 2061–2070.

[31] M. Rathinam and L.R. Petzold, A new look at proper orthogonal decomposition. *SIAM J. Numer. Anal.* **41** (2003) 1893–1925.

[32] R. Reimbayev, K. Daley and I.V. Belykh, When two wrongs make a right: synchronized neuronal bursting from combined electrical and inhibitory coupling. *Phil. Trans. R Soc. A* **375** (2017) 20160282.
[33] H.G. Rotstein and R. Kuske, Localized and asynchronous patterns via canards in coupled calcium oscillators. *Physica D* **215** (2006) 46–61.

[34] J.E. Rubin, J. Signerska-Rynkowska, J. Touboul and A. Vidal, Wild oscillations in a nonlinear neuron model with resets:(II) Mixed-mode oscillations. *Discrete Continu. Dyn. Syst.* **B** **22** (2017) 4003–4039.

[35] A. Stefański, Determining thresholds of complete synchronization, and application. *World Sci. Ser. Nonlinear Sci.* **A** **67** (2009).

[36] D.W. Storti and R.H. Rand, Dynamics of two strongly coupled relaxation oscillators. *SIAM J. Appl. Math.* **46** (1986) 56–67.

[37] D. Storti and R.H. Rand, A simplified model of coupled relaxation oscillators. *Int. J. Non-linear Mech.* **22** (1987) 283–289.

[38] C. Stosiek, O. Garaschuk, K. Holthoff and A. Konnerth, In vivo two-photon calcium imaging of neuronal networks. *Proc. Natl. Acad. Sci.* **100** (2003) 7319–7324.

[39] P. Szmolyan and M. Wechselberger, Canards in $\mathbb{R}^3$. *J. Differ. Equ.* **177** (2001) 419–453.

[40] D. Terman, E. Lee, J. Rinzel and T. Bem, Stability of anti-phase and in-phase locking by electrical coupling but not fast inhibition alone. *SIAM J. Appl. Dyn. Syst.* **10** (2011) 1127–1153.

[41] S. Volkwein, Proper Orthogonal Decomposition: Theory and Reduced-Order Modelling. Lecture Notes. University of Konstanz, Konstanz (2013).

[42] D.D. Wang and A. Bordey, The astrocyte odyssey. *Progr. Neurobiol.* **86** (2008) 342–367.

[43] M. Wechselberger, Existence and bifurcation of canards in $\mathbb{R}^3$ in the case of a folded node. *SIAM J. Appl. Dyn. Syst.* **4** (2005) 101–139.

[44] C. Wyart, C. Ybert, L. Bourdieu, C. Herr, C. Prinz and D. Chatenay, Constrained synaptic connectivity in functional mammalian neuronal networks grown on patterned surfaces. *J. Neurosci. Methods* **117** (2002) 123–131.

[45] A.M. Zhabotinsky, H.G. Rotstein, I.R. Epstein and N. Kopell, A canard mechanism for localization in systems of globally coupled oscillators. *SIAM J. Appl. Math.* **63** (2003) 1998–2019.

Subscribe to Open (S2O)
A fair and sustainable open access model

<figure><img src="image_26.png" /></figure>

This journal is currently published in open access under a **Subscribe-to-Open model (S2O)**. S2O is a transformative model that aims to move subscription journals to open access. Open access is the free, immediate, online availability of research articles combined with the rights to use these articles fully in the digital environment. We are thankful to our subscribers and sponsors for making it possible to publish this journal in open access, free of charge for authors.

**Please help to maintain this journal in open access!**

Check that your library subscribes to the journal, or make a personal donation to the S2O programme, by contacting
**subscribers@edpsciences.org**

More information, including a list of sponsors and a financial transparency report, available at:
**https://www.edpsciences.org/en/maths-s2o-programme**