# Awakened Oscillations in Coupled Consumer-Resource Pairs

Authors: Almaz Mustafin
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

The paper concerns two interacting consumer-resource pairs based on chemostat-like equations under the assumption that the dynamics of the resource is considerably slower than that of the consumer. The presence of two different time scales enables to carry out a fairly complete analysis of the problem. This is done by treating consumers and resources in the coupled system as fast-scale and slow-scale variables, respectively, and subsequently considering developments in phase planes of these variables, fast and slow, as if they are independent. When uncoupled, each pair has unique asymptotically stable steady state and no self-sustained oscillatory behavior (although damped oscillations about the equilibrium are admitted). When the consumer-resource pairs are weakly coupled through direct reciprocal inhibition of consumers, the whole system exhibits self-sustained relaxation oscillations with a period that can be significantly longer than intrinsic relaxation time of either pair. It is shown that the model equations adequately describe locally linked consumer-resource systems of quite different nature: living populations under interspecific interference competition and lasers coupled via their cavity losses.


---

## Body

## 1. Introduction

Recently, there has been a great deal of activity aimed at studying the synchronization of coupled oscillators of diverse nature [1–5]. The theory of synchronization implies that even in uncoupled state the individual elementary units exhibit self-sustained oscillations. However, no less interesting are the systems where local coupling is essential for the very generation of oscillations and not only for their modulation or phase adjustment.As far back as in early 1970s, Smale [6] constructed a counterintuitive mathematical example of a biological cell modeled by the chemical kinetics of four metabolites, x
1
,
…
,
x
4, such that the reaction equations d
x
/
d
t
=
R
(
x
) for the set of metabolites, x
=
(
x
1
,
…
,
x
4
), had a globally stable equilibrium. The cell is “dead,” in that the concentrations of its metabolites always tend to the same fixed levels. When two such cells are coupled by linear diffusion terms of the form M
(
x
2
-
x
1
), where M is a diagonal matrix with the elements μ
k
δ
k
l, however, the resulting equations are shown to have a globally stable limit cycle. The concentrations of the metabolites begin to oscillate, and the system becomes “alive.” In Smale’s words, “there is a paradoxical aspect to the example. One has two dead (mathematically dead) cells interacting by a diffusion process which has a tendency in itself to equalize the concentrations. Yet in interaction, a state continues to pulse indefinitely.”The reaction equations involved in Smale’s model were too general to appeal to any specific process. Since then, inspired by his pioneer work, a number of models have been proposed containing biologically plausible mechanisms by which coupling of identical nonoscillating cells could generate synchronous oscillations. Among them are a model of electrically coupled cells, characterized by an excitable membrane and calcium dynamics [7]; a model in which coupling a passive diffusive cytoplasmic bulk with an excitable membrane (having an activator-inhibitor dynamics) produces a self-sustained oscillatory behavior [8]; an analog operational-amplifier implementation of neural cells connected by passive coupling (where conductance of the resistive connection simulates the diffusion coefficient) [9], and so forth. The authors of the last model suggested a term “awakening dynamics” for the phenomenon.The subject of the present paper is an emergence of collective oscillations in a system of coupled nonoscillatory consumer-resource (CR) pairs. Owing to simplicity, this model in a sense may be considered minimal. Our choice of coupled CR equations as a matter of enquiry is dictated primarily by the ubiquity and importance of CR interactions.CR models are the fundamental building blocks used in mathematical description and simulation of ecosystems. Depending on a specific nature of the involved CR interactions, they can take the forms of predator-prey, plant-herbivore, parasite-host, and victim-exploiter systems [10]. However applications of the CR models extend far beyond the ecology and are found wherever one can speak of win-loss interactions. In its broad meaning, resource is any substance which can lead to increased growth rate of the consumer as its availability in the environment is increased. As this takes place, the resource is certainly consumed. Consuming the resource means tending to reduce its availability. When carefully examined, CR models are identified in the following fields: epidemiology (susceptible and infected [11, ch. 10]), laser dynamics (photons and electrons [12, ch. 6]), labor economics (share of labor and employment rate [13, p. 28]), theoretical immunology (antigens and B lymphocytes [14, p. 299]), kinetics of chain chemical reactions (lipid molecules and free radicals [15]), and in numerous other studies from diverse disciplines.We consider the situation when each of two consumer species exploits one respective resource only. As explained in the Appendix, terms “consumer” and “resource” in our model may bear not only their literal ecological meaning, but the physical meaning of photon density and population inversion in a laser cavity as well. Both resources are being supplied with constant rates like in a chemostat and consumed according to a simple mass-action kinetics. The resources are thought to be noninteractive. When uncoupled, self-inhibition of the individual consumer population is due to intraspecific interference. The coupling is assumed to originate solely from the interspecific interference competition between the consumers and quantitatively expressed by a bilinear term combining the competitor densities. Thus, the per individual loss rate of either consumer is proportional to the density of its counterpart.Representation of competition between species in terms of loss-coupling dates back to the classical model of Lotka-Volterra-Gause (LVG) [16]. The LVG model operates with carrying capacities of the species, rather than referring explicitly to any essential resources. As shown by MacArthur [17], LVG equations may be considered as a quasi-steady-state approximation to the CR equations accentuating resource-mediated nature of competition, under the assumption of relatively rapid dynamics of the involved resources. Thereafter trophic competition have developed into a major descriptor of competition in the ecological literature generally and in conceptualizing ecosystems as systems of coupled CR oscillators specifically [5]. In contrast to the prevailing models of competition we consider the case of pure interspecific interference competition between the consumers with no consumption-induced contribution. Actually, our model is nothing more nor less than LVG equations augmented with the rate equations for the resources. Another key assumption of the model is that the dynamics of the consumers is much faster than that of the resources.From a physical perspective, by and large similar equations with the like time hierarchy (fast consumer and slow resource) have been in use for coupled longitudinal modes in a semiconductor laser with an intracavity-doubling crystal since the work of Baer [18]. These equations have been treated mostly numerically. The notable analytical result belongs to Erneux and Mandel [19] who succeeded to show that the system admits antiphase periodic solutions by reducing it to the equations for quasi-conservative oscillator. However this result has to do with the onset of low-amplitude quasi-harmonic oscillations. Unlike their study, our approach deals with well-developed high-amplitude essentially nonlinear oscillations. Besides, we propose the model to be valid not only for competing laser modes, but for loss-coupled lasers as well.We analyze the model using geometric singular perturbation technique according to which the full system of equations is decomposed into fast and slow subsystems. As we shall see below, the model reveals qualitatively different behavior at intense and weak competition between the consumer species. If coupling is strong, one of the consumers wins and completely dominates. When coupling is weak, the model exhibits low-frequency antiphase relaxation oscillations with each species alternatively taking the dominant role.


## 2. The Model

The two-consumer, two-resource model we consider is the following nondimensional system of four ordinary differential equations:(1a)
u
˙
1
=
γ
1
-
(
u
1
+
1
)
v
1
-
u
1
,
(1b)
u
˙
2
=
γ
2
-
(
u
2
+
1
)
v
2
-
u
2
,
(1c)
ε
v
˙
1
=
(
u
1
-
δ
v
1
-
ϰ
2
v
2
)
v
1
,
(1d)
ε
v
˙
2
=
(
u
2
-
δ
v
2
-
ϰ
1
v
1
)
v
2
.Here dots indicate differentiation with respect to the nondimensional time variable t, u
i and v
i  (i
=
1,2) are quantities measuring the respective population sizes of ith resource and ith consumer, γ
i
>
0  (i
=
1,2) is the inflow rate of ith resource, 0
<
δ
≪
1 is a parameter representing consumer self-limitation, ϰ
j
>
0  (j
=
1,2 and j
≠
i) quantifies the inhibitory effect of jth consumer on the growth of ith consumer due to coupling, and 0
<
ε
≪
1 is a singular perturbation parameter indicating that the dynamics of the consumers is much faster than that of the resources.It should be mentioned that, being proportional to its dimensional prototype,v
i directly represents population density of consumer species and is always nonnegative. Quantity u
i, however, is not a population size in the true sense of the word. It is rather an affine transformation of a population size of the form N
→
a
N
+
b, where a and b are scaling constants. This is done for reasons of mathematical convenience. Unlike a purely linear transformation, an affine map does not preserve the zero point, so in (1a)–(1d)  u
i
=
-
1 corresponds to zero population size in reality. Nevertheless, from here on we shall apply the term “resource” to u
i for brevity.For more details and discussion on the derivation of the model (1a)–(1d) from different perspectives the reader is referred to the Appendix.


## 3. Model Analysis and Implications

### 3.1. A Single CR Pair

