The emergence of polyglot entrainment responses to periodic inputs in vicinities of Hopf bifurcations in slow-fast systems

Entrainment refers to a type of synchronization in which a stable phase relationship is maintained between the output of an oscillator and an external periodic forcing signal.^{13,29,32,36,39,47,55,56,59,77} Entrainment patterns are described in terms of the number of input oscillations (N) that are phase-locked to a number of output oscillations (M), referred to as N:M patterns. The properties of entrainment are characterized in terms of Arnold tongues⁵⁵ in the input period-amplitude parameter space. Within each tongue, the system exhibits the same entrainment pattern. Arnold tongues have stereotypical shapes, from which their name derives, that form connected sets in the input period space for any fixed amplitude.^{1,17,21,33,37,60,65} However, we have observed the presence of 1:1 Arnold tongues with a complex “mountain/valley” shape in the well-known FitzHugh–Nagumo model^20,50 (and modified versions of this model) in the vicinity of Hopf bifurcations. These Arnold tongue diagrams appear to have a set of cojoined 1:1 Arnold tongues, and not all pairs of points with 1:1 entrainment can be joined by a (straight) line without that line going through a part of parameter space that is not in the 1:1 tongue. We refer to this situation as polyglot entrainment, and classical Arnold tongue diagrams with a single 1:1 region as monoglot entrainment. A similar polyglot phenomenon, although with different properties, has been observed in the Selkov model for glycolytic oscillations.⁷² 

Traditionally, entrainment studies focus on systems with intrinsic oscillators that exhibit sustained oscillations in the absence of time-dependent inputs. However, the entrainment patterns we observed include systems that exhibit damped oscillations in the absence of such inputs. This type of entrainment has received much less attention (but see Refs. 5, 24, 45, 74, and 78) and its associated dynamical properties have not been thoroughly investigated.

Although Arnold tongues are often used in entrainment studies, the dynamic mechanisms that give rise to the structure of Arnold tongues are not well understood. In particular, it is unclear how the tongues depend on the nonlinearities present in the system and the time constants at which the system operates. The goal of this paper is to investigate the dynamic mechanisms underlying the structure of disconnected 1:1 Arnold tongues in the FitzHugh–Nagumo model and other models of biological and chemical oscillators. Our results will shed light onto the synchronization properties of oscillatory systems in networks including both sustained oscillators and damped intrinsic oscillators.

We use the following periodically forced model of FitzHugh–Nagumo (FHN) type²⁰ used previously in Ref. 61,

where $v$ and $w$ are dimensionless variables representing a biological cell’s membrane potential and the recovery variable, respectively. The function $f(v)=−2v3+3v2$ is cubic with the minimum and maximum occurring at $(0,0)$ and $(1,1)$⁠, respectively. The parameters $ϵ$⁠, $α$⁠, and $λ$ represent the time scale separation between the two dependent variables, the slope of the $w$-nullcline (the curve defined by the set of points for which $dw/dt=0$⁠), and the displacement of the $w$-nullcline with respect to the $v$-nullcline (the curve defined by the set of points for which $dv/dt=0$⁠). The last term in Eq. (1) is a time-dependent, periodic input with a constant amplitude $A$⁠. We use two different types of waveforms for $F(t)$⁠: square wave and sinusoidal, with period $T$⁠, duty cycle 50%, and minima and maxima equal to $F=0$ and $F=1$⁠, respectively. We refer to the time intervals where $F>0.5$ as the “on” state and $F<0.5$ as the “off” state.

One advantage of using a square-wave forcing is that one can decompose the periodically driven system into two autonomous systems, one with the forcing turned off and another one with the forcing turned on. In the latter case, a change of variables $W=w−A$ allows us to move the constant forcing term $A$ into the second equation, yielding

When the forcing is on, the forcing amplitude contributes to the displacement of the $W$-nullcline. Alternatively, without the change of variables, $A$ causes the $v$-nullcline to displace abruptly between the two regimes. For sinusoidal inputs, this displacement is continuous and gradual.

System (1) and (2) undergoes a Hopf bifurcation^28,31,46,68 as the fixed-point crosses a vicinity of the minimum of the (cubic) $v$-nullcline (or alternatively, its maximum). The Hopf bifurcation point and the criticality properties can be computed in terms of the model parameters under the assumption that $ϵ≪1$⁠.^40,61,80 The fixed-point is stable (unstable) for $λ<λH$ (⁠$λ>λH$⁠) as shown in Fig. 1. The Hopf bifurcation is subcritical (supercritical) if $α<3$ (⁠$α>3$⁠). The dynamics of the FHN model for $a=2$ and various $λ$ values are shown in Fig. 2.

Arnold tongues are areas in parameter space where the system exhibits different locking modes in response to periodic forcing.^33,55 Outside the Arnold tongues, the input and output do not exhibit N:M phase-locking patterns.

Here, we focus on investigating 1:1 entrainment patterns. For the Arnold tongues we use here, the horizontal axis corresponds to the forcing period $T$ and the vertical axis corresponds to the forcing amplitude $A$⁠. Typically, with increasing forcing amplitude the Arnold tongue widens, indicating that entrainment can occur for a larger range of periods.^52,55

In order to numerically compute the Arnold tongues, for each fixed value of $A$⁠, we performed simulations (with initial conditions $V=0.05,w=0.05$⁠) to find the minimum and maximum values of $T$ for which the system exhibits 1:1 entrainment. We defined 1:1 entrainment to occur when the period of the system response (⁠$Ts$⁠) matched the period of the forcing (⁠$T$⁠) within some tolerance $η$⁠, i.e., $|T−Ts|<η$⁠, with $η=0.001$⁠. We computed $Ts$ as the time between peaks of the system response, using the MATLAB (The Mathworks, Natick, MA) findpeaks function with the option MinPeakProminence to set a value of one half the amplitude of the largest peak in each simulation as the threshold for admissible peaks.

To verify that the initial conditions used for the simulations are not impacting our results, we repeated the computation of the Arnold tongues shown in Fig. 3(A4) four times with initial conditions chosen at different locations along the stable limit cycle of the unforced system. In all cases, the recomputed Arnold tongues were identical to the ones shown in Fig. 3(A4). We then checked further for the possibility of multistability when the forcing amplitude is large by conducting simulations with 50 random initial conditions chosen from a uniform distribution (with bounds of [$−$1,1] for $V$ and $w$⁠). Specifically, with parameter values used for the monoglot Arnold tongue in Fig. 3(A4), we set $A=0.06$ and ran simulations over a range of forcing periods. For all 50 random initial conditions, we found the same 1:1 entrainment region as shown in Fig. 3(A4).

Phase-plane diagrams, time course plots, and Arnold tongues were computed numerically using the modified Euler method⁹ (a Runge–Kutta method of order 2) with a time step $Δt=0.05$ in MATLAB. Bifurcation diagrams were computed using the AUTO feature of XPPAUT.¹⁵ 

The concept of 1:1 polyglot entrainment refers to the ability of a system, in response to periodic inputs, to display separated 1:1 Arnold tongues in the input period-amplitude parameter space. In general, one expects the entrainment to be of monoglot type (for a given value of the input amplitude, the 1:1 entrainment segment is connected), which has been the focus of previous studies (but see Refs. 22 and 72).

In the following sections, we identify the conditions under which 1:1 polyglot entrainment emerges in a number of case studies with representative dynamic properties. We begin by using the forced FHN model [Eqs. (1) and (2)] and then we extend our results to a number of additional models forced by the same (square-wave) input and to the FHN model forced by sinusoidal inputs. Square-wave and sinusoidal inputs having the same number of cycles differ in the mechanism of transition between peaks and troughs, which are abrupt for square-wave inputs and gradual for sinusoidal inputs. It is not clear a priori whether and how this additional frequency content in square-wave input with respect to sinusoidal input affects the 1:1 entrainment properties of the forced systems.

The fixed-points for the (autonomous) FHN system are determined by the parameters $λ$ and $α$⁠. We use values of the parameter $ϵ$ for which there is a time scale separation between the activator (⁠$v$⁠) and inhibitor (⁠$w$⁠) variables and the FHN system exhibits relaxation oscillations. In the phase-plane diagram [e.g., Figs. 2(A) and 2(D)], the limit cycle trajectory transitions in between the two slow manifolds, which are located in vicinities of the left and right branches of the v-nullcline, along fast fibers.

To aid in the analysis, we consider autonomous FHN models with a constant input equal to the periodic input’s amplitude. The fixed-points of the time-independent forced system are different and may have different stability properties from the unforced system. We note that for the FHN model, the addition of a constant forcing to the activator equations is equivalent to an increase in the parameter $λ$⁠.

In the following sections, we are going to restrict our focus to the cases of $α=2$⁠, for which we have a subcritical Hopf bifurcation upon varying $λ$⁠, and $α=4$⁠, for which we have a supercritical Hopf bifurcation upon varying $λ$⁠. We analyze the entrainment results for different $λ$ values. We split the $λ$ axis into regions where we have monoglot (single 1:1 region in the Arnold tongue diagram) or polyglot (multiple 1:1 regions in the Arnold tongue diagram) entrainment.

We consider periodic forcing with different frequencies and two types of waveforms, sinusoidal and square-wave, that are representative of gradual and abrupt transitions in the external forcing signal. Both are amenable to analysis using dynamical systems tools. In the latter case, the dynamics of the periodically forced system can be decomposed into two two-dimensional subsystems with constant forcing, which simplifies the analysis. We call the times at which the forcing turns on and off the switching times, and the corresponding points in the phase-plane diagrams the switching points. Switching points serve as initial conditions at the transitions between the on and off states of the forcing (once 1:1 entrainment occurs, the switching points always stay fixed).

Here, we consider the FHN model with $α=2$ and $ϵ=0.01$ and representative values of $λ$ using a square-wave input with amplitude $A$ and a 50% duty cycle. Figure 2 shows representative phase-plane diagrams (arranged in the order that will be used in the following sections). As $λ$ increases (panels B1, C1, D1, A1) the $w$-nullcline moves to the right and the fixed-point transitions from a stable node (away from the knee of the $v$-nullcline, on the left, not shown) to a stable focus [Figs. 2(B1) and 2(C1)], to an unstable focus [Figs. 2(D1) and 2(A1)], to an unstable node (away from the knee of the $v$-nullcline, on the right, not shown). The stable focus in Fig. 2(B1) is located close to the boundary between stable foci and nodes [further away from the knee of the $v$-nullcline than in Fig. 2(A1)] and, therefore, it exhibits strongly damped oscillations, but not sustained oscillations.

The effect of periodic forcing by positive square wave inputs (with amplitude $A$⁠) can be thought of as an abrupt transition between two autonomous FHN systems, one with $A=0$ (forcing “turned off”) and the other with a nonzero value of $A$ (forcing “turned on”). We refer to the latter as the tonically forced FHN system. In the absence of time-dependent forcing, increasing (constant) values of $A$ have the same effect as increasing values of $λ$ (increasing the constant forcing $A$ is equivalent to increasing $λ$ via a translation of the variable $w$⁠). Therefore, as $A$ increases, the fixed-point moves to the right and its stability properties change accordingly. If the unforced FHN model exhibits large amplitude oscillations (LAOs) for $A=0$⁠, then it will do so for $A>0$ [e.g., Figs. 2(A2) and 2(D2)]. If the unforced FHN exhibits small amplitude (damped) oscillations (SAOs), then the tonically forced FHN model may exhibit SAOs with a weaker damping [e.g., the transition from Figs. 2(B1) and 2(C1)] or LAOs [e.g., Figs. 2(B2) and 2(C2)]. We consider the response of the FHN model to periodic forcing in these three scenarios below. Finally, if the unforced FHN model has a stable node, then the tonically forced FHN model may have a stable node, exhibit SAOs, or exhibit LAOs. The results for this scenario is similar to the previous one and are presented in the supplementary material. We note that the existence of LAOs does not preclude the existence of damped SAOs for the small bistability range of values of $λ$⁠.

Here, we describe the monoglot entrainment properties of the FHN system when both the unforced (⁠$A=0$⁠) and the constantly forced system have an unstable focus with large amplitude oscillations (LAOs) with a natural (intrinsic) period (⁠$T0≈200$ for the unforced system). The trajectories of both the unforced and constantly forced system exhibit SAOs as they spiral out toward the LAO limit cycle [Fig. 2(A)].

Figure 3(A) shows representative entrained responses of the FHN model to periodic stimulation with the same amplitude as in Fig. 2(A2). The 1:1 entrainment reflects the similarity between periods of the forcing and the unforced system as the system dynamically transitions between the regimes described in Fig. 2(A). The region of 1:1 entrainment in the Arnold tongue [Fig. 3(A4), green] widens with increasing values of $A$⁠. Figures 3(A1) and 3(A3) show examples of entrained responses with different input–output ratios illustrating the ability of the FHN model to follow the corresponding inputs via a mixed-mode oscillations (MMOs) mechanism.

Figure 2(B) illustrates the behavior of the FHN model when the unforced system (⁠$A=0$⁠) is a stable focus with strongly damped SAOs [Fig. 2(B1)] and the constantly forced system has an unstable focus [as in Fig. 2(A2)] and displays LAOs [Fig. 2(B2)].

The entrainment of non-sustained oscillators has been significantly less studied than the entrainment of sustained oscillators. One important conceptual difference between the two protocols is the lack of a reference period in the cell to be entrained. In the previous section, we argued that the entrainment can be understood as the abrupt transition between two oscillatory regimes and entrainment occurs when the time scales of the two regimes are compatible. Here, there is no oscillatory time scale associated with the unforced system. Furthermore, while the constantly forced system is able to show LAOs for large enough values of $A$ [Fig. 2(B2)], for lower values of $A$⁠, the dynamics remains in a damped oscillatory regime without LAOs (not shown). Therefore, entrainment, when it occurs, cannot be explained simply by the compatibility of time scales, but requires a deeper explanation in terms of dynamical systems ideas.