Whenϰ
1,2
=
0, the communities are uncoupled and completely independent. An isolated CR pair is governed by equations:(2a)
u
˙
=
γ
-
(
u
+
1
)
v
-
u
,
(2b)
ε
v
˙
=
(
u
-
δ
v
)
v
.There exist two nonnegative steady states:(3a)
u
¯
=
γ
,
v
¯
=
0
;
(3b)
u
¯
=
1
2
(
1
+
(
4
γ
+
2
+
δ
)
δ
-
1
-
δ
)
=
γ
δ
+
O
(
δ
2
)
,
v
¯
=
1
2
δ
(
1
+
(
4
γ
+
2
+
δ
)
δ
-
1
-
δ
)
=
γ
-
γ
(
γ
+
1
)
δ
+
O
(
δ
2
)
.The linearization of (2a) and (2b) takes the form
(4)
J
=
(
-
v
¯
-
1
-
u
¯
-
1
v
¯
ε
u
¯
-
2
δ
v
¯
ε
)
,
where (u
¯, v
¯) is one of the above steady states (3a) and (3b).At (3a), one eigenvalue is negative and one is positive: λ
1
=
-
1, λ
2
=
γ
/
ε. Thus (3a) is a saddle point. At (3b), Tr
⁡
J
=
-
(
δ
/
ε
+
1
)
v
¯
-
1
<
0 and det
⁡
J
=
(
1
+
2
δ
v
¯
+
δ
)
v
¯
/
ε
>
0, so the steady state is a stable node/focus. Specifically, focus is the case for
(5)
(
Tr
⁡
J
)
2
-
4
det
⁡
J
=
(
γ
2
δ
2
+
O
(
δ
3
)
)
ε
-
2
+
(
-
4
γ
+
O
(
δ
)
)
ε
-
1
+
O
(
1
)
<
0
,
whence one obtains an asymptotic estimate
(6)
δ
=
o
(
ε
1
/
2
)
.
Figure 1(c) illustrates this focus-node bifurcation numerically. Damped oscillations is a well-known inherent feature of the photon-carrier dynamics in class-B lasers.Phase portraits of an isolated CR system (2a) and (2b) for two cases when the unique stable equilibrium (3b) is a node (a), and a focus (b). Shown are vertical (u
˙
=
0) and horizontal (v
˙
=
0) nullclines. Steady states are marked by open circles. Numerical simulations have been carried out with data from Table 2. The respective values of the second order loss parameter δ in (a) and (b) are 0.1 and 0.01. Typically, consumer scale (v) is enormous; the equilibrium point happens to lie very close to the origin. For better appearance, henceforward we display all the graphs using coordinate transformation u
→
arsinh
(
u
/
γ
δ
) and v
→
arsinh
(
v
/
γ
). (c) Eigenvalues of (4) at (3b) as a function of the second order loss δ. Critical value δ
cr indicates a boundary between the regions of damped oscillations and aperiodic damping.
(a)
(b)
(c)An isolated CR system (2a) and (2b) admits, therefore, the only solution which tends asymptotically towards the unique steady state. Periodic solutions are excluded. However the temporal dynamics of approaching this steady state essentially depends on the interplay between ε and δ, and is worth another look. There are in fact three timescales, O
(
ε
), O
(
1
), and O
(
δ
), involved in the CR system (2a) and (2b) when it is overdamped, and two—O
(
ε
) and O
(
1
)—when underdamped.System (2a) and (2b) is singularly perturbed because the derivative of one of its state variables, v, is multiplied by a small positive parameter ε. Singular perturbation cause two-time-scale behavior of the system characterized by the presence of slow and fast transients in the system’s response to external stimuli.Replacingt in (2a) and (2b) with a fast time variable τ
=
t
/
ε and setting ε
=
0 we obtain the fast subsystem:(7a)
u
′
=
0
,
(7b)
v
′
=
(
u
-
δ
v
)
v
,where prime means differentiation with respect to τ. In the stretched timescale τ the slow resource variable u, according to (7a), is replaced by its initial value and reckoned as constant parameter. Equation (7b) is of a logistic type which has the solution
(8)
v
(
τ
)
=
u
/
δ
1
+
(
u
/
δ
v
(
0
)
-
1
)
exp
⁡
⁡
(
-
u
τ
)
,
valid for t
=
O
(
ε
).After a lapse of considerable time (in the fast scale)v converges to either of two fixed points depending on a sign of u:
(9)
lim
⁡
τ
→
∞
v
(
τ
)
=
{
u
δ
,
if
u
>
0
;
0
,
if
u
<
0
.
It means that every single trajectory starting within the positive quadrant of (
u
,
v
) plane far enough from the stable steady state (3b) will run almost parallel to the vertical axis and hit the line u
-
δ
v
=
0 practically in a finite time of order O
(
ε
). This is shown in Figures 1(a) and 1(b).Now setε
=
0 in (2a) and (2b) to get the slow subsystem:(10a)
u
˙
=
γ
-
(
u
+
1
)
v
-
u
,
(10b)
0
=
(
u
-
δ
v
)
v
.Equation (10b) describes a slow manifold consisting of two lines in the (
u
,
v
) plane: v
=
u
/
δ and v
=
0. By (9), the former attracts all the trajectories in the first quadrant, while the latter—all those in the second. In as much as the quasi-steady state of (7b), v
¯
=
u
/
δ, is an isolated root of (10b) and v
¯ is a stable solution of (10b) for any u
>
0; the assumptions of Tikhonov’s theorem [20] are satisfied, and one may proceed to approximate u and v in terms of the solution of the reduced system (10a) and (10b) in the slow timescale. It means that after arriving at the slow manifold v
=
u
/
δ, the representing point of the full system (2a) and (2b) will move along the manifold toward the equilibrium point with a characteristic velocity of order O
(
1
).In the immediate proximity to equilibrium (3b) the behavior of the trajectory is determined by the type of the fixed point, whether it is a node or a focus. Eventually this depends on the value of δ. Namely, close to a stable node, the system has two distinct real negative eigenvalues, one fast (λ
1), and one slow (λ
2):
(11)
λ
1
=
(
-
γ
δ
+
O
(
δ
2
)
)
ε
-
1
+
(
δ
-
1
+
O
(
1
)
)
+
O
(
ε
)
,
λ
2
=
(
-
δ
-
1
+
O
(
1
)
)
+
O
(
ε
)
,
for
1
>
δ
≫
ε
.
Since |
λ
1
|
≫
|
λ
2
|, trajectories starting off the associated eigenvector 2 (which is tangent to the slow manifold v
=
u
/
δ) converge to that vector along lines almost parallel to eigenvector 1 (which is parallel to the vertical axis v). As they approach vector 2 they become tangent to it and move along it up to the very nodal point (Figure 1(a)). The characteristic time constant of this final stage is of order O
(
δ
).In the vicinity of a stable focus the motion is qualitatively different: trajectories still keep converging to the equilibrium, but no longer follow the slow manifoldv
=
u
/
δ (Figure 1(b)). The reason is that for a focus eigenvectors are complex. Recalling that condition (6) must be true for a focus, calculate the eigenvalues in a limit case of δ
→
0:
(12)
λ
1,2
=
-
1
2
(
γ
+
1
)
±
i
γ
ε
for
δ
=
0
.
Loosely speaking, one may think that near a focal point separation of state variables into slow and fast ones ceases to have its conventional meaning. There is no point to talk about motion along the slow manifold since there are no reduced one-dimensional systems corresponding to the neighborhood of a focus. The perturbation parameter ε affects only the frequency of damped oscillation, but not the damping rate. The radius of the focal spiral uniformly shrinks with a time constant of order O
(
1
).ToO
(
1
) for small δ there is a way to recast (2a) and (2b) near the focal equilibrium (
0
,
γ
) in a convenient form where the perturbation parameter ε does not multiply any right hand side. Namely, performing the scaling t
=
μ
2
s, u
=
μ
2
γ
ξ, and v
=
γ
(
1
+
η
), where μ
2
=
ε
/
γ, one obtains
(13)
ξ
˙
=
-
η
-
μ
2
ξ
(
1
+
γ
(
1
+
η
)
)
,
η
˙
=
ξ
(
1
+
η
)
.
Here dots stand for differentiation with respect to the time variable s.Equations (13) represent a weakly perturbed Hamiltonian system. For μ
=
0 the system is pure Hamiltonian and admits the first integral
(14)
H
=
1
2
ξ
2
+
η
-
ln
⁡
⁡
(
1
+
η
)
,
which is a conserved quantity (H
˙
=
0). The periodic solutions of the Hamiltonian system form a one-parameter family with the equilibrium (
0,0
) as center point. The condition 0
<
μ
≪
1 makes (13) quasi-conservative with phase trajectories slowly spiralling to the equilibrium (Figure 2).Figure 2
Phase-plane trajectory of the weakly perturbed Hamiltonian system (13) for γ
=
1.19375, ε
=
0.727273
×
1
0
-
3, and μ
=
0.157107. Near the equilibrium the oscillations are almost harmonic.System (13) can be rewritten as a second order differential equation for η only:
(15)
η
¨
+
η
+
η
2
-
η
˙
2
1
+
η
+
μ
2
(
1
+
γ
(
1
+
η
)
)
η
˙
=
0
.
This is the equation for a nonlinear quasi-conservative oscillator. As is seen from (15), for |
η
|
≪
1 and |
η
˙
|
≪
1 the oscillator is not only quasi-conservative, but also quasi-linear with the frequency ω
0
=
1. These conditions motivate introducing the new variable ζ
=
η
/
μ. As a result one obtains
(16)
ζ
¨
+
ζ
=
μ
(
-
ζ
2
+
ζ
˙
2
1
+
μ
ζ
-
μ
(
1
+
γ
+
μ
γ
ζ
)
ζ
˙
)
.
This equation of quasi-conservative and quasi-linear oscillator will be used later on when analyzing the onset of synchronous periodicity.


### 3.2. Coupled Communities: Stability

In the following, we will distinguish two main types of coupling: strong,ϰ
1,2
>
1, and weak, ϰ
1,2
<
1.Physically feasible equilibria, or fixed points,(
u
¯
1
,
u
¯
2
,
v
¯
1
,
v
¯
2
), of (1a)–(1d) are those for which v
¯
1,2
⩾
0. We denote the interior fixed point by F
12
=
(
u
¯
1
,
u
¯
2
,
v
¯
1
,
v
¯
2
), where the subscripts stand for the consumers. Lack of a certain index at a boundary fixed point means that the consumer concerned is not present (extinct). Thus F
1
=
(
u
¯
1
,
u
¯
2
,
v
¯
1
,
0
) and F
2
=
(
u
¯
1
,
u
¯
2
,
0
,
v
¯
2
) designate either of one-consumer equilibria corresponding to dominance, while F
=
(
u
¯
1
,
u
¯
2
,
0,0
) means both consumers having been washed out.Model (1a)–(1d) has four feasible steady states. To O
(
1
) for small δ
(17a)
F
:
u
¯
1
=
γ
1
,
u
¯
2
=
γ
2
,
v
¯
1
=
0
,
v
¯
2
=
0
;
(17b)
F
1
:
u
¯
1
=
0
,
u
¯
2
=
γ
2
,
v
¯
1
=
γ
1
,
v
¯
2
=
0
;
(17c)
F
2
:
u
¯
1
=
γ
1
,
u
¯
2
=
0
,
v
¯
1
=
0
,
v
¯
2
=
γ
2
;
F
12
:
u
¯
1
=
ϰ
1
γ
1
-
ϰ
2
γ
2
-
ϰ
1
ϰ
2
+
1
±
R
2
(
ϰ
1
-
1
)
,
u
¯
2
=
-
ϰ
1
γ
1
+
ϰ
2
γ
2
-
ϰ
1
ϰ
2
+
1
±
R
2
(
ϰ
2
-
1
)
,
(17d)
v
¯
1
=
u
¯
2
ϰ
1
,
v
¯
2
=
u
¯
1
ϰ
2
,where
(18)
R
=
(
ϰ
1
γ
1
-
ϰ
2
γ
2
-
ϰ
1
ϰ
2
+
1
)
2
+
4
ϰ
2
(
ϰ
1
-
1
)
(
ϰ
1
γ
1
-
γ
2
)
.Boundary equilibriaF, F
1, and F
2 always exist. F
12 exists if
(19)
1
ϰ
2
<
γ
2
γ
1
<
ϰ
1
for
ϰ
1
>
1
∧
ϰ
2
>
1
(
strong
coupling
)
,
or
(20)
ϰ
1
<
γ
2
γ
1
<
1
ϰ
2
for
ϰ
1
<
1
∧
ϰ
2
<
1
(
weak
coupling
)
.
In case (19) the expression with “+” in (17d) is realized, while in case (20) the one with “−” is feasible.F is always unstable, because two of the four associated eigenvalues are positive:
(21)
λ
(
F
)
:
γ
1
ε
,
γ
2
ε
,
-
1
,
-
1
.Correct toO
(
1
) in ε, the eigenvalues for F
1 and F
2 are(22a)
λ
(
F
1
)
:
-
1
,
-
1
2
(
γ
1
+
1
)
±
i
γ
1
ε
,
γ
2
-
ϰ
1
γ
1
ε
;
(22b)
λ
(
F
2
)
:
1
,
-
1
2
(
γ
2
+
1
)
±
i
γ
2
ε
,
γ
1
-
ϰ
2
γ
2
ε
,whence it follows that (17b) and (17c) are stable when(23a)
γ
2
γ
1
<
ϰ
1
,
(23b)
γ
1
γ
2
<
ϰ
2
,respectively.The necessary and sufficient conditions for all the eigenvalues of the Jacobian matrix,(24)
(
-
v
¯
1
-
1
0
-
ϰ
2
v
¯
2
-
1
0
0
-
v
¯
2
-
1
0
-
ϰ
1
v
¯
1
-
1
v
¯
1
ε
0
0
-
ϰ
2
v
¯
1
ε
0
v
¯
2
ε
-
ϰ
1
v
¯
2
ε
0
)
,
evaluated at F
12, to have negative real parts are, from the Routh-Hurwitz criterion, (25a)
c
0
>
0
,
(25b)
c
3
>
0
,
(25c)
c
2
c
3
-
c
1
>
0
,
(25d)
c
1
(
c
2
c
3
-
c
1
)
-
c
0
c
3
2
>
0
,where c
0, c
1, c
2, and c
3 are the coefficients of the characteristic polynomial λ
4
+
c
3
λ
3
+
c
2
λ
2
+
c
1
λ
+
c
0 of (24):
(26)
c
0
=
v
¯
1
v
¯
2
(
ϰ
1
(
v
¯
1
-
ϰ
2
(
v
¯
1
+
v
¯
2
+
1
)
)
+
ϰ
2
v
¯
2
+
1
)
ε
-
2
,
1
c
1
=
-
ϰ
1
ϰ
2
v
¯
1
v
¯
2
(
v
¯
1
+
v
¯
2
+
2
)
ε
-
2
+
(
v
¯
1
(
v
¯
2
(
ϰ
1

+

ϰ
2
+
ϰ
1
v
¯
1

+

ϰ
2
v
¯
2

+

2
)
+
1
)
+
v
¯
2
)
ε
-
1
,
c
2
=
-
ϰ
1
ϰ
2
v
¯
1
v
¯
2
ε
-
2
+
(
v
¯
1
(
v
¯
2
(
ϰ
1
+
ϰ
2
)
+
1
)
+
v
¯
2
)
ε
-
1
+
(
v
¯
1
+
1
)
(
v
¯
2
+
1
)
,
c
3
=
v
¯
1
+
v
¯
2
+
2
.
11111111111111111111111111111To analyze the validity of (25a)–(25d) we put ϰ
2
=
μ
ϰ
1, where μ
=
O
(
1
). Then for (25a) to be true, ϰ
1
=
o
(
1
) should be met, which is incompatible with strong coupling, yet possible for weak coupling. Conditions (25b) and (25c) are always valid. As is known, (25d) guarantees a simple complex conjugate pair of eigenvalues corresponding to a linearization about steady state F
12 to have negative real part. For ϰ
1
=
o
(
1
) it boils down to
(27)
(
-
v
¯
1
v
¯
2
(
v
¯
1
+
v
¯
2
+
2
)
(
v
¯
1
+
v
¯
2
+
v
¯
1
2
+
v
¯
2
2
)
μ
ϰ
1
2
+
O
(
ϰ
1
3
)
)
ε
-
3
+
(
(
v
¯
1
+
1
)
(
v
¯
2
+
1
)
(
v
¯
1
-
v
¯
2
)
2
+
O
(
ϰ
1
)
)
ε
-
2
+
O
(
ε
-
1
)
>
0
,
whence it follows that
(28)
ϰ
1
=
o
(
ε
1
/
2
)
.
With regard to a fairly small value of ε, (28) may be thought to be broken under most physically meaningful conditions unless coupling is infinitesimally weak. Hence, normally, condition (25d) of the Routh-Hurwitz criterion is never fulfilled and the interior fixed point F
12—if it exists—is always unstable by growing oscillations.Existence and stability conditions of possible nonnegative equilibrium points are summarized in Table1.Table 1
Existence and stability conditions of nonnegative equilibria of system (1a)–(1d).
Equilibrium
Existence
Stability
F
Always
Unstable
F
1
Always
γ
2
/
γ
1
<
ϰ
1
F
2
Always
γ
2
/
γ
1
>
1
/
ϰ
2
F
12
1
/
ϰ
2
<
γ
2
/
γ
1
<
ϰ
1 for ϰ
1,2
>
1 (strong coupling)
Unstable
ϰ
1
<
γ
2
/
γ
1
<
1
/
ϰ
2 for ϰ
1,2
<
1 (weak coupling)
ϰ
1,2
=
o
(
ε
1
/
2
)


### 3.3. Strong Coupling: Bistability and Hysteresis

Strong coupling in (1a)–(1d) makes bistability of boundary equilibria possible, as evident from Table 1. When both F
1 and F
2 are stable with an unstable coexistence steady state F
12, the system being studied is able to exhibit a hysteresis effect. Given the strong coupling, suppose that the inflow of resource 2 is kept at some constant level, γ
2
=
γ
2
*, while the inflow of resource 1, γ
1, steadily increases from a value less than γ
2
*
/
ϰ
1 along the line γ
2
=
γ
2
* in the (
γ
1
,
γ
2
) parameter plane, as shown in Figure 3(a). Referring to (17a)–(17d) and (23a) and (23b), one sees that initially F
1 is unstable; F
2 is stable and F
12 does not exist with the complete dominance of consumer 2. Within the interval γ
2
*
/
ϰ
1
<
γ
1
<
ϰ
2
γ
2
* steady state F
1 becomes stable yet empty and F
12 becomes existent (according to (19)) yet unstable, so the situation remains unchanged until point (
ϰ
2
γ
2
*
,
γ
2
*
) in Figure 3(a) has been reached from the left. For a larger γ
1, state F
2 gives up its stability and the system jumps to F
1. Consumer 2 gets washed out, while consumer 1 takes over. If now we start reducing γ
1, the system remains in F
1 until γ
1 drops to the lower critical value γ
2
*
/
ϰ
1, beyond which F
1 is no longer stable and there is a reverse jump to F
2. In other words, as γ
1 progresses along γ
2
=
γ
2
* there is a discontinuous switch from consumer 2 to consumer 1 at ϰ
2
γ
2
*, while as γ
1 retraces its steps, there is a discontinuous switch from consumer 1 to consumer 2 at γ
2
*
/
ϰ
1. Figure 3(b) illustrates how the steady-state level of consumer 1 responds to infinitesimally slow changes of the inflow rate of its own resource. The hysteresis is made possible thanks to the concurrent stability of both boundary equilibria and instability of the interior fixed point for γ
1
∈
(
γ
2
*
/
ϰ
1
,
ϰ
2
γ
2
*
). In terms of electronics, such a situation would describe a flip-flop circuit—bistable multivibrator—having two stable conditions, each corresponding to one of two alternative input signals.Hysteresis in the two-consumer, two-resource system (1a)–(1d). It takes place in the case of strong coupling ϰ
1,2
>
1. (a) The resource-supply parametric plane (
γ
1
,
γ
2
). Consumer 1 completely dominates below the line γ
1
-
ϰ
2
γ
2
=
0, whereas consumer 2—above the line ϰ
1
γ
1
-
γ
2
=
0. Both boundary steady states are stable in the region of bistability confined by the two aforementioned lines and marked off by gray, with the dominance of either consumer being a matter of path-dependency. (b) Equilibrium level of consumer 1 as a function of the resource supply. Numerical values of the coupling strengths are chosen to be ϰ
1
=
1.2 and ϰ
2
=
1.5.
(a)
(b)


### 3.4. Weak Coupling: Antiphase Relaxation Oscillations

As seen from Table1, the very existence of interior equilibrium F
12 in a case of weak coupling (condition (20)) implies instability of both boundary fixed points, F
1 and F
2. System (1a)–(1d) happens to possess four nonnegative steady states, none of them being stable. As we have found, F
12 is unstable through growing oscillations. In such a case, the model would thus be expected to have a limit cycle in its four-dimensional phase space corresponding to self-sustained oscillations.As is known, Hopf bifurcations come in both super- and subcritical types. If a small, attracting limit cycle appears immediately after the fixed point goes unstable, and if its amplitude shrinks back to zero as the parameter is reversed, the bifurcation is supercritical; otherwise, it’s probably subcritical, in which case the nearest attractor might be far from the fixed point, and the system may exhibit hysteresis as the parameter is reversed.To check whether the bifurcation is supercritical or subcritical, we employ the averaging method [21]. For negligible δ the problem reduces to weakly coupled quasi-conservative oscillators. Introduce the new time s
=
t
/
μ
2, where μ
2
=
ε
/
γ
1. The new dynamic variables ξ
1,2 and η
1,2 are defined by formulas ξ
1,2
=
u
1,2
/
μ
2
γ
1,2 and η
1,2
=
(
v
1,2
-
γ
1,2
)
/
γ
1,2. Thus the variables ξ
1,2 and η
1,2 are the respective deviations of u
1,2 and v
1,2 from their standalone steady-state values. With these new variables equations (1a)–(1d) become(29a)
ξ
˙
1
=
-
η
1
-
μ
2
ξ
1
(
1
+
γ
1
(
1
+
η
1
)
)
,
(29b)
η
˙
1
=
ξ
1
(
1
+
η
1
)
-
(
γ
2
γ
1
)
ρ
2
(
1
+
η
1
)
(
1
+
η
2
)
,
(29c)
ξ
˙
2
=
-
η
2
-
μ
2
ξ
2
(
1
+
γ
2
(
1
+
η
2
)
)
,
(29d)
η
˙
2
=
(
γ
2
γ
1
)
ξ
2
(
1
+
η
2
)
-
ρ
1
(
1
+
η
1
)
(
1
+
η
2
)
,where dots mean differentiation with respect to s and ρ
1,2
=
ϰ
1,2
/
μ
2 are the scaled coupling strengths.For brevity, we restrict our consideration to the case of identical CR pairs and symmetric coupling settingγ
1,2
=
γ and ρ
1,2
=
ρ. With these assumptions system (29a)–(29d) can be transformed to two coupled second-order equations for quasi-conservative, quasi-linear oscillators. This can be done as follows. First, (29b) and (29d) are solved for ξ
1 and ξ
2. Secondly, (29b) and (29d) are differentiated with respect to time to get η
¨
1 and η
¨
2. Thirdly, ξ
1, ξ
2, ξ
˙
1, and ξ
˙
2 are plugged into the equations for η
¨
1 and η
¨
2. And finally, η
1,2 are replaced by the new variables ζ
1,2: η
1,2
=
μ
ζ
1,2.As a consequence, we arrive at the following coupled equations:(30)
ζ
¨
1,2
+
ζ
1,2
=
-
ρ
ζ
˙
2,1
+
μ
Z
1
(
1,2
)
+
μ
2
Z
2
(
1,2
)
+
O
(
μ
3
)
,
where
(31)
Z
1
(
1,2
)
=
-
ρ
(
γ
+
1
)
-
ζ
1,2
2
+
ζ
˙
1,2
2
-
ρ
ζ
1,2
ζ
˙
2,1
,
Z
2
(
1,2
)
=
-
(
γ
+
1
)
(
ζ
˙
1,2
+
ρ
ζ
2,1
)
-
ζ
1
(
ρ
(
2
γ
+
1
)
+
ζ
˙
1,2
2
)
.We seek the solution of (30) in the quasi-harmonic form
(32)
ζ
1,2
=
a
1,2
cos
⁡
⁡
(
ω
s
+
φ
1,2
)
,
ζ
˙
1,2
=
-
a
1,2
ω
sin
⁡
(
ω
s
+
φ
1,2
)
,
where a
1,2 and φ
1,2 are slowly varying amplitudes and phases, and ω is the frequency of the synchronous oscillations. Differentiating the assumed form of ζ
1,2 and equating the result to the assumed form of ζ
˙
1,2 yields the first pair of relationships between a
1,2 and φ
1,2:
(33)
a
˙
1,2
cos
⁡
⁡
(
ω
s
+
φ
1,2
)
-
a
1,2
φ
˙
1,2
sin
⁡
(
ω
s
+
φ
1,2
)
=
0
.Then, differentiatingζ
˙
1,2 and substituting the resulting expression for ζ
¨
1,2 as well as the assumed forms for ζ
1,2 and ζ
˙
1,2 into (30) yields the second pair of equations relating a
1,2 and φ
1,2. Separating into equations for the rate of change of a
1,2 and φ
1,2 one obtains
(34)
a
˙
1,2
=
ω
-
1
sin
⁡
(
ω
s
+
φ
1,2
)
(
A
0
(
1,2
)
+
μ
A
1
(
1,2
)
+
μ
2
A
2
(
1,2
)
)
,
φ
˙
1,2
=
a
1,2
-
1
ω
-
1
cos
⁡
⁡
(
ω
s
+
φ
1,2
)
(
A
0
(
1,2
)
+
μ
A
1
(
1,2
)
+
μ
2
A
2
(
1,2
)
)
,
where
(35)
A
0
(
1,2
)
=
-
a
1,2
(
ω
2
-
1
)
cos
⁡
⁡
(
ω
s
+
φ
1,2
)
-
ρ
a
2,1
ω
sin
⁡
(
ω
s
+
φ
2,1
)
,
A
1
(
1,2
)
=
a
1,2
(
-
a
1,2
ω
2
sin
2
(
ω
s
+
φ
1,2
)
+
a
1,2
cos
⁡
2
(
ω
s
+
φ
1,2
)
-
ρ
a
2,1
ω
cos
⁡
⁡
(
ω
s
+
φ
1,2
)
sin
⁡
(
ω
s
+
φ
2,1
)
-
a
1,2
ω
2
sin
2
)
+
ρ
(
γ
+
1
)
,
A
2
(
1,2
)
=
a
1,2
cos
⁡
(
ω
s
+
φ
1,2
)
×
(
a
1,2
2
ω
2
sin
2
(
ω
s
+
φ
1,2
)
+
ρ
(
2
γ
+
1
)
)
+
(
γ
+
1
)
×
(
ρ
a
2,1
cos
⁡
⁡
(
ω
s
+
φ
2,1
)
-
a
1,2
ω
sin
⁡
(
ω
s
+
φ
1,2
)
)
.To this point no approximations have been made except of the expansion in powers ofμ in (30). Averaging equations (34) over the period 2
π
/
ω and considering a
1,2, φ
1,2, a
˙
1,2 and φ
˙
1,2 to be constants while performing the averaging, one obtains the following equations describing the slow variations of a
1,2 and φ
1,2: (36a)
a
˙
1,2
=
1
2
π
∫
0
2
π
/
ω
‍
(
A
0
(
1,2
)
+
μ
A
1
(
1,2
)
+
μ
2
A
2
(
1,2
)
)
×
sin
⁡
(
ω
s
+
φ
1,2
)
d
s
=
-
ρ
a
2,1
ω
cos
⁡
⁡
ϕ
+
μ
2
(
γ
+
1
)
(
a
1,2
ω
-
ρ
a
2,1
sin
⁡
ϕ
)
2
ω
,
(36b)
φ
˙
1,2
=
1
2
π
a
1,2
∫
0
2
π
/
ω
‍
(
A
0
(
1,2
)
+
μ
A
1
(
1,2
)
+
μ
2
A
2
(
1,2
)
)
×
cos
⁡
⁡
(
ω
s
+
φ
1,2
)
d
s
=
ρ
(
a
1,2
2
+
a
2,1
2
)
sin
⁡
ϕ
2
a
1,2
a
2,1
+
O
(
μ
2
)
,where ϕ
=
φ
1,2
-
φ
2,1.Any stable fixed points of (36a) and (36b) could mean that the phase difference between the coupled oscillators do not change in time (ϕ
=
const), and the oscillations are periodic with constant amplitudes a
1,2. Thus, finding the conditions when these fixed points are stable would mean finding the conditions at which the synchronization occurs. Equations (36a) and (36b) show that the first approximation of the averaging method does not predict any nonzero steady-state amplitudes. The reason for this seems to lie in the fact that the even terms in (35) do not contribute to the value of integral (36a). To take account of them we have to employ in what follows the second approximation of the averaging method. As to the phase difference, ϕ, we see that stable is only antiphase steady-state regime with ϕ
=
π.LetR be the vector of the right-hand sides of (34). Also denote the operation of averaging by the angle brackets 〈
·
〉. Then 〈
R
〉 will mean the right-hand sides of (36a) and (36b). Following the conventional procedure we retain only vibrational terms in R,
(37)
R
~
=
R
-
〈
R
〉
,
and integrate R
~ up to an arbitrary function of the amplitudes chosen for simplicity to be equal to zero:
(38)
R
^
=
∫
‍
R
~
d
s
.
Now, upon calculation of the Jacobian matrix ∂
R
/
∂
(
a
1
,
a
2
,
φ
1
,
φ
2
)
T, we can write down the second approximation symbolically as
(39)
(
a
˙
1
,
a
˙
2
,
φ
˙
1
,
φ
˙
2
)
T
=
〈
R
〉
+
〈
∂
R
∂
(
a
1
,
a
2
,
φ
1
,
φ
2
)
T
R
^
〉
.
(The terms now neglected by averaging are of a higher order of magnitude with respect to the small parameters than the terms neglected in the first approximation.)Due to the symmetry of the case we may seta
1
=
a
2
=
a enabling us to present the resulting equations of the second approximation as compact as possible: (40a)
a
˙
=
-
1
2
a
ρ
cos
⁡
⁡
ϕ
+
μ
2
24
ω
3
a
(
ρ
(
ω
cos
⁡
⁡
ϕ
(
a
2
(
2
ω
2
+
5
)
+
12
ρ
(
γ
+
1
)
)
111111111111
+
2
ω
2
sin
⁡
ϕ
(
3
a
2
ρ
cos
⁡
⁡
ϕ
-
ρ
a
2
+
6
(
γ
+
1
)
)
111111111111
-
2
a
2
ω
(
ω
2
+
1
)
cos
⁡
⁡
2
ϕ
111111111111
-
6
(
γ
+
1
)
sin
⁡
ϕ
ω
3
a
2
(
2
ω
2
+
5
)
)
-
12
(
γ
+
1
)
ω
3
a
2
(
2
ω
2
+
5
)
)
,
(40b)
ϕ
˙
=
ρ
sin
⁡
ϕ
+
μ
2
24
ω
2
ρ
sin
⁡
ϕ
(
3
a
2
(
3
ω
2
-
10
)
+
8
a
2
(
ω
2
+
1
)
cos
⁡
⁡
ϕ
1111111111111
11
-
24
ρ
(
γ
+
1
)
(
ω
2
+
1
)
)
.The stable steady-state solution for the phase difference isϕ
=
π. Inserting this value into (40a) and putting a
˙
=
0 we get three different steady-state solutions for a (including the trivial one), of which only one,
(41)
a
=
2
3
ρ
ω
2
-
μ
2
(
γ
+
1
)
(
ρ
2
+
ω
2
)
μ
ρ
(
4
ω
2
+
7
)
,
being stable. Formula (41) indicates a soft transition to sustained oscillations as the coupling parameter ρ is progressively increased from a critical value. Thus we expect that the Hopf bifurcation is supercritical.Figure4 plots (
u
1
,
v
1
) phase planes for coupling strength below and above the Hopf bifurcation. When ϰ
1,2
=
0 (consumers are decoupled) the internal equilibrium is a stable focus (Figure 4(a)). For ϰ
1,2
=
0.011 (very weak coupling; just below the bifurcation) F
12 is still a stable focus, though a very gently winding one: the decay is slow (Figure 4(b)). For ϰ
1,2
=
0.012 (very weak coupling; just above the bifurcation) there is an unstable focus at F
12 and a stable oval limit cycle of small size representing low-amplitude quasi-harmonic oscillations (Figure 4(c)). On further increasing of coupling strength the limit cycle continuously grows in size and takes an irregular shape indicating nonlinearity of the synchronous oscillations (Figures 4(d)–4(f)).The development of a limit cycle in model (1a)–(1d) in relation to the coupling strength. The projections of phase trajectories and fixed points on (
u
1
,
v
1
) plane are presented. “•” and “∘” mark the respective stable and unstable internal steady state F
12; “×” stands for the boundary steady state F
1. (a) ϰ
1,2
=
0: CR pairs are decoupled; each has unique stable steady state. (b) ϰ
1,2
=
0.011: very weak coupling; the spiral winds only algebraically fast. The system is on the verge of the Hopf bifurcation. (c) ϰ
1,2
=
0.012: very weak coupling; the stable limit cycle is just born. (d) ϰ
1,2
=
0.025: weak coupling; nonlinear oscillations. (e) ϰ
1,2
=
0.05: moderate coupling; the limit cycle passes near the basin of the saddle point. (f) ϰ
1,2
=
0.6: moderate coupling; the limit cycle is about to merge with a heteroclinic cycle.
(a)
(b)
(c)
(d)
(e)
(f)As a practical matter, the range of very weak coupling not too far away from the Hopf bifurcation, where oscillations are quasi-linear and quasi-harmonic, is of less concern to us than is the range of far more feasible not-too-weak coupling, corresponding to well-developed substantially nonlinear oscillations. We are going to demonstrate that given conditions (20), system (1a)–(1d) exhibits relaxation oscillatory behavior, with the two coupled CR pairs being antiphase locked.By the assumption,0
<
ε
≪
1, meaning that system (1a)–(1d) is singularly perturbed. The slow variables are resources, u
1 and u
2, and the fast variables are consumers, v
1 and v
2. The standard practice of reducing such systems is the adiabatic elimination of the fast variables, when the left-hand side in the fast equation is replaced by zero, thus turning this differential equation into an algebraic equation. It is assumed that the fast variables quickly relax to their momentary equilibrium, quasi-steady-state, values obtained from the algebraic equations, in which the slow variables are treated as parameters. “Frozen” slow variables do not move substantially in this short adaptation time of the fast variables. Quasi-steady-state values of the fast variables can then be expressed by values of the slow variables. The fast variables hastily adapt to the motion of the slow variables. The former are entrained by the latter. The utility of quasi-steady-state approximation is that it allows us to reduce the dimension of the system by retaining only slow variables in the model. One has to establish the validity of the adiabatic elimination in each specific case by using the recommendations of the singular perturbation theory [22]. In particular, Tikhonov’s theorem [20] requires the quasi-steady state of the fast equations to be stable.To decompose the full system (1a)–(1d) into fast and slow subsystems, introduce fast time variable τ
=
t
/
ε. Now rescale (1a)–(1d) by replacing t with τ
ε and, after taking ε
=
0, it becomes(42a)
u
1
′
=
u
2
′
=
0
,
(42b)
v
1
′
=
(
u
1
-
δ
v
1
-
ϰ
2
v
2
)
v
1
,
(42c)
v
2
′
=
(
u
2
-
δ
v
2
-
ϰ
1
v
1
)
v
2
,where prime means differentiation with respect to τ. This is the fast subsystem, where u
1 and u
2 are replaced by their initial values and treated as parameters. It yields inner solution, valid for t
=
O
(
ε
).Settingε
=
0 in (1a)–(1d) leads to the slow subsystem(43a)
u
˙
1
=
γ
1
-
(
u
1
+
1
)
v
1
-
u
1
,
(43b)
u
˙
2
=
γ
2
-
(
u
2
+
1
)
v
2
-
u
2
,
(43c)
0
=
(
u
1
-
δ
v
1
-
ϰ
2
v
2
)
v
1
,
(43d)
0
=
(
u
2
-
δ
v
2
-
ϰ
1
v
1
)
v
2
,which produces outer solution, valid for t
=
O
(
1
). In this singular limit as ε
→
0, the subsystem defines a slow flow on the surface (slow manifold) given by (43c) and (43d). Outer solution is valid for those u
1 and u
2, for which the quasi-steady states of the fast subsystem are stable.We anticipate the dynamics of the full system (1a)–(1d) in its four-dimensional phase space (
u
1
,
u
2
,
v
1
,
v
2
) to consist of two typical motions: quickly approaching the slow manifold (43c) and (43d), and slowly sliding over it until a leave point (where the solution disappears) is reached. After that, the representing point may possibly jump to another local solution of (43c) and (43d).Thus, we ought to find all quasi-steady states of the fast subsystem (42a)–(42c), map the domains of their stability onto the slow phase plane (
u
1
,
u
2
), and then investigate the dynamics of the slow subsystem (43a)–(43d) with piecewise continuous functions.Fast subsystem (42a)–(42c), which is nothing but the conventional LVG model, has four quasisteady-states—three boundary and one interior—denoted by Q (the slow variables are deemed to be frozen):(44a)
Q
:
v
~
1
=
0
,
v
~
2
=
0
;
(44b)
Q
1
:
v
~
1
=
u
1
δ
,
v
~
2
=
0
;
(44c)
Q
2
:
v
~
1
=
0
,
v
~
2
=
u
2
δ
;
(44d)
Q
12
:
v
~
1
=
ϰ
2
u
2
-
δ
u
1
ϰ
1
ϰ
2
-
δ
2
,
v
~
2
=
ϰ
1
u
1
-
δ
u
2
ϰ
1
ϰ
2
-
δ
2
.Existence and stability of these quasi-steady states are determined by (
u
1
,
u
2
)—position of a representing point in the phase plane of the slow subsystem, shown in Figure 5(a).(a) Phase plane of the slow subsystem (43a)–(43d) sectored (by dashed lines) into stability domains of the corresponding quasi-steady states of the fast subsystem (42a)–(42c). Both boundary quasi-steady states, each corresponding to the situation when either of the two consumers completely dominates, are stable within the opening of the angle formed by the dashed lines. Lines u
1
=
γ
1
δ and u
2
=
γ
2 are the respective nullclines u
˙
1
=
0 and u
˙
2
=
0 of the piecewise subsystem (47a)-(47b). Lines u
1
=
γ
1 and u
2
=
γ
2
δ mean the same for the piecewise subsystem (52a)-(52b). Intersections of the nullclines (marked by the open circles) are equilibria of the associated piecewise slow subsystems, and they must lie outside the above-mentioned opening to allow for the relaxation oscillations. (b), (c), and (d) are the respective phase portraits of the fast subsystem (42a)–(42c) generated by points “b,” “c,” and “d” in the slow phase plane (a).
(a)
(b)
(c)
(d)Q, Q
1 and Q
2 always exist for all u
1 and u
2 from the positive quadrant of the slow phase plane. For not-too-weak coupling, such that ϰ
1
ϰ
2
>
δ
2, Q
12 exists for all u
1 and u
2 satisfying the condition δ
u
1
/
ϰ
2
<
u
2
<
u
1
ϰ
1
/
δ, that is, within the opening of the angle formed by lines δ
u
1
-
ϰ
2
u
2
=
0 and ϰ
1
u
1
-
δ
u
2
=
0 in Figure 5(a). The opening shrinks as the coupling strengths get weaker.Jacobian matrix of the fast subsystem,(45)
(
u
1
-
2
δ
v
~
1
-
ϰ
2
v
~
2
ε
0
-
ϰ
2
v
~
1
ε
-
ϰ
1
v
~
2
ε
u
2
-
2
δ
v
~
2
-
ϰ
1
v
~
1
ε
)
,
has the following sets of eigenvalues at (44a)–(44d):(46a)
λ
(
Q
)
:
u
1
ε
,
u
2
ε
;
(46b)
λ
(
Q
1
)
:
-
ϰ
1
u
1
ε
δ
+
u
2
ε
,
-
u
1
ε
;
(46c)
λ
(
Q
2
)
:
-
ϰ
2
u
2
ε
δ
+
u
1
ε
,
-
u
2
ε
;
(46d)
λ
(
Q
12
)
:
±
u
1
u
2
+
O
(
δ
)
ε
.Based on (46a)–(46d) one concludes that Q (the origin) is always an unstable node for all u
1 and u
2 from the positive quadrant of the slow phase plane. Q
1 is a stable node for δ
u
2
<
ϰ
1
u
1, that is, below the line ϰ
1
u
1
-
δ
u
2
=
0 in Figure 5(a), otherwise it is a saddle. Similarly, Q
2 is a stable node for δ
u
1
<
ϰ
2
u
2, that is, above the line δ
u
1
-
ϰ
2
u
2
=
0 in the plane of slow variables, otherwise it is a saddle. The interior quasi-steady state, Q
12, is always a saddle.The performed typology of fixed points of the fast subsystem (42a)–(42c) leads to three qualitatively different phase portraits depicted by Figures 5(b)–5(d).Suppose initiallyQ
1 is stable and Q
2 is not. Consumer 1 completely dominates. This corresponds to slow variables u
1 and u
2 being somewhere below the line δ
u
1
-
ϰ
2
u
2
=
0 of Figure 5(a). Fast subsystem (42a)–(42c) has phase portrait of a type shown in Figure 5(b). While u
1 remains much greater than γ
1
δ, the dynamics of the resources (treated as bifurcation parameters in reference to the consumers) is described by a system of two independent equations:(47a)
u
˙
1
=
γ
1
-
(
u
1
+
1
δ
+
1
)
u
1
,
(47b)
u
˙
2
=
γ
2
-
u
2
,which is a piecewise version of the slow subsystem (43a)–(43d) for v
1
=
u
1
/
δ and v
2
=
0. System (47a) and (47b) has stable steady state
(48)
u
^
1
(
1
)
=
1
2
(
r
1
-
δ
-
1
)
=
γ
1
δ
+
O
(
δ
2
)
,
u
^
2
(
1
)
=
γ
2
,
where r
1
=
1
+
δ
(
4
γ
1
+
2
+
δ
). This equilibrium lies in the upper left corner of Figure 5(a).While heading to (48), the trajectory crosses the line δ
u
1
-
ϰ
2
u
2
=
0 and enters the domain of bistability of both Q
1 and Q
2. Fast subsystem (42a)–(42c) takes new phase portrait of a type presented by Figure 5(c). However the dominance of consumer 1 persists.Upon introducing the deviationsξ
1 and η
2 from the respective steady states (48), the system (47a)-(47b) become(49a)
ξ
˙
1
=
-
ξ
1
(
r
1
+
ξ
1
)
δ
,
(49b)
η
˙
2
=
-
η
2
,whence one finds
(50)
ξ
1
(
t
)
=
r
1
(
r
1
/
ξ
1
(
0
)
+
1
)
exp
⁡
⁡
(
r
1
t
/
δ
)
-
1
,
η
2
(
t
)
=
η
2
(
0
)
exp
⁡
(
-
t
)
.It follows from (50), that the dynamics of variable u
1 is faster than that of u
2 due to small δ. Clearly, the representing point must have relaxed to the vertical line u
1
=
γ
1
δ well before approaching the horizontal line u
2
=
γ
2. Further developments depend on whether the involved individual CR systems are overdamped or underdamped.Ifδ is great enough to damp intrinsic oscillations in the constituent CR pairs, the representing point will slide along the nullcline u
˙
1
=
0 (slow manifold of the system (47a)-(47b)), steadily tending to point (48) (Figure 6(a)).The limit cycle of system (1a)–(1d) projected onto phase plane of the slow variables in the cases of (a) strong second-order damping, and (b) underdamping. The parameters of the model chosen for numerical simulation are ε
=
0.727273
×
1
0
-
3, γ
1
=
1.19375, γ
2
=
1, ϰ
1
=
0.5, ϰ
2
=
0.8, δ
=
0.1 (a), and δ
=
0.01 (b).
(a)
(b)Ifδ is small in terms of the condition (6), then the individual CR pairs are underdamped. In the immediate vicinity of u
^
1
(
1
) the division of variables on slow and fast ones loses its meaning and reduced (47a) is no longer valid. Instead of (47a) one should write down two equations:(51a)
u
˙
1
=
γ
1
-
(
u
1
+
1
)
v
1
-
u
1
,
(51b)
ε
v
˙
1
=
(
u
1
-
δ
v
1
)
v
1
.However this system is identical to (3a) and (3b) for an uncoupled CR system and, in essence, describes convergence to a focal point via dying oscillations. In the plane of slow variables (
u
1
,
u
2
) these oscillations manifest themselves in damped transverse fluctuations superimposed on the independent vertical motion along the nullcline u
˙
1
=
0 (u
1
=
u
^
1
(
1
)
≈
γ
1
δ) toward point (48) (Figure 6(b)).By virtue of the condition (20), the trajectory has to cross the line ϰ
1
u
1
-
δ
u
2
=
0 on its way toward the neighborhood of steady state (48). As soon as this has happened, node Q
1 in the plane (
v
1
,
v
2
) will be absorbed by saddle Q
12. A new phase portrait of the fast subsystem (42a)–(42c) takes on the appearance of Figure 5(d). Consumer 1 rapidly washes out, and the alternative boundary quasi-steady state Q
2 becomes stable, with consumer 2 dominating.In terms of the four-dimensional phase space of full system (1a)–(1d), the representing point is now in the other stable branch of the slow manifold (43c) and (43d). The motion over this alternative branch obeys the piecewise subsystem:(52a)
u
˙
1
=
γ
1
-
u
1
,
(52b)
u
˙
2
=
γ
2
-
(
u
2
+
1
δ
+
1
)
u
2
,with the initial conditions u
1
(
0
)
=
u
^
1
(
1
)
≈
γ
1
δ and u
2
(
0
)
=
γ
1
ϰ
1.The dynamics of (52a) and (52b) is basically similar to that of (47a) and (47b) analyzed above. System (52a) and (52b) has a stable steady state:
(53)
u
^
1
(
2
)
=
γ
1
,
u
^
2
(
2
)
=
1
2
(
r
2
-
δ
-
1
)
=
γ
2
δ
+
O
(
δ
2
)
,
where r
2
=
1
+
δ
(
4
γ
2
+
2
+
δ
). This equilibrium lies in the lower right corner of Figure 5(a).Variableu
2, being more rapid in comparison to u
1, quickly enters the neighborhood of the nullcline u
˙
2
=
0 given by u
2
=
u
^
2
(
2
)
≈
γ
2
δ and then—depending on the value of δ—finally approaches the nullcline either monotonically (Figure 6(a)) or via damped oscillations according to equations(54a)
u
˙
2
=
γ
2
-
(
u
2
+
1
)
v
2
-
u
2
,
(54b)
ε
v
˙
2
=
(
u
2
-
δ
v
2
)
v
2
.(Figure 6(b)). System (54a)-(54b) describes underdamped intrinsic oscillations of uncoupled consumer 2 for small δ.At the same time,u
1 steadily and independently tends to u
^
1
(
2
)
=
γ
1. Again, because point (53) is located below the line δ
u
1
-
ϰ
2
u
2
=
0 (on the strengths of the condition (20)), the trajectory would certainly cross that line at a point (
γ
2
ϰ
2
,
γ
2
δ
), whereupon node Q
2 would be absorbed by saddle Q
12. The system returns to the first branch of the slow manifold, and thereby the oscillatory cycle gets closed.


## 3.1. A Single CR Pair

Whenϰ
1,2
=
0, the communities are uncoupled and completely independent. An isolated CR pair is governed by equations:(2a)
u
˙
=
γ
-
(
u
+
1
)
v
-
u
,
(2b)
ε
v
˙
=
(
u
-
δ
v
)
v
.There exist two nonnegative steady states:(3a)
u
¯
=
γ
,
v
¯
=
0
;
(3b)
u
¯
=
1
2
(
1
+
(
4
γ
+
2
+
δ
)
δ
-
1
-
δ
)
=
γ
δ
+
O
(
δ
2
)
,
v
¯
=
1
2
δ
(
1
+
(
4
γ
+
2
+
δ
)
δ
-
1
-
δ
)
=
γ
-
γ
(
γ
+
1
)
δ
+
O
(
δ
2
)
.The linearization of (2a) and (2b) takes the form
(4)
J
=
(
-
v
¯
-
1
-
u
¯
-
1
v
¯
ε
u
¯
-
2
δ
v
¯
ε
)
,
where (u
¯, v
¯) is one of the above steady states (3a) and (3b).At (3a), one eigenvalue is negative and one is positive: λ
1
=
-
1, λ
2
=
γ
/
ε. Thus (3a) is a saddle point. At (3b), Tr
⁡
J
=
-
(
δ
/
ε
+
1
)
v
¯
-
1
<
0 and det
⁡
J
=
(
1
+
2
δ
v
¯
+
δ
)
v
¯
/
ε
>
0, so the steady state is a stable node/focus. Specifically, focus is the case for
(5)
(
Tr
⁡
J
)
2
-
4
det
⁡
J
=
(
γ
2
δ
2
+
O
(
δ
3
)
)
ε
-
2
+
(
-
4
γ
+
O
(
δ
)
)
ε
-
1
+
O
(
1
)
<
0
,
whence one obtains an asymptotic estimate
(6)
δ
=
o
(
ε
1
/
2
)
.
Figure 1(c) illustrates this focus-node bifurcation numerically. Damped oscillations is a well-known inherent feature of the photon-carrier dynamics in class-B lasers.Phase portraits of an isolated CR system (2a) and (2b) for two cases when the unique stable equilibrium (3b) is a node (a), and a focus (b). Shown are vertical (u
˙
=
0) and horizontal (v
˙
=
0) nullclines. Steady states are marked by open circles. Numerical simulations have been carried out with data from Table 2. The respective values of the second order loss parameter δ in (a) and (b) are 0.1 and 0.01. Typically, consumer scale (v) is enormous; the equilibrium point happens to lie very close to the origin. For better appearance, henceforward we display all the graphs using coordinate transformation u
→
arsinh
(
u
/
γ
δ
) and v
→
arsinh
(
v
/
γ
). (c) Eigenvalues of (4) at (3b) as a function of the second order loss δ. Critical value δ
cr indicates a boundary between the regions of damped oscillations and aperiodic damping.
(a)
(b)
(c)An isolated CR system (2a) and (2b) admits, therefore, the only solution which tends asymptotically towards the unique steady state. Periodic solutions are excluded. However the temporal dynamics of approaching this steady state essentially depends on the interplay between ε and δ, and is worth another look. There are in fact three timescales, O
(
ε
), O
(
1
), and O
(
δ
), involved in the CR system (2a) and (2b) when it is overdamped, and two—O
(
ε
) and O
(
1
)—when underdamped.System (2a) and (2b) is singularly perturbed because the derivative of one of its state variables, v, is multiplied by a small positive parameter ε. Singular perturbation cause two-time-scale behavior of the system characterized by the presence of slow and fast transients in the system’s response to external stimuli.Replacingt in (2a) and (2b) with a fast time variable τ
=
t
/
ε and setting ε
=
0 we obtain the fast subsystem:(7a)
u
′
=
0
,
(7b)
v
′
=
(
u
-
δ
v
)
v
,where prime means differentiation with respect to τ. In the stretched timescale τ the slow resource variable u, according to (7a), is replaced by its initial value and reckoned as constant parameter. Equation (7b) is of a logistic type which has the solution
(8)
v
(
τ
)
=
u
/
δ
1
+
(
u
/
δ
v
(
0
)
-
1
)
exp
⁡
⁡
(
-
u
τ
)
,
valid for t
=
O
(
ε
).After a lapse of considerable time (in the fast scale)v converges to either of two fixed points depending on a sign of u:
(9)
lim
⁡
τ
→
∞
v
(
τ
)
=
{
u
δ
,
if
u
>
0
;
0
,
if
u
<
0
.
It means that every single trajectory starting within the positive quadrant of (
u
,
v
) plane far enough from the stable steady state (3b) will run almost parallel to the vertical axis and hit the line u
-
δ
v
=
0 practically in a finite time of order O
(
ε
). This is shown in Figures 1(a) and 1(b).Now setε
=
0 in (2a) and (2b) to get the slow subsystem:(10a)
u
˙
=
γ
-
(
u
+
1
)
v
-
u
,
(10b)
0
=
(
u
-
δ
v
)
v
.Equation (10b) describes a slow manifold consisting of two lines in the (
u
,
v
) plane: v
=
u
/
δ and v
=
0. By (9), the former attracts all the trajectories in the first quadrant, while the latter—all those in the second. In as much as the quasi-steady state of (7b), v
¯
=
u
/
δ, is an isolated root of (10b) and v
¯ is a stable solution of (10b) for any u
>
0; the assumptions of Tikhonov’s theorem [20] are satisfied, and one may proceed to approximate u and v in terms of the solution of the reduced system (10a) and (10b) in the slow timescale. It means that after arriving at the slow manifold v
=
u
/
δ, the representing point of the full system (2a) and (2b) will move along the manifold toward the equilibrium point with a characteristic velocity of order O
(
1
).In the immediate proximity to equilibrium (3b) the behavior of the trajectory is determined by the type of the fixed point, whether it is a node or a focus. Eventually this depends on the value of δ. Namely, close to a stable node, the system has two distinct real negative eigenvalues, one fast (λ
1), and one slow (λ
2):
(11)
λ
1
=
(
-
γ
δ
+
O
(
δ
2
)
)
ε
-
1
+
(
δ
-
1
+
O
(
1
)
)
+
O
(
ε
)
,
λ
2
=
(
-
δ
-
1
+
O
(
1
)
)
+
O
(
ε
)
,
for
1
>
δ
≫
ε
.
Since |
λ
1
|
≫
|
λ
2
|, trajectories starting off the associated eigenvector 2 (which is tangent to the slow manifold v
=
u
/
δ) converge to that vector along lines almost parallel to eigenvector 1 (which is parallel to the vertical axis v). As they approach vector 2 they become tangent to it and move along it up to the very nodal point (Figure 1(a)). The characteristic time constant of this final stage is of order O
(
δ
).In the vicinity of a stable focus the motion is qualitatively different: trajectories still keep converging to the equilibrium, but no longer follow the slow manifoldv
=
u
/
δ (Figure 1(b)). The reason is that for a focus eigenvectors are complex. Recalling that condition (6) must be true for a focus, calculate the eigenvalues in a limit case of δ
→
0:
(12)
λ
1,2
=
-
1
2
(
γ
+
1
)
±
i
γ
ε
for
δ
=
0
.
Loosely speaking, one may think that near a focal point separation of state variables into slow and fast ones ceases to have its conventional meaning. There is no point to talk about motion along the slow manifold since there are no reduced one-dimensional systems corresponding to the neighborhood of a focus. The perturbation parameter ε affects only the frequency of damped oscillation, but not the damping rate. The radius of the focal spiral uniformly shrinks with a time constant of order O
(
1
).ToO
(
1
) for small δ there is a way to recast (2a) and (2b) near the focal equilibrium (
0
,
γ
) in a convenient form where the perturbation parameter ε does not multiply any right hand side. Namely, performing the scaling t
=
μ
2
s, u
=
μ
2
γ
ξ, and v
=
γ
(
1
+
η
), where μ
2
=
ε
/
γ, one obtains
(13)
ξ
˙
=
-
η
-
μ
2
ξ
(
1
+
γ
(
1
+
η
)
)
,
η
˙
=
ξ
(
1
+
η
)
.
Here dots stand for differentiation with respect to the time variable s.Equations (13) represent a weakly perturbed Hamiltonian system. For μ
=
0 the system is pure Hamiltonian and admits the first integral
(14)
H
=
1
2
ξ
2
+
η
-
ln
⁡
⁡
(
1
+
η
)
,
which is a conserved quantity (H
˙
=
0). The periodic solutions of the Hamiltonian system form a one-parameter family with the equilibrium (
0,0
) as center point. The condition 0
<
μ
≪
1 makes (13) quasi-conservative with phase trajectories slowly spiralling to the equilibrium (Figure 2).Figure 2
Phase-plane trajectory of the weakly perturbed Hamiltonian system (13) for γ
=
1.19375, ε
=
0.727273
×
1
0
-
3, and μ
=
0.157107. Near the equilibrium the oscillations are almost harmonic.System (13) can be rewritten as a second order differential equation for η only:
(15)
η
¨
+
η
+
η
2
-
η
˙
2
1
+
η
+
μ
2
(
1
+
γ
(
1
+
η
)
)
η
˙
=
0
.
This is the equation for a nonlinear quasi-conservative oscillator. As is seen from (15), for |
η
|
≪
1 and |
η
˙
|
≪
1 the oscillator is not only quasi-conservative, but also quasi-linear with the frequency ω
0
=
1. These conditions motivate introducing the new variable ζ
=
η
/
μ. As a result one obtains
(16)
ζ
¨
+
ζ
=
μ
(
-
ζ
2
+
ζ
˙
2
1
+
μ
ζ
-
μ
(
1
+
γ
+
μ
γ
ζ
)
ζ
˙
)
.
This equation of quasi-conservative and quasi-linear oscillator will be used later on when analyzing the onset of synchronous periodicity.


## 3.2. Coupled Communities: Stability

In the following, we will distinguish two main types of coupling: strong,ϰ
1,2
>
1, and weak, ϰ
1,2
<
1.Physically feasible equilibria, or fixed points,(
u
¯
1
,
u
¯
2
,
v
¯
1
,
v
¯
2
), of (1a)–(1d) are those for which v
¯
1,2
⩾
0. We denote the interior fixed point by F
12
=
(
u
¯
1
,
u
¯
2
,
v
¯
1
,
v
¯
2
), where the subscripts stand for the consumers. Lack of a certain index at a boundary fixed point means that the consumer concerned is not present (extinct). Thus F
1
=
(
u
¯
1
,
u
¯
2
,
v
¯
1
,
0
) and F
2
=
(
u
¯
1
,
u
¯
2
,
0
,
v
¯
2
) designate either of one-consumer equilibria corresponding to dominance, while F
=
(
u
¯
1
,
u
¯
2
,
0,0
) means both consumers having been washed out.Model (1a)–(1d) has four feasible steady states. To O
(
1
) for small δ
(17a)
F
:
u
¯
1
=
γ
1
,
u
¯
2
=
γ
2
,
v
¯
1
=
0
,
v
¯
2
=
0
;
(17b)
F
1
:
u
¯
1
=
0
,
u
¯
2
=
γ
2
,
v
¯
1
=
γ
1
,
v
¯
2
=
0
;
(17c)
F
2
:
u
¯
1
=
γ
1
,
u
¯
2
=
0
,
v
¯
1
=
0
,
v
¯
2
=
γ
2
;
F
12
:
u
¯
1
=
ϰ
1
γ
1
-
ϰ
2
γ
2
-
ϰ
1
ϰ
2
+
1
±
R
2
(
ϰ
1
-
1
)
,
u
¯
2
=
-
ϰ
1
γ
1
+
ϰ
2
γ
2
-
ϰ
1
ϰ
2
+
1
±
R
2
(
ϰ
2
-
1
)
,
(17d)
v
¯
1
=
u
¯
2
ϰ
1
,
v
¯
2
=
u
¯
1
ϰ
2
,where
(18)
R
=
(
ϰ
1
γ
1
-
ϰ
2
γ
2
-
ϰ
1
ϰ
2
+
1
)
2
+
4
ϰ
2
(
ϰ
1
-
1
)
(
ϰ
1
γ
1
-
γ
2
)
.Boundary equilibriaF, F
1, and F
2 always exist. F
12 exists if
(19)
1
ϰ
2
<
γ
2
γ
1
<
ϰ
1
for
ϰ
1
>
1
∧
ϰ
2
>
1
(
strong
coupling
)
,
or
(20)
ϰ
1
<
γ
2
γ
1
<
1
ϰ
2
for
ϰ
1
<
1
∧
ϰ
2
<
1
(
weak
coupling
)
.
In case (19) the expression with “+” in (17d) is realized, while in case (20) the one with “−” is feasible.F is always unstable, because two of the four associated eigenvalues are positive:
(21)
λ
(
F
)
:
γ
1
ε
,
γ
2
ε
,
-
1
,
-
1
.Correct toO
(
1
) in ε, the eigenvalues for F
1 and F
2 are(22a)
λ
(
F
1
)
:
-
1
,
-
1
2
(
γ
1
+
1
)
±
i
γ
1
ε
,
γ
2
-
ϰ
1
γ
1
ε
;
(22b)
λ
(
F
2
)
:
1
,
-
1
2
(
γ
2
+
1
)
±
i
γ
2
ε
,
γ
1
-
ϰ
2
γ
2
ε
,whence it follows that (17b) and (17c) are stable when(23a)
γ
2
γ
1
<
ϰ
1
,
(23b)
γ
1
γ
2
<
ϰ
2
,respectively.The necessary and sufficient conditions for all the eigenvalues of the Jacobian matrix,(24)
(
-
v
¯
1
-
1
0
-
ϰ
2
v
¯
2
-
1
0
0
-
v
¯
2
-
1
0
-
ϰ
1
v
¯
1
-
1
v
¯
1
ε
0
0
-
ϰ
2
v
¯
1
ε
0
v
¯
2
ε
-
ϰ
1
v
¯
2
ε
0
)
,
evaluated at F
12, to have negative real parts are, from the Routh-Hurwitz criterion, (25a)
c
0
>
0
,
(25b)
c
3
>
0
,
(25c)
c
2
c
3
-
c
1
>
0
,
(25d)
c
1
(
c
2
c
3
-
c
1
)
-
c
0
c
3
2
>
0
,where c
0, c
1, c
2, and c
3 are the coefficients of the characteristic polynomial λ
4
+
c
3
λ
3
+
c
2
λ
2
+
c
1
λ
+
c
0 of (24):
(26)
c
0
=
v
¯
1
v
¯
2
(
ϰ
1
(
v
¯
1
-
ϰ
2
(
v
¯
1
+
v
¯
2
+
1
)
)
+
ϰ
2
v
¯
2
+
1
)
ε
-
2
,
1
c
1
=
-
ϰ
1
ϰ
2
v
¯
1
v
¯
2
(
v
¯
1
+
v
¯
2
+
2
)
ε
-
2
+
(
v
¯
1
(
v
¯
2
(
ϰ
1

+

ϰ
2
+
ϰ
1
v
¯
1

+

ϰ
2
v
¯
2

+

2
)
+
1
)
+
v
¯
2
)
ε
-
1
,
c
2
=
-
ϰ
1
ϰ
2
v
¯
1
v
¯
2
ε
-
2
+
(
v
¯
1
(
v
¯
2
(
ϰ
1
+
ϰ
2
)
+
1
)
+
v
¯
2
)
ε
-
1
+
(
v
¯
1
+
1
)
(
v
¯
2
+
1
)
,
c
3
=
v
¯
1
+
v
¯
2
+
2
.
11111111111111111111111111111To analyze the validity of (25a)–(25d) we put ϰ
2
=
μ
ϰ
1, where μ
=
O
(
1
). Then for (25a) to be true, ϰ
1
=
o
(
1
) should be met, which is incompatible with strong coupling, yet possible for weak coupling. Conditions (25b) and (25c) are always valid. As is known, (25d) guarantees a simple complex conjugate pair of eigenvalues corresponding to a linearization about steady state F
12 to have negative real part. For ϰ
1
=
o
(
1
) it boils down to
(27)
(
-
v
¯
1
v
¯
2
(
v
¯
1
+
v
¯
2
+
2
)
(
v
¯
1
+
v
¯
2
+
v
¯
1
2
+
v
¯
2
2
)
μ
ϰ
1
2
+
O
(
ϰ
1
3
)
)
ε
-
3
+
(
(
v
¯
1
+
1
)
(
v
¯
2
+
1
)
(
v
¯
1
-
v
¯
2
)
2
+
O
(
ϰ
1
)
)
ε
-
2
+
O
(
ε
-
1
)
>
0
,
whence it follows that
(28)
ϰ
1
=
o
(
ε
1
/
2
)
.
With regard to a fairly small value of ε, (28) may be thought to be broken under most physically meaningful conditions unless coupling is infinitesimally weak. Hence, normally, condition (25d) of the Routh-Hurwitz criterion is never fulfilled and the interior fixed point F
12—if it exists—is always unstable by growing oscillations.Existence and stability conditions of possible nonnegative equilibrium points are summarized in Table1.Table 1
Existence and stability conditions of nonnegative equilibria of system (1a)–(1d).
Equilibrium
Existence
Stability
F
Always
Unstable
F
1
Always
γ
2
/
γ
1
<
ϰ
1
F
2
Always
γ
2
/
γ
1
>
1
/
ϰ
2
F
12
1
/
ϰ
2
<
γ
2
/
γ
1
<
ϰ
1 for ϰ
1,2
>
1 (strong coupling)
Unstable
ϰ
1
<
γ
2
/
γ
1
<
1
/
ϰ
2 for ϰ
1,2
<
1 (weak coupling)
ϰ
1,2
=
o
(
ε
1
/
2
)


## 3.3. Strong Coupling: Bistability and Hysteresis

Strong coupling in (1a)–(1d) makes bistability of boundary equilibria possible, as evident from Table 1. When both F
1 and F
2 are stable with an unstable coexistence steady state F
12, the system being studied is able to exhibit a hysteresis effect. Given the strong coupling, suppose that the inflow of resource 2 is kept at some constant level, γ
2
=
γ
2
*, while the inflow of resource 1, γ
1, steadily increases from a value less than γ
2
*
/
ϰ
1 along the line γ
2
=
γ
2
* in the (
γ
1
,
γ
2
) parameter plane, as shown in Figure 3(a). Referring to (17a)–(17d) and (23a) and (23b), one sees that initially F
1 is unstable; F
2 is stable and F
12 does not exist with the complete dominance of consumer 2. Within the interval γ
2
*
/
ϰ
1
<
γ
1
<
ϰ
2
γ
2
* steady state F
1 becomes stable yet empty and F
12 becomes existent (according to (19)) yet unstable, so the situation remains unchanged until point (
ϰ
2
γ
2
*
,
γ
2
*
) in Figure 3(a) has been reached from the left. For a larger γ
1, state F
2 gives up its stability and the system jumps to F
1. Consumer 2 gets washed out, while consumer 1 takes over. If now we start reducing γ
1, the system remains in F
1 until γ
1 drops to the lower critical value γ
2
*
/
ϰ
1, beyond which F
1 is no longer stable and there is a reverse jump to F
2. In other words, as γ
1 progresses along γ
2
=
γ
2
* there is a discontinuous switch from consumer 2 to consumer 1 at ϰ
2
γ
2
*, while as γ
1 retraces its steps, there is a discontinuous switch from consumer 1 to consumer 2 at γ
2
*
/
ϰ
1. Figure 3(b) illustrates how the steady-state level of consumer 1 responds to infinitesimally slow changes of the inflow rate of its own resource. The hysteresis is made possible thanks to the concurrent stability of both boundary equilibria and instability of the interior fixed point for γ
1
∈
(
γ
2
*
/
ϰ
1
,
ϰ
2
γ
2
*
). In terms of electronics, such a situation would describe a flip-flop circuit—bistable multivibrator—having two stable conditions, each corresponding to one of two alternative input signals.Hysteresis in the two-consumer, two-resource system (1a)–(1d). It takes place in the case of strong coupling ϰ
1,2
>
1. (a) The resource-supply parametric plane (
γ
1
,
γ
2
). Consumer 1 completely dominates below the line γ
1
-
ϰ
2
γ
2
=
0, whereas consumer 2—above the line ϰ
1
γ
1
-
γ
2
=
0. Both boundary steady states are stable in the region of bistability confined by the two aforementioned lines and marked off by gray, with the dominance of either consumer being a matter of path-dependency. (b) Equilibrium level of consumer 1 as a function of the resource supply. Numerical values of the coupling strengths are chosen to be ϰ
1
=
1.2 and ϰ
2
=
1.5.
(a)
(b)


## 3.4. Weak Coupling: Antiphase Relaxation Oscillations

As seen from Table1, the very existence of interior equilibrium F
12 in a case of weak coupling (condition (20)) implies instability of both boundary fixed points, F
1 and F
2. System (1a)–(1d) happens to possess four nonnegative steady states, none of them being stable. As we have found, F
12 is unstable through growing oscillations. In such a case, the model would thus be expected to have a limit cycle in its four-dimensional phase space corresponding to self-sustained oscillations.As is known, Hopf bifurcations come in both super- and subcritical types. If a small, attracting limit cycle appears immediately after the fixed point goes unstable, and if its amplitude shrinks back to zero as the parameter is reversed, the bifurcation is supercritical; otherwise, it’s probably subcritical, in which case the nearest attractor might be far from the fixed point, and the system may exhibit hysteresis as the parameter is reversed.To check whether the bifurcation is supercritical or subcritical, we employ the averaging method [21]. For negligible δ the problem reduces to weakly coupled quasi-conservative oscillators. Introduce the new time s
=
t
/
μ
2, where μ
2
=
ε
/
γ
1. The new dynamic variables ξ
1,2 and η
1,2 are defined by formulas ξ
1,2
=
u
1,2
/
μ
2
γ
1,2 and η
1,2
=
(
v
1,2
-
γ
1,2
)
/
γ
1,2. Thus the variables ξ
1,2 and η
1,2 are the respective deviations of u
1,2 and v
1,2 from their standalone steady-state values. With these new variables equations (1a)–(1d) become(29a)
ξ
˙
1
=
-
η
1
-
μ
2
ξ
1
(
1
+
γ
1
(
1
+
η
1
)
)
,
(29b)
η
˙
1
=
ξ
1
(
1
+
η
1
)
-
(
γ
2
γ
1
)
ρ
2
(
1
+
η
1
)
(
1
+
η
2
)
,
(29c)
ξ
˙
2
=
-
η
2
-
μ
2
ξ
2
(
1
+
γ
2
(
1
+
η
2
)
)
,
(29d)
η
˙
2
=
(
γ
2
γ
1
)
ξ
2
(
1
+
η
2
)
-
ρ
1
(
1
+
η
1
)
(
1
+
η
2
)
,where dots mean differentiation with respect to s and ρ
1,2
=
ϰ
1,2
/
μ
2 are the scaled coupling strengths.For brevity, we restrict our consideration to the case of identical CR pairs and symmetric coupling settingγ
1,2
=
γ and ρ
1,2
=
ρ. With these assumptions system (29a)–(29d) can be transformed to two coupled second-order equations for quasi-conservative, quasi-linear oscillators. This can be done as follows. First, (29b) and (29d) are solved for ξ
1 and ξ
2. Secondly, (29b) and (29d) are differentiated with respect to time to get η
¨
1 and η
¨
2. Thirdly, ξ
1, ξ
2, ξ
˙
1, and ξ
˙
2 are plugged into the equations for η
¨
1 and η
¨
2. And finally, η
1,2 are replaced by the new variables ζ
1,2: η
1,2
=
μ
ζ
1,2.As a consequence, we arrive at the following coupled equations:(30)
ζ
¨
1,2
+
ζ
1,2
=
-
ρ
ζ
˙
2,1
+
μ
Z
1
(
1,2
)
+
μ
2
Z
2
(
1,2
)
+
O
(
μ
3
)
,
where
(31)
Z
1
(
1,2
)
=
-
ρ
(
γ
+
1
)
-
ζ
1,2
2
+
ζ
˙
1,2
2
-
ρ
ζ
1,2
ζ
˙
2,1
,
Z
2
(
1,2
)
=
-
(
γ
+
1
)
(
ζ
˙
1,2
+
ρ
ζ
2,1
)
-
ζ
1
(
ρ
(
2
γ
+
1
)
+
ζ
˙
1,2
2
)
.We seek the solution of (30) in the quasi-harmonic form
(32)
ζ
1,2
=
a
1,2
cos
⁡
⁡
(
ω
s
+
φ
1,2
)
,
ζ
˙
1,2
=
-
a
1,2
ω
sin
⁡
(
ω
s
+
φ
1,2
)
,
where a
1,2 and φ
1,2 are slowly varying amplitudes and phases, and ω is the frequency of the synchronous oscillations. Differentiating the assumed form of ζ
1,2 and equating the result to the assumed form of ζ
˙
1,2 yields the first pair of relationships between a
1,2 and φ
1,2:
(33)
a
˙
1,2
cos
⁡
⁡
(
ω
s
+
φ
1,2
)
-
a
1,2
φ
˙
1,2
sin
⁡
(
ω
s
+
φ
1,2
)
=
0
.Then, differentiatingζ
˙
1,2 and substituting the resulting expression for ζ
¨
1,2 as well as the assumed forms for ζ
1,2 and ζ
˙
1,2 into (30) yields the second pair of equations relating a
1,2 and φ
1,2. Separating into equations for the rate of change of a
1,2 and φ
1,2 one obtains
(34)
a
˙
1,2
=
ω
-
1
sin
⁡
(
ω
s
+
φ
1,2
)
(
A
0
(
1,2
)
+
μ
A
1
(
1,2
)
+
μ
2
A
2
(
1,2
)
)
,
φ
˙
1,2
=
a
1,2
-
1
ω
-
1
cos
⁡
⁡
(
ω
s
+
φ
1,2
)
(
A
0
(
1,2
)
+
μ
A
1
(
1,2
)
+
μ
2
A
2
(
1,2
)
)
,
where
(35)
A
0
(
1,2
)
=
-
a
1,2
(
ω
2
-
1
)
cos
⁡
⁡
(
ω
s
+
φ
1,2
)
-
ρ
a
2,1
ω
sin
⁡
(
ω
s
+
φ
2,1
)
,
A
1
(
1,2
)
=
a
1,2
(
-
a
1,2
ω
2
sin
2
(
ω
s
+
φ
1,2
)
+
a
1,2
cos
⁡
2
(
ω
s
+
φ
1,2
)
-
ρ
a
2,1
ω
cos
⁡
⁡
(
ω
s
+
φ
1,2
)
sin
⁡
(
ω
s
+
φ
2,1
)
-
a
1,2
ω
2
sin
2
)
+
ρ
(
γ
+
1
)
,
A
2
(
1,2
)
=
a
1,2
cos
⁡
(
ω
s
+
φ
1,2
)
×
(
a
1,2
2
ω
2
sin
2
(
ω
s
+
φ
1,2
)
+
ρ
(
2
γ
+
1
)
)
+
(
γ
+
1
)
×
(
ρ
a
2,1
cos
⁡
⁡
(
ω
s
+
φ
2,1
)
-
a
1,2
ω
sin
⁡
(
ω
s
+
φ
1,2
)
)
.To this point no approximations have been made except of the expansion in powers ofμ in (30). Averaging equations (34) over the period 2
π
/
ω and considering a
1,2, φ
1,2, a
˙
1,2 and φ
˙
1,2 to be constants while performing the averaging, one obtains the following equations describing the slow variations of a
1,2 and φ
1,2: (36a)
a
˙
1,2
=
1
2
π
∫
0
2
π
/
ω
‍
(
A
0
(
1,2
)
+
μ
A
1
(
1,2
)
+
μ
2
A
2
(
1,2
)
)
×
sin
⁡
(
ω
s
+
φ
1,2
)
d
s
=
-
ρ
a
2,1
ω
cos
⁡
⁡
ϕ
+
μ
2
(
γ
+
1
)
(
a
1,2
ω
-
ρ
a
2,1
sin
⁡
ϕ
)
2
ω
,
(36b)
φ
˙
1,2
=
1
2
π
a
1,2
∫
0
2
π
/
ω
‍
(
A
0
(
1,2
)
+
μ
A
1
(
1,2
)
+
μ
2
A
2
(
1,2
)
)
×
cos
⁡
⁡
(
ω
s
+
φ
1,2
)
d
s
=
ρ
(
a
1,2
2
+
a
2,1
2
)
sin
⁡
ϕ
2
a
1,2
a
2,1
+
O
(
μ
2
)
,where ϕ
=
φ
1,2
-
φ
2,1.Any stable fixed points of (36a) and (36b) could mean that the phase difference between the coupled oscillators do not change in time (ϕ
=
const), and the oscillations are periodic with constant amplitudes a
1,2. Thus, finding the conditions when these fixed points are stable would mean finding the conditions at which the synchronization occurs. Equations (36a) and (36b) show that the first approximation of the averaging method does not predict any nonzero steady-state amplitudes. The reason for this seems to lie in the fact that the even terms in (35) do not contribute to the value of integral (36a). To take account of them we have to employ in what follows the second approximation of the averaging method. As to the phase difference, ϕ, we see that stable is only antiphase steady-state regime with ϕ
=
π.LetR be the vector of the right-hand sides of (34). Also denote the operation of averaging by the angle brackets 〈
·
〉. Then 〈
R
〉 will mean the right-hand sides of (36a) and (36b). Following the conventional procedure we retain only vibrational terms in R,
(37)
R
~
=
R
-
〈
R
〉
,
and integrate R
~ up to an arbitrary function of the amplitudes chosen for simplicity to be equal to zero:
(38)
R
^
=
∫
‍
R
~
d
s
.
Now, upon calculation of the Jacobian matrix ∂
R
/
∂
(
a
1
,
a
2
,
φ
1
,
φ
2
)
T, we can write down the second approximation symbolically as
(39)
(
a
˙
1
,
a
˙
2
,
φ
˙
1
,
φ
˙
2
)
T
=
〈
R
〉
+
〈
∂
R
∂
(
a
1
,
a
2
,
φ
1
,
φ
2
)
T
R
^
〉
.
(The terms now neglected by averaging are of a higher order of magnitude with respect to the small parameters than the terms neglected in the first approximation.)Due to the symmetry of the case we may seta
1
=
a
2
=
a enabling us to present the resulting equations of the second approximation as compact as possible: (40a)
a
˙
=
-
1
2
a
ρ
cos
⁡
⁡
ϕ
+
μ
2
24
ω
3
a
(
ρ
(
ω
cos
⁡
⁡
ϕ
(
a
2
(
2
ω
2
+
5
)
+
12
ρ
(
γ
+
1
)
)
111111111111
+
2
ω
2
sin
⁡
ϕ
(
3
a
2
ρ
cos
⁡
⁡
ϕ
-
ρ
a
2
+
6
(
γ
+
1
)
)
111111111111
-
2
a
2
ω
(
ω
2
+
1
)
cos
⁡
⁡
2
ϕ
111111111111
-
6
(
γ
+
1
)
sin
⁡
ϕ
ω
3
a
2
(
2
ω
2
+
5
)
)
-
12
(
γ
+
1
)
ω
3
a
2
(
2
ω
2
+
5
)
)
,
(40b)
ϕ
˙
=
ρ
sin
⁡
ϕ
+
μ
2
24
ω
2
ρ
sin
⁡
ϕ
(
3
a
2
(
3
ω
2
-
10
)
+
8
a
2
(
ω
2
+
1
)
cos
⁡
⁡
ϕ
1111111111111
11
-
24
ρ
(
γ
+
1
)
(
ω
2
+
1
)
)
.The stable steady-state solution for the phase difference isϕ
=
π. Inserting this value into (40a) and putting a
˙
=
0 we get three different steady-state solutions for a (including the trivial one), of which only one,
(41)
a
=
2
3
ρ
ω
2
-
μ
2
(
γ
+
1
)
(
ρ
2
+
ω
2
)
μ
ρ
(
4
ω
2
+
7
)
,
being stable. Formula (41) indicates a soft transition to sustained oscillations as the coupling parameter ρ is progressively increased from a critical value. Thus we expect that the Hopf bifurcation is supercritical.Figure4 plots (
u
1
,
v
1
) phase planes for coupling strength below and above the Hopf bifurcation. When ϰ
1,2
=
0 (consumers are decoupled) the internal equilibrium is a stable focus (Figure 4(a)). For ϰ
1,2
=
0.011 (very weak coupling; just below the bifurcation) F
12 is still a stable focus, though a very gently winding one: the decay is slow (Figure 4(b)). For ϰ
1,2
=
0.012 (very weak coupling; just above the bifurcation) there is an unstable focus at F
12 and a stable oval limit cycle of small size representing low-amplitude quasi-harmonic oscillations (Figure 4(c)). On further increasing of coupling strength the limit cycle continuously grows in size and takes an irregular shape indicating nonlinearity of the synchronous oscillations (Figures 4(d)–4(f)).The development of a limit cycle in model (1a)–(1d) in relation to the coupling strength. The projections of phase trajectories and fixed points on (
u
1
,
v
1
) plane are presented. “•” and “∘” mark the respective stable and unstable internal steady state F
12; “×” stands for the boundary steady state F
1. (a) ϰ
1,2
=
0: CR pairs are decoupled; each has unique stable steady state. (b) ϰ
1,2
=
0.011: very weak coupling; the spiral winds only algebraically fast. The system is on the verge of the Hopf bifurcation. (c) ϰ
1,2
=
0.012: very weak coupling; the stable limit cycle is just born. (d) ϰ
1,2
=
0.025: weak coupling; nonlinear oscillations. (e) ϰ
1,2
=
0.05: moderate coupling; the limit cycle passes near the basin of the saddle point. (f) ϰ
1,2
=
0.6: moderate coupling; the limit cycle is about to merge with a heteroclinic cycle.
(a)
(b)
(c)
(d)
(e)
(f)As a practical matter, the range of very weak coupling not too far away from the Hopf bifurcation, where oscillations are quasi-linear and quasi-harmonic, is of less concern to us than is the range of far more feasible not-too-weak coupling, corresponding to well-developed substantially nonlinear oscillations. We are going to demonstrate that given conditions (20), system (1a)–(1d) exhibits relaxation oscillatory behavior, with the two coupled CR pairs being antiphase locked.By the assumption,0
<
ε
≪
1, meaning that system (1a)–(1d) is singularly perturbed. The slow variables are resources, u
1 and u
2, and the fast variables are consumers, v
1 and v
2. The standard practice of reducing such systems is the adiabatic elimination of the fast variables, when the left-hand side in the fast equation is replaced by zero, thus turning this differential equation into an algebraic equation. It is assumed that the fast variables quickly relax to their momentary equilibrium, quasi-steady-state, values obtained from the algebraic equations, in which the slow variables are treated as parameters. “Frozen” slow variables do not move substantially in this short adaptation time of the fast variables. Quasi-steady-state values of the fast variables can then be expressed by values of the slow variables. The fast variables hastily adapt to the motion of the slow variables. The former are entrained by the latter. The utility of quasi-steady-state approximation is that it allows us to reduce the dimension of the system by retaining only slow variables in the model. One has to establish the validity of the adiabatic elimination in each specific case by using the recommendations of the singular perturbation theory [22]. In particular, Tikhonov’s theorem [20] requires the quasi-steady state of the fast equations to be stable.To decompose the full system (1a)–(1d) into fast and slow subsystems, introduce fast time variable τ
=
t
/
ε. Now rescale (1a)–(1d) by replacing t with τ
ε and, after taking ε
=
0, it becomes(42a)
u
1
′
=
u
2
′
=
0
,
(42b)
v
1
′
=
(
u
1
-
δ
v
1
-
ϰ
2
v
2
)
v
1
,
(42c)
v
2
′
=
(
u
2
-
δ
v
2
-
ϰ
1
v
1
)
v
2
,where prime means differentiation with respect to τ. This is the fast subsystem, where u
1 and u
2 are replaced by their initial values and treated as parameters. It yields inner solution, valid for t
=
O
(
ε
).Settingε
=
0 in (1a)–(1d) leads to the slow subsystem(43a)
u
˙
1
=
γ
1
-
(
u
1
+
1
)
v
1
-
u
1
,
(43b)
u
˙
2
=
γ
2
-
(
u
2
+
1
)
v
2
-
u
2
,
(43c)
0
=
(
u
1
-
δ
v
1
-
ϰ
2
v
2
)
v
1
,
(43d)
0
=
(
u
2
-
δ
v
2
-
ϰ
1
v
1
)
v
2
,which produces outer solution, valid for t
=
O
(
1
). In this singular limit as ε
→
0, the subsystem defines a slow flow on the surface (slow manifold) given by (43c) and (43d). Outer solution is valid for those u
1 and u
2, for which the quasi-steady states of the fast subsystem are stable.We anticipate the dynamics of the full system (1a)–(1d) in its four-dimensional phase space (
u
1
,
u
2
,
v
1
,
v
2
) to consist of two typical motions: quickly approaching the slow manifold (43c) and (43d), and slowly sliding over it until a leave point (where the solution disappears) is reached. After that, the representing point may possibly jump to another local solution of (43c) and (43d).Thus, we ought to find all quasi-steady states of the fast subsystem (42a)–(42c), map the domains of their stability onto the slow phase plane (
u
1
,
u
2
), and then investigate the dynamics of the slow subsystem (43a)–(43d) with piecewise continuous functions.Fast subsystem (42a)–(42c), which is nothing but the conventional LVG model, has four quasisteady-states—three boundary and one interior—denoted by Q (the slow variables are deemed to be frozen):(44a)
Q
:
v
~
1
=
0
,
v
~
2
=
0
;
(44b)
Q
1
:
v
~
1
=
u
1
δ
,
v
~
2
=
0
;
(44c)
Q
2
:
v
~
1
=
0
,
v
~
2
=
u
2
δ
;
(44d)
Q
12
:
v
~
1
=
ϰ
2
u
2
-
δ
u
1
ϰ
1
ϰ
2
-
δ
2
,
v
~
2
=
ϰ
1
u
1
-
δ
u
2
ϰ
1
ϰ
2
-
δ
2
.Existence and stability of these quasi-steady states are determined by (
u
1
,
u
2
)—position of a representing point in the phase plane of the slow subsystem, shown in Figure 5(a).(a) Phase plane of the slow subsystem (43a)–(43d) sectored (by dashed lines) into stability domains of the corresponding quasi-steady states of the fast subsystem (42a)–(42c). Both boundary quasi-steady states, each corresponding to the situation when either of the two consumers completely dominates, are stable within the opening of the angle formed by the dashed lines. Lines u
1
=
γ
1
δ and u
2
=
γ
2 are the respective nullclines u
˙
1
=
0 and u
˙
2
=
0 of the piecewise subsystem (47a)-(47b). Lines u
1
=
γ
1 and u
2
=
γ
2
δ mean the same for the piecewise subsystem (52a)-(52b). Intersections of the nullclines (marked by the open circles) are equilibria of the associated piecewise slow subsystems, and they must lie outside the above-mentioned opening to allow for the relaxation oscillations. (b), (c), and (d) are the respective phase portraits of the fast subsystem (42a)–(42c) generated by points “b,” “c,” and “d” in the slow phase plane (a).
(a)
(b)
(c)
(d)Q, Q
1 and Q
2 always exist for all u
1 and u
2 from the positive quadrant of the slow phase plane. For not-too-weak coupling, such that ϰ
1
ϰ
2
>
δ
2, Q
12 exists for all u
1 and u
2 satisfying the condition δ
u
1
/
ϰ
2
<
u
2
<
u
1
ϰ
1
/
δ, that is, within the opening of the angle formed by lines δ
u
1
-
ϰ
2
u
2
=
0 and ϰ
1
u
1
-
δ
u
2
=
0 in Figure 5(a). The opening shrinks as the coupling strengths get weaker.Jacobian matrix of the fast subsystem,(45)
(
u
1
-
2
δ
v
~
1
-
ϰ
2
v
~
2
ε
0
-
ϰ
2
v
~
1
ε
-
ϰ
1
v
~
2
ε
u
2
-
2
δ
v
~
2
-
ϰ
1
v
~
1
ε
)
,
has the following sets of eigenvalues at (44a)–(44d):(46a)
λ
(
Q
)
:
u
1
ε
,
u
2
ε
;
(46b)
λ
(
Q
1
)
:
-
ϰ
1
u
1
ε
δ
+
u
2
ε
,
-
u
1
ε
;
(46c)
λ
(
Q
2
)
:
-
ϰ
2
u
2
ε
δ
+
u
1
ε
,
-
u
2
ε
;
(46d)
λ
(
Q
12
)
:
±
u
1
u
2
+
O
(
δ
)
ε
.Based on (46a)–(46d) one concludes that Q (the origin) is always an unstable node for all u
1 and u
2 from the positive quadrant of the slow phase plane. Q
1 is a stable node for δ
u
2
<
ϰ
1
u
1, that is, below the line ϰ
1
u
1
-
δ
u
2
=
0 in Figure 5(a), otherwise it is a saddle. Similarly, Q
2 is a stable node for δ
u
1
<
ϰ
2
u
2, that is, above the line δ
u
1
-
ϰ
2
u
2
=
0 in the plane of slow variables, otherwise it is a saddle. The interior quasi-steady state, Q
12, is always a saddle.The performed typology of fixed points of the fast subsystem (42a)–(42c) leads to three qualitatively different phase portraits depicted by Figures 5(b)–5(d).Suppose initiallyQ
1 is stable and Q
2 is not. Consumer 1 completely dominates. This corresponds to slow variables u
1 and u
2 being somewhere below the line δ
u
1
-
ϰ
2
u
2
=
0 of Figure 5(a). Fast subsystem (42a)–(42c) has phase portrait of a type shown in Figure 5(b). While u
1 remains much greater than γ
1
δ, the dynamics of the resources (treated as bifurcation parameters in reference to the consumers) is described by a system of two independent equations:(47a)
u
˙
1
=
γ
1
-
(
u
1
+
1
δ
+
1
)
u
1
,
(47b)
u
˙
2
=
γ
2
-
u
2
,which is a piecewise version of the slow subsystem (43a)–(43d) for v
1
=
u
1
/
δ and v
2
=
0. System (47a) and (47b) has stable steady state
(48)
u
^
1
(
1
)
=
1
2
(
r
1
-
δ
-
1
)
=
γ
1
δ
+
O
(
δ
2
)
,
u
^
2
(
1
)
=
γ
2
,
where r
1
=
1
+
δ
(
4
γ
1
+
2
+
δ
). This equilibrium lies in the upper left corner of Figure 5(a).While heading to (48), the trajectory crosses the line δ
u
1
-
ϰ
2
u
2
=
0 and enters the domain of bistability of both Q
1 and Q
2. Fast subsystem (42a)–(42c) takes new phase portrait of a type presented by Figure 5(c). However the dominance of consumer 1 persists.Upon introducing the deviationsξ
1 and η
2 from the respective steady states (48), the system (47a)-(47b) become(49a)
ξ
˙
1
=
-
ξ
1
(
r
1
+
ξ
1
)
δ
,
(49b)
η
˙
2
=
-
η
2
,whence one finds
(50)
ξ
1
(
t
)
=
r
1
(
r
1
/
ξ
1
(
0
)
+
1
)
exp
⁡
⁡
(
r
1
t
/
δ
)
-
1
,
η
2
(
t
)
=
η
2
(
0
)
exp
⁡
(
-
t
)
.It follows from (50), that the dynamics of variable u
1 is faster than that of u
2 due to small δ. Clearly, the representing point must have relaxed to the vertical line u
1
=
γ
1
δ well before approaching the horizontal line u
2
=
γ
2. Further developments depend on whether the involved individual CR systems are overdamped or underdamped.Ifδ is great enough to damp intrinsic oscillations in the constituent CR pairs, the representing point will slide along the nullcline u
˙
1
=
0 (slow manifold of the system (47a)-(47b)), steadily tending to point (48) (Figure 6(a)).The limit cycle of system (1a)–(1d) projected onto phase plane of the slow variables in the cases of (a) strong second-order damping, and (b) underdamping. The parameters of the model chosen for numerical simulation are ε
=
0.727273
×
1
0
-
3, γ
1
=
1.19375, γ
2
=
1, ϰ
1
=
0.5, ϰ
2
=
0.8, δ
=
0.1 (a), and δ
=
0.01 (b).
(a)
(b)Ifδ is small in terms of the condition (6), then the individual CR pairs are underdamped. In the immediate vicinity of u
^
1
(
1
) the division of variables on slow and fast ones loses its meaning and reduced (47a) is no longer valid. Instead of (47a) one should write down two equations:(51a)
u
˙
1
=
γ
1
-
(
u
1
+
1
)
v
1
-
u
1
,
(51b)
ε
v
˙
1
=
(
u
1
-
δ
v
1
)
v
1
.However this system is identical to (3a) and (3b) for an uncoupled CR system and, in essence, describes convergence to a focal point via dying oscillations. In the plane of slow variables (
u
1
,
u
2
) these oscillations manifest themselves in damped transverse fluctuations superimposed on the independent vertical motion along the nullcline u
˙
1
=
0 (u
1
=
u
^
1
(
1
)
≈
γ
1
δ) toward point (48) (Figure 6(b)).By virtue of the condition (20), the trajectory has to cross the line ϰ
1
u
1
-
δ
u
2
=
0 on its way toward the neighborhood of steady state (48). As soon as this has happened, node Q
1 in the plane (
v
1
,
v
2
) will be absorbed by saddle Q
12. A new phase portrait of the fast subsystem (42a)–(42c) takes on the appearance of Figure 5(d). Consumer 1 rapidly washes out, and the alternative boundary quasi-steady state Q
2 becomes stable, with consumer 2 dominating.In terms of the four-dimensional phase space of full system (1a)–(1d), the representing point is now in the other stable branch of the slow manifold (43c) and (43d). The motion over this alternative branch obeys the piecewise subsystem:(52a)
u
˙
1
=
γ
1
-
u
1
,
(52b)
u
˙
2
=
γ
2
-
(
u
2
+
1
δ
+
1
)
u
2
,with the initial conditions u
1
(
0
)
=
u
^
1
(
1
)
≈
γ
1
δ and u
2
(
0
)
=
γ
1
ϰ
1.The dynamics of (52a) and (52b) is basically similar to that of (47a) and (47b) analyzed above. System (52a) and (52b) has a stable steady state:
(53)
u
^
1
(
2
)
=
γ
1
,
u
^
2
(
2
)
=
1
2
(
r
2
-
δ
-
1
)
=
γ
2
δ
+
O
(
δ
2
)
,
where r
2
=
1
+
δ
(
4
γ
2
+
2
+
δ
). This equilibrium lies in the lower right corner of Figure 5(a).Variableu
2, being more rapid in comparison to u
1, quickly enters the neighborhood of the nullcline u
˙
2
=
0 given by u
2
=
u
^
2
(
2
)
≈
γ
2
δ and then—depending on the value of δ—finally approaches the nullcline either monotonically (Figure 6(a)) or via damped oscillations according to equations(54a)
u
˙
2
=
γ
2
-
(
u
2
+
1
)
v
2
-
u
2
,
(54b)
ε
v
˙
2
=
(
u
2
-
δ
v
2
)
v
2
.(Figure 6(b)). System (54a)-(54b) describes underdamped intrinsic oscillations of uncoupled consumer 2 for small δ.At the same time,u
1 steadily and independently tends to u
^
1
(
2
)
=
γ
1. Again, because point (53) is located below the line δ
u
1
-
ϰ
2
u
2
=
0 (on the strengths of the condition (20)), the trajectory would certainly cross that line at a point (
γ
2
ϰ
2
,
γ
2
δ
), whereupon node Q
2 would be absorbed by saddle Q
12. The system returns to the first branch of the slow manifold, and thereby the oscillatory cycle gets closed.


## 4. Results and Discussion

Figure7 shows the results of a numerical integration of system (1a)–(1d) for the case of underdamped individual consumer-research pairs. The two coupled communities execute self-sustained relaxation oscillations which are antiphase-locked.Time profiles of antiphase relaxation oscillations in two coupled CR pairs modeled by (1a)–(1d). (a) and (b) show the, respective, resource 1 and consumer 1; (c) and (d) display the, respective, resource 2 and consumer 2. Numerical values of the parameters are those mentioned in the caption of Figure 6(b).
(a)
(b)
(c)
(d)The resourcesu
1 and u
2 demonstrate sawtooth periodic pulses. The oscillation range for the resource levels remains finite and, what is important, it does not depend on the intrinsic second order loss δ (measuring intraspecific interference).The times of motion over either branch of the slow manifold (43c) and (43d) add up to give a predominant contribution to the period of oscillations, T. These times are determined mainly by the dynamics of the slow resource variables u
1 and u
2 and, to a zeroth approximation in ε and δ, can be found as solutions to the equations of motion (52a) and (47b) with respective boundary conditions (
0
,
γ
2
ϰ
2
) and (
0
,
γ
1
ϰ
1
). In this way one obtains a quite simple estimate for the period:
(55)
T
=
∫
0
γ
2
ϰ
2
‍
d
z
γ
1
-
z
+
∫
0
γ
1
ϰ
1
‍
d
z
γ
2
-
z
=
ln
⁡
1
1
-
ϰ
2
(
γ
2
/
γ
1
)
+
ln
⁡
1
1
-
ϰ
1
(
γ
1
/
γ
2
)
.
It is interesting that, according to (55), the period depends on the ratio of the two resource inflows, γ
1 and γ
2, rather than on each of them individually, and does not depend completely on concrete value of δ.The consumersv
1 and v
2 change periodically between extinction and respective constant levels γ
1 and γ
2. Very brief transient from zero to flat nonzero level within each cycle is accompanied by a highly pronounced spiky overshoot. The magnitude of the spike tends to infinity as δ
→
0, in view of (44b) and (44c). Depending on the intensity of intraspecific interference, the overshoot may or may not be followed by a tail of fading high-frequency oscillations, when a consumer variable falls below its steady-state value and then bounces back above, taking some time to settle close to its steady-state value. In signal processing, such a kind of transient oscillation is known as “ringing.” There is no ringing if the involved CR pair does not oscillate due to significant intraspecific interference. Ringing takes place if intraspecific interference is negligible and therefore the involved CR pair is characterized by underdamped intrinsic oscillations. The “pitch” of ringing is just the frequency of these intrinsic oscillations.One notices that when one consumer is very scanty, the coupled system behaves like an isolated CR pair (2a) and (2b). Another essential feature of the dynamics is the role of the resource variables in determining when the consumers emerge and wash out. When u
1, for example, rises above a threshold value (determined by the amount of losses experienced by v
1), then v
1 comes into dominance causing v
2 in turn to disappear. So except for their transient spiking and ringing, the consumer levels, either flat nonzero or essentially zero, are determined by the hysteretic cycling of the respective resources.Scrutinizing a cycle of consumer oscillations one may distinguish four parts within it:(1)
v
1 is essentially zero, while v
2 is approximately equal to its uncoupled steady-state value, γ
2. u
1 increases due to its inflow until it overcomes losses for consumer 1;
(2)
with a sufficient resource stock,v
1 now emerges. The population exhibits a spike due to the fast time scale of the consumer equations. The sharp increase in population saturates the available resource level, so u
1 drops. Cross-losses cause v
2 to wash out;
(3)
v
1 and u
1 relax to quasi-steady-state values, as if there were only one isolated CR pair. v
2 is essentially zero. u
2 is increasing, like u
1 did in part 1;
(4)
u
2 surpasses the losses, v
2 emerges and the subsequent cross-losses cause v
1 to wash out. The spiking v
2 also causes a substantial decrease in the available stock of the associated resource. The sequence begins again.Presenting his famous model Smale remarked that “it is more difficult to reduce the number of chemicals to two or even three” [6]. As distinct from Smale’s example, coupling in our case makes self-sustained synchronous oscillations possible for just two variables.As we have seen, phase trajectory of the system constantly moves from the neighborhood of unstable boundary equilibriumF
1, where only consumer 1 is present, to the neighborhood of F
2, where consumer 2 completely dominates, back to F
1, and so on in cyclic alternation. This kind of trajectory was termed “heteroclinic cycle” by Kirlinger [23]. A heteroclinic cycle occurs when the outflow (unstable manifold) from one saddle point is directly connected to the inflow (stable manifold) of another saddle point, and vice versa. It is closely related to another notion of the nonlinear dynamics, a homoclinic cycle, which emerges when the unstable and the stable manifolds of the same saddle coincide and form a closed loop.Homo- and heteroclinic cycles are not robust structures in the sense that infinitesimally small change of system parameters destroy them. However in the practical sense, any limit cycle passing in close proximity to saddle points will be indistinguishable from a heteroclinic cycle (Figure8). The only difference is strict periodicity, although the period of the limit cycle in a neighborhood of the heteroclinic cycle may be long. Besides, at the threshold of homo-/heteroclinic bifurcation the period is susceptible to external noise.Figure 8
A 3D-projection of the limit cycle in system (1a)–(1d) for parameters chosen in a neighborhood of the heteroclinic cycle. See caption to Figure 6(b) for the parameters.In the context of our model, as coupling becomes stronger, the stable limit cycle swells and passes closer and closer to boundary fixed points which are node-saddles or focus-saddles (Figures4(d)–4(f)). Depending on the interplay between the parameters, eventually it may bang into one or both of these equilibria creating either a homoclinic or heteroclinic cycle, respectively. This corresponds to γ
2
/
γ
1
=
ϰ
1 and γ
1
/
γ
2
=
ϰ
2. On further increasing the coupling, the saddle connection breaks and the loop is destroyed.It is worth noting that heteroclinic cycles were first found by May and Leonard [24] in a classical LVG system with competing three species. However in their model the cycle is not truly periodic: as time goes on, the system tends to stay in the neighborhood of any one boundary equilibrium ever longer, so that the “total time spent in completing one cycle is likewise proportional to the length of time the system has been running.” Moreover, May and Leonard state that “the phenomenon clearly requires at least three competitors, which is why it cannot occur in models with two competitors.” This statement is echoed by Vandermeer [25] who extended their theory on higher dimensions: “It appears to be the case that all cases of an odd number of species follow this basic pattern, whereas all cases of even number of species result in extinction of half of the components, leaving the other half living independently at their carrying capacities.” In view of our results, the above conclusion is by far and away true providing one stays within the framework of classical LVG equations, which in fact imply a high rapidity of the resource dynamics. In our model of just two competitors the slowness of the resource relative to the consumer is essential for the oscillations to occur, because it provides the necessary inertia to the system.Physically, our model is most likely feasible because it is based on the well-established rate equations (A.8b) for semiconductor lasers, and therefore should be considered as a model of antiphase synchronization of two lasers via their loss-coupling.Ecologically, the feasibility of the model is tightly bound to justification of the adopted time hierarchy in system (1a)–(1d). Time scales are usually inverted in ecosystems, as opposed to lasers, the most common case being rapid consumption of food by species. However it seems reasonable to propose that our model may describe the first level of an ecosystem, at which the consumers are autotrophs and the resources are mineral nutrients. The ability to exploit different substrates leads to a possibility of stable coexistence of different organisms descending from the common ancestor. Divergent evolution is just the formation of new species: due to mutations two populations emerge with the same genetic code but having proteins able to process different substrates. Providing the environmental conditions are quite stable on the evolutionary timescale, the inflows of inorganic substrates from the surroundings may be considered constant and the washout time of a substrate may occur much longer than the life expectancy of a species (recall the definition of ε from (A.6)).


## 5. Conclusion

We proposed a model of two CR pairs linked by interspecific interference competition. When uncoupled, an individual CR pair has a unique stable steady state and does not admit periodic solutions. If intraspecific interference within the species is strong enough, the equilibrium is nonoscillatory (stable node), otherwise the steadying occurs by decaying oscillations (stable focus).When coupled, the model behaves differently at strong and weak competitive interaction between the consumers. When coupling is strong, one of the consumers wins. Which consumer wins or loses depends critically on the relative intensities of the resource inflows and coupling strengths. In the case of bistability, when the system acts like a bistable multivibrator (flip-flop circuit), the winner may be determined by the initial conditions. Any static coexistence of competing consumers is not possible.When coupling is moderately weak, the model reveals low-frequency antiphase relaxation oscillations represented by a continuous flow of rectangular pulses. The system works as an astable multivibrator continually switching between its two states, neither of which is stable. The consumers cannot coexist even dynamically: in each of two alternating states one consumer completely dominates and the other is on the verge of extinction. The most intriguing feature of the model is that each of the participating CR pairs taken separately does not oscillate; both communities are completely quiescent, however, in interaction, when coupled in a nonlinear way, the resulting system turns into a relaxation oscillator.One way or the other, it is believed that the proposed model for coupling-induced oscillations in nonoscillatory CR pairs can be considered as a minimal in that class of population-dynamical systems and its mechanism can be applied to networks with large numbers of nonoscillatory elements and complex architecture.



---

*Source: 561958-2014-09-10.xml*