Figures 3(B4) and 3(C4) show the 1:1 entrainment regions (green) for these two cases. Notably, they have a different shape than the standard Arnold tongues [e.g., Fig. 3(A4)]. In particular, 1:1 entrainment is broken as $A$ increases for constant values of $T$⁠, but maintained as $T$ increases for constant values of $A$ below some critical value (dashed horizontal line). This critical value corresponds to the transition between two stability regimes for the constantly forced system: unstable focus with LAOs [above, Fig. 3(B4)] and stable focus [below, Fig. 3(C4)].

Comparison between Figs. 3(B2) and 3(C2) shows relatively similar patterns for the entrainment in the two regimes for the same oscillatory input. Comparison between Figs. 3(B3) and 3(C3) shows that the increased entrainment robustness as $T$ increases in Fig. 3(C3) is associated with the presence of MMOs in the forced system reflecting the dynamic interaction between the forcing and the cell.

The results for this scenario are similar to the results for a stable focus that is not in the vicinity of a Hopf bifurcation and are presented in Fig. S11.

The 1:1 entrainment for the FHN model in the parameter regimes discussed above breaks in different ways depending on whether $T$ is larger or smaller than its values in the 1:1 entrainment region, and whether the unforced system is in a LAO regime or stable focus regime (Fig. 3). In all cases, the N:M patterns have $N>M$ for smaller values of $T$ and $N<M$ for larger values of $T$⁠. In the first case, the breaking of 1:1 entrainment is due to a cycle skipping mechanism where the input turns off before the response (⁠$v$⁠) jumps up, and therefore the LAO fails to be generated. In the second case, the breaking of 1:1 entrainment is due to a cycle adding mechanism. In Fig. 3(A3), in contrast to the 1:1 entrainment case [e.g., Fig. 3(A2)], $v$ succeeds to jump up in spite of the fact that the input turns off. This seeming disconnect between the input and the response leads to cycle adding and indicates a more complex interaction of effective time scales. In Fig. 3(B3), the cycle adding mechanism results from the ability of the input to produce two output cycles per input cycle while it is on. Note that the inter-oscillation interval (IOI) is larger when the input is off than when it is on. None of these effects break 1:1 entrainment for larger values of $T$ in Fig. 3(C). In fact, Fig. 3(C3) shows that 1:1 entrainment is maintained for values of $T$ for which it is broken in the panels described above. This result is due to $v$ remaining roughly constant when the input turns off, instead of $v$ increasing as in the other panels.

Here, we use dynamical systems tools (extended phase-plane analysis) to provide a more detailed explanation of the results discussed above. As we mentioned above, the dynamics of the forced system can be interpreted as the abrupt transition between two autonomous FHN systems, the off (⁠$A=0$⁠) and on (⁠$A>0$⁠) states. Dynamically, the evolution of the response trajectory is governed by an abrupt alternation between the phase-plane diagrams for the corresponding values of $A$⁠. Figure 4 shows the superimposed phase-plane diagrams for $A=0$ (black) and $A>0$ (red). We use “numbered arrows” (in black and red) to denote the switching points between the unforced and forced systems.

Roughly speaking, the entrainability of the FHN system requires compatibility between its effective time scale and the time scale of the input. The former is primarily determined by the slow manifolds located in vicinities of the $v$-nullcline. In the unforced system (⁠$A=0$⁠), the limit cycle trajectory evolves along these manifolds and jumps up and down toward the right and left branches of the $v$-nullcline, respectively, as time progresses [e.g., Fig. 2(A)]. The interaction between the two time scales is interpreted as the $v$-nullcline rising and shifting down as the input turns on and off, respectively. This interaction process is affected by other factors, including the presence of fixed-points in one of the two regimes [e.g., Figs. 2(B2) and 2(B3)] that may transiently create slower dynamics. In addition, the presence of nonlinearities may disrupt the evolution of the trajectory as the $v$-nullcline moves and adds additional time scales by creating small amplitude oscillations (SAOs), for example. These may favor or oppose the entrainability of the system. While oscillatory time scales are usually associated with sustained oscillators, they are still present in damped FHN oscillators (close, but away from the sustained oscillations regime) and they become apparent (and functional) as they interact with the input. The presence of oscillatory time scales may also favor or oppose the entrainability of the FHN system.

Below, we discuss a number of representative cases in detail, both to highlight the basic aspects of entrainability in the FHN system and to develop language for subsequent sections.

When the forcing turns on (2), the $v$-nullcline raises and the trajectory, located further away from the lower knee, jumps up toward the right branch, moves up (slowly) along it, reaches the upper knee, and jumps down toward the left branch. The trajectory then moves down (slowly) along the left branch until the forcing turns off (1). When this happens, the $v$-nullcline shifts down, affecting the direction of motion of the trajectory. In Fig. 4(A2), the trajectory is about halfway to the lower knee and is affected only slightly by the shift in the $v$-nullcline. In Fig. 4(B2), the trajectory is affected significantly by the shift as it jumps to the left, and then moves down toward the lower knee at which point the forcing turns on (2). In Fig. 4(C2), the size of the effect that the shift in the nullcline has on the trajectory is in between the two cases described above.

The dynamics in Fig. 4(C3) are similar to those described above, except that the on-off transition of the forcing (1) occurs very close to the lower knee. The trajectory jumps to the left near the lower knee. Because of that, in contrast to the previous figures, the trajectory evolves on a much slower time scale than the ones at higher points on the slow manifold (in vicinities of the left branch). Therefore, it does not reach the lower knee, while the forcing is off and only jumps up when the forcing turns on (2).

When the forcing turns off (1), the $v$-nullcline shifts down and the trajectory is left closer to the upper knee, thus accelerating the jump down toward the left branch. The trajectory then evolves slowly along the left branch until it reaches the lower knee, at which point the forcing turns on (2). When this happens, the $v$-nullcline raises and trajectory jumps up to the right branch. The trajectory evolves along the right branch, reaches the upper knee, and jumps down. Because $T$ is relatively large (forcing relatively slow), the trajectory is able to jump up once more before the forcing turns off (1), and, therefore, can produce two LAOs per input cycle.

The (very) slow time scale operating in Fig. 4(C3), where 1:1 entrainment is maintained, is due to the fact that the tonically forced system is in a stable focus regime. The switch to an unstable focus regime (for the tonically forced system) in Fig. 4(B3) “destroys” this time scale and therefore the maintenance of 1:1 entrainment is no longer possible. More specifically, the increase in the input amplitude (⁠$A$⁠) as the tonically forced system transitions from having a stable fixed-point [Fig. 3(C3)] to an unstable fixed-point [Fig. 3(B3)], allows the forced oscillator corresponding to the latter to jump up just before the forcing is turned off. In contrast, for the lower value of $A$ [Fig. 3(C3)], the forced oscillator is still in a vicinity of the left branch of the $v$-nullcline when the oscillator is turned off and remains there until the forcing is turned on again. The resulting longer time scale is associated with the 1:1 entrainment observed for this value of $A$⁠.

These dynamics are more complex than those for Fig. 4(B3) since the first and third oscillations in the shaded region in Fig. 3(A3) occur while the forcing is off. Starting from the beginning of the shaded region, when the forcing turns off (1) the $v$-nullcline shifts down, the trajectory jumps down, it moves slowly along the lower knee and jumps up toward the right branch creating an oscillation, which is induced by the forcing ceasing to be active rather than by the forcing turning on. A second oscillation is generated after the forcing turns on (2) and this oscillation remains after the forcing turns off (3) since the trajectory needs to reach the upper knee in order for the oscillation to be terminated. In other words, the cell’s intrinsic dynamics dominate over the forcing here since the forcing does not induce or terminate this oscillation. The third oscillation is induced by the forcing (4). When the forcing turns off (5), the trajectory is very close to the lower knee, and therefore the trajectory jumps up well in advance of the forcing turning on (6). In fact, the trajectory is able to reach the upper knee when the forcing turns on and jumps down roughly when this happens.

The 2:1 patterns in both figures are generated by a cycle skipping mechanism. When the forcing turns on (2), the trajectory is approaching the stable focus. The $v$-nullcline raises and releases the trajectory that jumps up toward the right branch. The turning off of the forcing (3) roughly coincides with the passing of the trajectory to the upper knee. The trajectory jumps down toward the left branch (it would do that even in the absence of the forcing remained on, but the shifting down of the $v$-nullcline accelerates the process). The forcing turns on again (4) when the trajectory is evolving along the left branch. The ability of the cell in panel B1 to generate an oscillation in response to this depends on the competition between two time scales determined by the slow manifold and the input period. Cycle skipping results because the forcing turns off (1) before the trajectory manages to pass through the lower knee and jump up. The trajectory then needs to wait for the forcing to turn on (2) in order to jump up. In the absence of the forcing, the trajectory would converge to the stable focus. The dynamics in panel C1 is similar to those in panel B1, except that the forcing turning off (1, currently 3) does not prevent the trajectory to jump up. Even if the forcing remained on the trajectory would not jump up, but converge to the stable focus. From that perspective, the successive forcing turning off (1, currently 3) and on (2, currently 4) are necessary to maintain the 2:1 pattern.

The mechanism of generation of 4:3 patterns is different than the ones for the 2:1 patterns described above, resulting from the fact that both the unforced and tonically forced systems are oscillators and have intrinsic oscillatory time scales that operate independently of the input. The trajectory jumps up toward the right branch even if the forcing is turned off (1) and is already moving along the left branch when the forcing turns off (3). Moreover, the forcing turning on (2) has little effect on the evolution of the trajectory, which is already moving along the right branch. The stronger effect the input has on the response oscillatory dynamics, and the one that determines the 4:3 pattern, is the forcing turning off (1) that delays the jumping up of the trajectory.

Both the standard and non-standard 1:1 Arnold tongues discussed in the previous section consisted of 1:1 entrainment regions where horizontal segments of 1:1 entrainment for constant values of $A$ (as well as vertical segments for constant values of $T$⁠) are connected. Here we discuss a more complex type of 1:1 entrainment region structure arising in the FHN model, consisting of disconnected horizontal segments for low values of $A$⁠, while the segments remain connected for higher values of $A$ [e.g., Figs. 5(A1) and 5(B1)]. These regions have a multi-lobed 1:1 Arnold tongue. We refer to this phenomenon where 1:1 entrainment is interrupted as $T$ changes (for constant values of $A$⁠), and is restored as $T$ changes further, as polyglot entrainment. We use the same values of $α=2$ and $ϵ=0.01$ as in the previous section as well as square-wave input with amplitude $A$ (50% duty cycle).

We have observed polyglot entrainment in two parameter regimes: when the unforced system has a stable focus and is a damped oscillator [Fig. 5(A1)], or when it has an unstable focus and is a sustained oscillator [Fig. 5(B1)]. We discuss these two cases below. The main differences between them and the ones described in the previous sections (monoglot entrainment) is the location of the fixed-point, which is closer to the minimum of the $v$-nullcline in Fig. 2(C1) than in Fig. 2(B1). This creates a region of sensitivity to perturbations that is responsible for the generation of the multiple 1:1 tongues.

We observed polyglot entrainment responses when the unforced FHN system has a stable focus in the vicinity of a Hopf bifurcation located near the lower knee of the cubic nullcline (Fig. 5). Here, we take two representative cases. First, when the forced system has unstable focus (along with sustained oscillations) as shown in Fig. 5(A1). Figures 5(A3), 5(A5), and 5(A7) show representative traces for the 1:1 entrainment in the three tongues shown in Fig. 5(A1) (above the horizontal line). They are generated by an MMO mechanism that controls the time scale by adding SAOs as the period increases and therefore the oscillatory output adapts to the input. The corresponding extended phase-plane diagrams are presented in Fig. 6.

For the first (left) tongue [Figs. 5(A3) and 6(A)], there are no MMOs and the entrainment is standard. The trajectory arrives at the lower knee of the $v$-nullcline and jumps up toward the right branch. The motion of the $v$-nullcline following the input does not interfere with this process.

For the second (middle) tongue [Figs. 5(A5) and 6(C)], the trajectory crosses the $v$-nullcline in a vicinity of the lower knee, to the right, and turns back instead of jumping up toward the right branch. The forcing turns on when the trajectory is below the lower knee. This raises the $v$-nullcline leaving the trajectory in a region of fast fibers and therefore the trajectory jumps up. For the third (right) tongue [Fig. 5(A7)], the mechanism is similar, but the larger forcing period allows the trajectory to evolve around the lower knee twice before jumping up.

These mechanisms of 1:1 entrainment can be disrupted in different ways [Figs. 5(A2), 5(A4), 5(A6), and 5(A8)]. In Fig. 5(A2), the forcing is too fast, and activates the SAO mechanism only once every three cycles. In Fig. 5(A4) [see Fig. 6(B)], when the forcing turns on (2), the $v$-nullcline raises and the fixed-point raises too. As a result, the trajectory spirals out around this (forced) unstable fixed-point until the forcing turns off (3). When this happens, the trajectory spirals into the (unforced) stable fixed-point until the forcing turns on again (4) and the trajectory is able to jump up toward the right branch as the result of the $v$-nullcline raising and leaving the trajectory in a region of fast fibers (below the knee). These additional oscillations disrupt the ability of the cell to follow the input. The disruption mechanisms in Figs. 5(A6) and 5(A8) are similar to this one.

Second, when the forced system has a stable focus (and damped oscillations) as shown in Fig. 5(B1). Figures 5(B3), 5(B5), and 5(B7) show representative traces for the 1:1 entrainment in the three tongues shown in Fig. 5(B1) (below the horizontal line). They are also generated by MMO mechanisms similar to the ones described above as illustrated in the Figs. 7(A) and 7(C). The breaking of entrainment mechanisms is also similar to those discussed above, except that in some cases, the timing of the output is not good enough to cause the cell to jump up toward the right branch and produce oscillations [Figs. 5(B2) and 5(B6)]. By the time the forcing is turned on, the trajectory is “inside” the $v$-nullcline, and by the time the trajectory arrives in the region of fast fibers (below the $v$-nullcline), the forcing is off.

When the unforced system has an unstable focus, the forced system also has an unstable focus. Polyglot entrainment [Fig. 8(A1)] has also been observed when the unforced system has an unstable focus (and hence self-sustained oscillations) in the vicinity of a Hopf bifurcation. Figures 8(A3) and 8(A5) show representative traces for the 1:1 entrainment in the two tongues shown in Fig. 8(A1). In the polyglot entrainment observed here, the tongues do not merge as the forcing amplitude increases. The corresponding extended phase-plane diagrams are presented in Fig. 9.

For the first (left) tongue [Figs. 8(A1) and 8(A)], there are no MMOs and the 1:1 entrainment is standard. The trajectory when reaches the lower knee of the $v$-nullcline jumps up toward the right branch. For the second (right) tongue, 1:1 entrainmemt arise as a result of MMOs. When the forcing turns off, the trajectory is located in the very close vicinity of an unstable fixed-point (focus), thus, larger forcing period allows the trajectory to spiral out around the lower knee. When the forcing turns on, the trajectory is located far away from the vicinity of an unstable fixed-point, thus, jumps up toward the right branch and the large forcing period allow the trajectory to traverse around the phase-plane, and completes the round [Fig. 9(C)].

The disruption of 1:1 entrainment is shown in Fig. 8(A4) in which 2:2 patterns are observed for the corresponding forcing period. When the forcing turns off for the first forcing cycle, the trajectory is close to the unstable fixed-point near the lower knee, hence, it spirals out and when the forcing turns on, the trajectory is far from fixed-point of the forced system, hence the trajectory shoots up to the right creating a spike. For the second input, the forcing turns off and on far away from the lower knee thus a spike is observed with no SAOs [Fig. 9(B)].

In the above sections, we have focused on a scenario in which the unforced system has some degree of time scale separation (⁠$ϵ=0.01$⁠), undergoes a subcritical Hopf bifurcation as $λ$ is varied (⁠$α=2$⁠), and the periodic forcing $F(t)$ takes the form of a square wave with a 50% duty cycle. Since the location and criticality of the Hopf bifurcation depends on model parameters, we also explored other parameter regimes. We have observed both monoglot and polyglot entrainment in systems with less time scale separation (⁠$ϵ=0.1$⁠), a supercritical Hopf bifurcation as $λ$ is varied (⁠$α=4$⁠), square-wave forcing with various duty cycles, and sinusoidal forcing. Here, we briefly discuss the results in these other parameter regimes.

First, we consider time scale separation. As $ϵ$ increases, the separation of time scales becomes less pronounced and trajectories that would leave the lower knee region along fast fibers may no longer do so. Nevertheless, for $ϵ=0.1$⁠, we still observe both standard and nonstandard monoglot Arnold tongues similar to the case with $ϵ=0.01$ (cf. Fig. 3 and Fig. S16). For $ϵ=0.1$⁠, we also observe polyglot entrainment responses (cf. Fig. 5 and Fig. S17). However, for $ϵ=1$⁠, the forced system only exhibits subthreshold responses (moving around the lower knee) and thus we do not consider them to be entrained according to our previous description of entrainment. In these cases, the stable fixed-points are located along the middle branch of the cubic nullcline.

Second, we consider the criticality of the Hopf bifurcation. When the unforced system undergoes a supercritical Hopf bifurcation (⁠$α=4$⁠), we see a monoglot entrainment structure with standard and nonstandard Arnold tongues similar to the subcritical Hopf case (⁠$α=2$⁠, cf. Figs. 3 and S1). For $α=4$⁠, we also observed polyglot entrainment responses (cf. Figs. 5 and S2, see also Figs. S3–S5). The difference we observed is that the region of polyglot entrainment for $α=4$ is larger than for $α=2$ (cf. Figs. 1 and S6).

Third, we consider the effect that the duty cycle of the square-wave forcing has on monoglot and polyglot entrainment. We checked cases where the duration of the “on” period is either longer or shorter than the duration of the “off” period (75% on, 25% off and 25% on, 75% off, respectively). We found that varying the duty cycle had very little effect on the Arnold tongue in the standard monoglot entrainment case [cf. Figs. 3(A1) and S18]. However, in the case of nonstandard monoglot entrainment, reducing the relative duration of the “on” period (a 25% duty cycle) leads to a wider Arnold tongue [cf. Figs. 3(B1) and S18]. Similarly, in the case of polyglot entrainment, the size of the 1:1 entrainment regions in the input period-amplitude plane increases as the duty cycle is decreased [Figs. 5(A1) and S19].

Finally, we consider the form of the periodic forcing. With sinusoidal forcing, we observe standard monoglot Arnold tongues but we do not find the type of nonstandard monoglot Arnold tongues that were observed with square-wave forcing (cf. Figs. 3 and S12). In addition to monoglot entrainment, with sinusoidal forcing we also observed polyglot entrainment (Fig. S13). We note that the range of $λ$ values exhibiting polyglot entrainment with sinusoidal forcing is larger than the corresponding range for square-wave forcing (cf. Figs. 1 and S15). To explain this observation, we focused on a case (⁠$λ=−0.065$⁠, $α=2$⁠, $ϵ=0.01$⁠) for which the unforced system has a stable focus and we have monoglot entrainment with square-wave forcing but polyglot entrainment with sinusoidal forcing (see Fig. S25). With $A=0.03$ and a forcing period of $T=230$⁠, we have 1:1 entrainment with both types of forcing. In both situations, when the forcing is either off (in the square-wave case) or gradually approaching zero (in the sinusoidal case), then the trajectory is moving in the direction of converging to the stable fixed-point. With sinusoidal forcing, once the forcing reaches zero then it begins to increase and the $v$-nullcline gradually moves up. The trajectory stays in the vicinity of the moving fixed-point, and depending on the speed at which the trajectory is moving relative to the speed at which the $v$-nullcline is moving, the trajectory may cross the nullcline and change direction. When this occurs, it can disrupt 1:1 entrainment (see Fig. S26 and Supplemental Movie 1). As the forcing period increases, the $v$-nullcline moves more slowly, affecting the relative speeds of the trajectory and nullcline movements. Thus, with sinusoidal forcing, 1:1 entrainment can be lost and then regained depending on these relative speeds as the forcing period is varied, leading to a polyglot entrainment structure. In contrast, with square-wave forcing, when the forcing turns on the $v$-nullcline abruptly moves away from the trajectory located near the stable fixed point of the unforced system (see Fig. S27 and Supplemental Movie 2). Thus, the trajectory does not have an opportunity to cross the $v$-nullcline and change direction. This leads to a monoglot entrainment structure as the forcing period is increased.

In the previous sections, we analyzed the phenomenon of polyglot entrainment in the FHN model and uncovered the underlying mechanisms. The presence of polyglot entrainment (without using this name) has been reported before in the Selkov model.⁷² The question arises whether this type of entrainment is a general phenomenon for a certain class of oscillatory systems. Here, we briefly explore the 1:1 response patterns in models of neuronal, circadian, chemical and glycolytic oscillations in vicinities of a Hopf bifurcation and compare them with a representative example of another bifurcation type. In contrast to the FHN model, in these additional models, the parameters have a direct biological or chemical interpretation.

We investigated polyglot entrainment responses in two parameter regimes corresponding to excitability properties being governed by a Hopf bifurcation (type II) or a SNIC (saddle-node on an invariant cycle) bifurcation (type I) (Fig. S21).

In type I excitability, a SNIC bifurcation occurs at an $Iapp$ value near the lower knee of the $v$-nullcline and is responsible for the onset of oscillations. A Hopf bifurcation occurs at a higher $Iapp$ value, near the upper knee of the $v$-nullcline, and is responsible for the termination of the oscillations. In the vicinity of the Hopf bifurcation, we found polyglot entrainment responses. We did not observe polyglot entrainment near the SNIC bifurcation.

For type II excitability, we have Hopf bifurcations near the lower and upper knees that are responsible for onset and termination of oscillations, respectively. We found polyglot entrainment in the vicinity of both Hopf bifurcations [Figs. 10(A) and S22].

In both the 2D Oregonator and Novak–Tyson models, oscillations are generated via a Hopf bifurcation when the fixed-point is located in the vicinity of a knee of a cubic-like nullcline. We observed polyglot entrainment in both cases [see Figs. 10(B) and S23 for oregonator; Figs. 10(C) and S24 for Novak–Tyson].

Our results above suggest that the polyglot entrainment is generally found for models where the generation of oscillations is governed by a Hopf bifurcation. Here, we show that polyglot entrainment is not always found in such models. More specifically, the Lengyel–Epstein model^42,43 for chemical oscillations is a counterexample to the putative principle that models having a Hopf bifurcation as the mechanism for generation of oscillations exhibit polyglot entrainment [Fig. 11(A4)].

However, in contrast to the previous examples, the Hopf bifurcation fixed-point in the Lengyel-Epstein system is located far from the knee of the cubic-like nullcline [Figs. 11(A1) and 11(A2)].

We have identified and explained the dynamical mechanisms underlying a novel type of entrainment when the unforced system, not necessarily a sustained oscillator, is in the vicinity of a Hopf bifurcation. The 1:1 entrainment area in the frequency-amplitude parameter space has a different shape than that of a classical Arnold tongue diagram and has a complex boundary with multiple regions of 1:1 entrained behavior. We termed this polyglot entrainment and referred to the type of entrainment characterized by the classical Arnold tongue structure with a single 1:1 region as monoglot entrainment. Split 1:1 entrainment has been observed in models⁷² and in experiments with a saline oscillator²² in response to periodic inputs. To our knowledge, this phenomenon has not been fully characterized and the mechanisms underlying polyglot entrainment are not well understood. In this paper, we set out to address these issues.

The unforced (or constantly forced) FHN model is two-dimensional and amenable to phase-plane analysis. When the FHN model is subjected to time-dependent forcing, the system becomes three-dimensional. However, if the time-dependent forcing is periodic with abrupt transitions between on and off states, such as square-wave forcing, then we can decompose the system into a pair of two-dimensional systems (one with the forcing on, and one with the forcing off) and use phase-plane analysis to study how changes in the vector field in response to the forcing turning on and turning off affect the evolution of the trajectory leading to 1:1 entrainment or the lack thereof. The changes in the $v$-nullcline and in the location and stability of the fixed-points in response to the forcing parameters aids in our analysis. Using this extended phase-plane analysis approach, we characterized how 1:1 entrainment can be broken and then regained as the forcing period is varied, resulting in polyglot entrainment. This type of entrainment is only possible when the unforced system is in the vicinity of a Hopf bifurcation and exhibits either weakly damped oscillations (near the Hopf point with $λ<λH$⁠) or weak self-sustained oscillations (near the Hopf point with $λ>λH$⁠), and the forcing is weak enough. For stronger forcing, the 1:1 tongues merge and the polyglot entrainment transitions into monoglot entrainment (e.g., Fig. 5).

1:1 entrainment reflects the compatible interaction between two time scales, the input period (or, alternatively, the duration of the on phase) and the time it takes the trajectory to evolve between two upstroke times. In terms of the phase-plane diagram, the latter is the time it takes the trajectory to evolve from the upstroke point [e.g., the red dot in Figs. 6(A3) and 6(C3)] back to it. 1:1 entrainment occurs if and only if these two characteristic times match. The forcing activation raises the $v$-nullcline and facilitates the upstroke by leaving the trajectory below the $v$-nullcline and away from the basin of attraction of the stable fixed-point, if it exists. The stronger the input amplitude $A$⁠, the larger the distance between the trajectory and the $v$-nullcline, and the easier it is for the cell to produce an upstroke.

In the simplest 1:1 entrainment situation, the cellular response time scale is controlled by the slow manifolds. For low enough values of $ϵ$⁠, the trajectory evolves in vicinities of the right and left branches of the $v$-nullcline, and the upstroke and downstroke times are negligible. This type of 1:1 entrainment can be lost for various reasons as the period $T$ increases. One reason is that the trajectory may arrive to the upstroke region in the phase-plane diagram too early, and cross the unstable branch of the $v$-nullcline engaging in small amplitude oscillations [e.g., Fig. 6(A3), in Figs. 6(B3) and 6(C3) the black trajectories cross the unstable, black, branches of the corresponding $v$-nullclines]. In this case, the rise of the $v$-nullcline due to the forcing may not be enough to leave the trajectory below the $v$-nullcline, and therefore the upstroke is missed. Another reason is that trajectories arriving too early to the upstroke region may produce an upstroke while the forcing is still off.

Polyglot entrainment occurs if the 1:1 entrainment can be recovered as $T$ continues to increase. This requires an additional component to the cellular time scale that delays the arrival of the trajectory to the upstroke region until the forcing is turned on. Our analysis demonstrates that this additional component is provided by the time it takes the trajectory to move around the lower knee of the $v$-nullcline roughly an integer number of times, and be back in the upstroke region when the forcing is turned on [e.g., Fig. 6(C3)]. This is reflected as small amplitude oscillations in the cell’s responses during the forcing’s off phase.

Key to these mechanisms is the development of a slow time scale by trajectories moving around the lower knee of the $v$-nullcline, which allows for the cell’s response to match the forcing period. This mechanism can operate for systems in which excitability occurs via Hopf bifurcations, but not saddle-node on invariant circle bifurcations. These mechanisms also require that $ϵ$ is small enough so that the system produces relaxation oscillations. For larger values of $ϵ$⁠, trajectories are able to cross the middle branch of $v$-nullcline further away from the knee at the expense of the upstroke. By moving the location of the fixed-point toward the center of this middle branch, the cell recovers the ability to produce an upstroke, but we have not observed polyglot entrainment for the same parameter settings that produce polyglot entrainment with smaller values of $ϵ$⁠. In this paper we focused on polyglot entrainment in systems exhibiting oscillations of relaxation type. Further research is needed to investigate whether and under what conditions polyglot entrainment can be obtained for other types of oscillators.

Additionally, these mechanisms also require that $A$ is small enough. For larger values of $A$⁠, the $v$-nullcline raises to higher levels in response to the forcing turning on, thus facilitating the upstroke for larger ranges of $T$⁠. Therefore, the mechanisms of disruption of 1:1 entrainment discussed above cease to operate.

The literature on nonlinear dynamics and entrainment is vast with many classical texts and a wide variety of applications to natural and engineered systems.⁵⁵ Two studies that report results closely related to polyglot entrainment are González et al.²² and Verveyko et al.⁷² In Ref. 22, the authors perform experiments with a saline oscillator. This simple hydrodynamical oscillator consists of a container with salt water submerged in a larger container of fresh water. A tiny hole in the container with salt water (which is plugged until the start of the experiment) allows salt water to flow downward out of the smaller container due to the higher density of salt water compared to fresh water. However, after a few minutes, the flow reverses and fresh water flows upward into the smaller container for several tens of seconds.

This cycle then repeats thousands of times until the saline gradient dissipates. The oscillator can be subjecting to time-dependent forcing by periodically infusing a fixed volume of fresh water and then withdrawing it using a syringe, referred to as a volume pulse protocol. González et al. found that when the period of the volume pulse forcing was increased from 15 to 20, 30, 50, and 60% of the intrinsic period of the saline oscillator, the observed entrainment patterns went from 1:1 to 2:2, 2:1, 2:2, and back to 1:1. Although an Arnold tongue diagram was not plotted for these experimental results, we anticipate such a diagram would show a polyglot entrainment structure similar to what we have found in the FHN and other models.

In Ref. 72, the Selkov model of glycolytic oscillations was studied under periodic forcing. Glycolysis is a metabolic process that converts glucose into pyruvic acid, thereby releasing energy that is used to form ATP. The Selkov model⁶⁶ is a two-dimensional ODE system describing a key step of the glycolytic reaction. Verveyko et al. modeled periodic forcing of a model parameter to represent oscillations in the substrate for this enzymatic reaction. They characterized the entrainment patterns using an Arnold tongue diagram (see Fig. 2 of their paper) and found that inside the 1:1 tongue there is an area where trajectories diverge to infinity (see also Ref. 8) and therefore are not entrained. This leads to a split 1:1 tongue, with interruptions between regions of 1:1 entrainment, similar to our Arnold tongue diagrams for polyglot entrainment. However, the mechanism for polyglot entrainment is different, as the breaking of 1:1 entrainment is due to the existence of other $N:M$ patterns rather than unbounded solutions.

Although the entrainment of damped oscillators has been less explored, there have been some experimental and theoretical studies on this topic. Experimental work shows the presence of damped oscillations in the circadian clock of insects.⁴ In mammals, the central circadian clock (the suprachiasmatic nucleus, or SCN) is comprised of a heterogeneous cell population with a large proportion of cells that show weak (damped) oscillations.⁷³ In a theoretical work, Woller et al. studied a classic circadian oscillator model (the Goodwin model) and compared the entrainment properties when the unforced system shows either damped oscillations or limit cycle oscillations.⁷⁸ They found that the range of entrainment is wider (entrainment is achieved over a wider range of forcing periods) when the unforced system exhibits damped oscillations. With periodic forcing applied to limit cycle oscillators, they found a richer set of behaviors, including non-entrained chaotic and quasiperiodic solutions, which are not desirable biologically. A modeling study by Westermark et al. investigated the effect of noise on damped oscillators.⁷⁶ They showed that noise can generate sustained oscillations in a damped oscillator. They compared their simulations to experimental data from individual cells, and concluded that whether the circadian clock is a damped or self-sustained oscillator at the single-cell level could not be determined.

Further work by Gonze et al.²⁴ considered synchronization of populations of coupled circadian oscillators. In their work, they showed that when the coupled oscillators being considered are damped, high synchronization efficiency is achieved. Work by Bernard et al.⁵ also shows that efficient and robust synchronization is achieved in coupled damped oscillators. Komin et al.³⁸ studied entrainment of coupled oscillators where each oscillator has a different intrinsic period. They found that cells having damped oscillators are entrained more efficiently to external forcing than cells with sustained oscillations with different intrinsic periods. Gu et al.²⁷ investigated the range of entrainment for damped oscillators in different regions of the SCN. They found that the entrainment range widens when the proportion of damped oscillators in the light insensitive region of the SCN increases, whereas the entrainment range narrows when the proportion of damped oscillators in the light sensitive region increases.

To our knowledge, previous studies have not reported polyglot entrainment of damped oscillators.

Our results open a number of questions. First, 1:1 entrainment implies the presence of an input–output phase relationship. Classical studies on entrainment have paid little attention to how the phase of entrainment depends on the input frequency. In Perez-Cervera et al.,⁵⁴ it is shown that the timing of the input in relation to the phase of the output oscillation (i.e., whether it arrives before or after the crest of the oscillation) has a significant effect on the entrainment properties of the system. Motivated by these results, we find that in the FHN model with square-wave forcing, the phase of entrainment can either be positive (the oscillation peaks after the forcing turns on) or negative (the oscillation peaks before the forcing turns on) when the unforced system is a self-sustained oscillator and we have a monoglot entrainment structure [Fig. S20(A)]. However, when the unforced system is a damped oscillator, we have only observed positive phases of entrainment for both monoglot and polyglot entrainment structures [Figs. S20(B)–S20(D)]. Furthermore, the phase-shift structure (shapes of the phase of entrainment vs forcing period curves) depends on the type of 1:1 entrainment. Future work should address these issues in more detail by characterizing the phase-shift structure and understanding the underlying mechanisms.

Second, we have focused on models in which the unforced system is planar (two-dimensional). However, models with implications to realistic systems found in nature are expected to be higher-dimensional. Although analyzing the phase space of higher-dimensional models is more challenging, if one finds polyglot entrainment in a three-dimensional model it would be interesting to determine whether the geometric structure of the null surfaces in the phase-space diagram play a similar role to the knees of the cubic nullclines in two-dimensional models. Natural candidates to exhibit polyglot entrainment in three-dimensional models are these that have been shown to exhibit the canard phenomenon,^2,62,69,75 particularly extensions of the FHN to three-dimensions.

Third, polyglot entrainment is expected to have implications for circadian and neuronal systems. For example, in the context of the circadian clock, the period of external forcing is fixed at 24 hours but intrinsic periods of the circadian clock vary across individuals in a population. If the intrinsic period is closer to 24 hours, which is the period of forcing, then one may expect 1:1 entrainment. As the intrinsic period gets further away from 24 hours, we expect to lose 1:1 entrainment. However, if we have polyglot entrainment responses for a population of circadian clocks, then as the intrinsic period gets even further from 24 then 1:1 entrainment may be regained for some of the oscillators. Although there is much previous work on entrainment of circadian and related oscillators,^{1,6,7,10–12,14,16,23,25,26,30,33–35,41,48,51,54,63–65,67,70,81} we are not aware of any studies focused on polyglot entrainment of these oscillators.

Finally, within complex organisms, as well as some simple ones, cellular oscillators interact with each other to form systems of coupled oscillators.³⁵ For example, coupled circadian oscillations exist in multicellular cyanobacteria.³ In mammals, central circadian oscillators located in the SCN receive direct information about light–dark cycles via photic input pathways. These central oscillators in turn interact with circadian oscillators in peripheral tissues located throughout the body. Such systems are described as hierarchical coupled oscillators and how they entrain to light–dark cycles is an active area of research.⁴⁴ Modeling studies have shown that entrainment is achieved efficiently when damped oscillators are coupled to one another.^{5,24,38,45,74,76,78} In future work, it would be interesting to see if hierarchical coupled systems exhibit polyglot entrainment. Moreover, one could investigate polyglot entrainment responses in a network consisting of different types of oscillators. For example, does one observe polyglot entrainment with two coupled oscillators when one is a damped oscillator and the second is a self-sustained oscillator, or when one is a self-sustained oscillator and the second is non-oscillatory?

See the supplementary material which includes 27 figures and 2 movies illustrating a number of scenarios related to the main results of the paper.

This paper was partially funded by the NSF Grant Nos. CRCNS-DMS 1608077 (H.G.R.), NSF-DMS 1313861 (H.G.R.), and DMS-1555237 (C.O.D.).

The authors have no conflicts to disclose.

Data sharing is not applicable to this article as no new data were created or analyzed in this study. The codes used to generate the graphs in this model can be obtained from the corresponding authors.

In Secs. III B and III C, we observed polyglot entrainment in the periodically forced Fitzhugh–Nagumo model. Although the FHN model is a model of neuronal oscillations, the $w$ variable does not have a direct biophysical interpretation. This motivated us to explore, in Sec. III D, whether polyglot entrainment is observed in other two-dimensional models, for which the variables do have a direct biological or chemical interpretation. A brief description of these additional models of neuronal, circadian, and chemical oscillations is provided below.

The two-dimensional Morris–Lecar model⁵⁷ is a reduction of a model of oscillations in the barnacle giant muscle fiber,⁴⁹ and it has been used as a prototypical model of neuronal oscillations. The model has a cubic-like $V$-nullcline and a sigmoidal $w$-nullcline. This model can exhibit both type I and II neuronal excitability. Type I excitability occurs via a saddle-node on an invariant circle (SNIC) bifurcation, while type II excitability occurs via a Hopf bifurcation.⁵⁸ 

The periodically forced Morris–Lecar model is given by

In this model, $Cm$ is the membrane capacitance, $Iapp$ is the input current, $gL$ is the conductance for leak current, $EL$ is the reversal potential for leak current, $gK$ is the maximal conductance for $K+$ current, $EK$ is the reversal potential for $K+$ current, $gCa$ is the maximal conductance for $Ca2+$ current, $ECa$ is the reversal potential for $Ca2+$ current, and $V1$⁠, $V2$⁠, $V3$⁠, $V4$ are gating variable kinetic parameters, and $ϕ$ is a scaling parameter.

We used the following parameter values for all simulations: $Cm=20$ $μF/cm2$⁠, $EL=−60$ mV, $ECa=120$ mV, $EK=−84$ mV, $gL=2$ $μS/cm2$⁠, $gK=8$ $μS/cm2$⁠, $gCa=4$ $μS/cm2$⁠, $V1=−1.2$ mV, $V2=18$ mV, and $Iapp=95μA/cm2$⁠. Additionally, for type I excitability (SNIC), we set $ϕ=0.067$ $1/ms$⁠, $V3=12$ mV, and $V4=17.4$ mV. For type II excitability (Hopf), we set $ϕ=0.04$ 1/ms, $V3=2$ mV, and $V4=30$ mV.

The two-dimensional oregonator is a simplified version of the three-dimensional oregonator model, which represents chemical oscillations in the Belousov–Zhabotinsky reaction. The 2D oregonator serves as a reduced version of the Field–Koros–Noyes mechanism¹⁸ using a quasi-steady approximation.^19,79 This model features a cubic-like $v$-nullcline and a linear $w$-nullcline.

The periodically forced 2D oregonator is given by

The variables $v$ and $w$ represent the dimensionless concentrations of $HBRO2$ and Ce(IV), respectively. The parameter $η$ denotes a stoichiometric factor, $q$ is a rate parameter, and $ϵ$ is the time scale separation constant. We used the following parameters for all simulations: $η=2.28$ (fixed-point is a stable focus), $q=0.01$⁠, and $ϵ=0.025$⁠.

The NT model^53,71 represents a mechanism for circadian oscillations in Drosophila. It features a cubic-like $v$-nullcline and a linear $w$-nullcline.

The periodically forced Novak–Tyson model is given by

The variables $M$ and $Pt$ represent the dimensionless concentrations of mRNA concentration and total PER protein, respectively. In this model, $vm$ is the maximum rate of synthesis for $M$⁠, $km$ is the decay constant for $M$⁠, $kp1$ is the maximal rate for monomer phosphorylation, $kp2$ is the maximal rate for dimer phosphorylation, $kp3$ is the first order rate constant, $Keq$ is the equilibrium constant, $Pcrit$ is the dimer concern at the half-maximum transcription rate, and $Jp$ is the Michaelis constant for protein kinase.

We used the following parameters for all simulations: $vm=1$⁠, $km=0.1$⁠,$vp=0.5$⁠, $kp1=10$⁠, $kp2=0.03$⁠, $kp3=0.1$⁠, $Keq=3.3$⁠, $Pcrit=0.1$⁠, and $Jp=0.05$⁠.

The Lengyel–Epstein model^42,43 is a model for the chlorite-iodide-malonic acid (CIMA) reaction. The periodically forced Lengyel–Epstein model is given by

The variables $x$ and $y$ represent the concentrations of iodide and chlorite, respectively. $a$ and $b$ are the feed concentration parameters. We used the following parameter values: $σ=8$⁠, $b=1$⁠, and $a=15$ or $16$⁠.