Awakened oscillations in coupled consumer-resource pairs.
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  constructed a counterintuitive mathematical example of a biological cell modeled by the chemical kinetics of four metabolites, [x.sub.1], ..., [x.sub.4], such that the reaction equations dx/dt = R(x) for the set of metabolites, x = ([x.sub.1], ..., [x.sub.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.sub.2] - [x.sub.1]), where M is a diagonal matrix with the elements [[mu].sub.k][[delta].sub.kl], 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 ; 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 ; an analog operational-amplifier implementation of neural cells connected by passive coupling (where conductance of the resistive connection simulates the diffusion coefficient) , 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 . 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 ), 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) . The LVG model operates with carrying capacities of the species, rather than referring explicitly to any essential resources. As shown by MacArthur , LVG equations maybe 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 . 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 . These equations have been treated mostly numerically. The notable analytical result belongs to Erneux and Mandel  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:
[[??].sub.1] = [[gamma].sub.1] - ([u.sub.1] + 1) [v.sub.1] - [u.sub.1], (1a)
[[??].sub.2] = [[gamma].sub.2] - ([u.sub.2] + 1) [v.sub.2] - [u.sub.2], (1b)
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (1c)
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (1d)
Here dots indicate differentiation with respect to the nondimensional time variable t, [u.sub.i] and [v.sub.i] (i = 1,2) are quantities measuring the respective population sizes of ith resource and ith consumer, [[gamma].sub.i] > 0 (i = 1,2) is the inflow rate of ith resource, 0 < [delta] [much less than] 1 is a parameter representing consumer self-limitation, [[??].sub.j] > 0 (j = 1,2 and j [not equal to] i) quantifies the inhibitory effect of jth consumer on the growth of ith consumer due to coupling, and 0 < [epsilon] [much less than] 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.sub.i] directly represents population density of consumer species and is always nonnegative. Quantity [u.sub.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 [right arrow] aN + 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.sub.i] = -1 corresponds to zero population size in reality. Nevertheless, from here on we shall apply the term "resource" to [u.sub.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 [[??].sub.1,2] = 0, the communities are uncoupled and completely independent. An isolated CR pair is governed by equations:
[??] = [gamma] - (u + 1) v - u, (2a)
[epsilon][??] = (u - [delta]v)v. (2b)
There exist two nonnegative steady states:
[bar.u] = [gamma], [bar.v] = 0; (3a)
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (3b)
The linearization of (2a) and (2b) takes the form
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (4)
where ([bar.u], [bar.v]) is one of the above steady states (3a) and 3b).
At (3a), one eigenvalue is negative and one is positive: [[lambda].sub.1] = -1, [[lambda].sub.2] = [gamma]/[epsilon]. Thus (3a) is a saddle point. At (3b), Tr J = -([delta]/[epsilon] + 1)[bar.v] - 1 < 0 and det J = (1 + 2[delta][bar.v] + [delta])[bar.v]/[epsilon] > 0, so the steady state is a stable node/focus. Specifically, focus is the case for
[(Tr J).sup.2] - 4 det J = ([[gamma].sup.2][[delta].sup.2] + O([[delta].sup.3]))[[epsilon].sup.-2] + (-4[gamma] + O([delta])) [[epsilon].sup.-1] + O(1) < 0, (5)
whence one obtains an asymptotic estimate
[delta] = o([[epsilon].sup.1/2]). (6)
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.
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 [epsilon] and [delta], and is worth another look. There are in fact three timescales, O([epsilon]), O(1), and O([delta]), involved in the CR system (2a) and (2b) when it is overdamped, and two--O([epsilon]) 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 [epsilon]. 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.
Replacing t in (2a) and 2b) with a fast time variable [tau] = t/[epsilon] and setting [epsilon] = 0 we obtain the fast subsystem:
u' = 0, (7a)
v' = (u - [delta]v)v, (7b)
where prime means differentiation with respect to [tau]. In the stretched timescale [tau] 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
v([tau]) = [u/[delta]]/[1 + (u/[delta]v(0) - 1)exp(-u[tau])], (8)
valid for t = O([epsilon]).
After a lapse of considerable time (in the fast scale) v converges to either of two fixed points depending on a sign of u:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (9)
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 - [delta]v =0 practically in a finite time of order O([epsilon]). This is shown in Figures 1(a) and 1(b).
Now set [epsilon] = 0 in (2a) and (2b) to get the slow subsystem:
[??] = [gamma] - (u + 1) v - u, (10a)
0 = (u - [delta]v) v. (10b)
Equation (10b) describes a slow manifold consisting of two lines in the (u, v) plane: v = u/[delta] 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), [bar.v] = u/[delta], is an isolated root of 10b) and [bar.v] is a stable solution of (10b) for any u > 0; the assumptions of Tikhonov's theorem  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/[delta], 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 [delta]. Namely, close to a stable node, the system has two distinct real negative eigenvalues, one fast ([[lambda].sub.1]), and one slow ([[lambda].sub.2]):
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (11)
Since [absolute value of [[lambda].sub.1]] [much greater than] [absolute value of [[lambda].sub.2]], trajectories starting off the associated eigenvector 2 (which is tangent to the slow manifold v = u/[delta]) 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([delta]).
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 manifold v = u/[delta] (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 [delta] [right arrow] 0:
[[lambda].sub.1,2] = -[1/2]([gamma] + 1) [+ or -] i[square root of [gamma]/[epsilon]] for [delta] = 0. (12)
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 [epsilon] 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).
To O(1) for small [delta] there is a way to recast (2a) and 2b) near the focal equilibrium (0, [gamma]) in a convenient form where the perturbation parameter [epsilon] does not multiply any right hand side. Namely, performing the scaling t = [[mu].sup.2]s, u = [[mu].sup.2][gamma][xi], and v = [gamma](1 + [eta]), where [[mu].sup.2] = [square root of [epsilon]/[gamma]], one obtains
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (13)
Here dots stand for differentiation with respect to the time variable s.
Equations (13) represent a weakly perturbed Hamiltonian system. For [mu] = 0 the system is pure Hamiltonian and admits the first integral
H = [1/2][[xi].sup.2] + [eta] - ln (1 + [eta]), (14)
which is a conserved quantity ([??] = 0). The periodic solutions of the Hamiltonian system form a one-parameter family with the equilibrium (0, 0) as center point. The condition 0 < [mu] [much less than] 1 makes (13) quasi-conservative with phase trajectories slowly spiralling to the equilibrium (Figure 2).
System (13) can be rewritten as a second order differential equation for [eta] only:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (15)
This is the equation for a nonlinear quasi-conservative oscillator. As is seen from (15), for [absolute value of [eta]] [much less than] 1 and [absolute value of [??]] [much less than] 1 the oscillator is not only quasi-conservative, but also quasilinear with the frequency [[omega].sub.0] =1. These conditions motivate introducing the new variable [zeta] = [eta]/[mu]. As a result one obtains
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (16)
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, [[??].sub.1,2] > 1, and weak, [[??].sub.1,2] < 1.
Physically feasible equilibria, or fixed points, ([[bar.u].sub.1], [[bar.u].sub.2], [[bar.v].sub.1], [[bar.v].sub.2]), of (1a)-(1d) are those for which [[bar.v].sub.1,2] [greater than or equal to] 0. We denote the interior fixed point by [F.sub.12] = ([[bar.u].sub.1], [[bar.u].sub.2], [[bar.v].sub.1], [[bar.v].sub.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.sub.1] = ([[bar.u].sub.1], [[bar.u].sub.2], [[bar.v].sub.1], 0) and [F.sub.2] = ([[bar.u].sub.1], [[bar.u].sub.2], 0, [[bar.v].sub.2]) designate either of one-consumer equilibria corresponding to dominance, while F = ([[bar.u].sub.1], [[bar.u].sub.2], 0, 0) means both consumers having been washed out.
Model (1a)-(1d) has four feasible steady states. To O(1) for small [delta][right arrow]
F: [[bar.u].sub.1] = [[gamma].sub.1], [[bar.u].sub.2] = [[gamma].sub.2], [[bar.v].sub.1] = 0, [[bar.v].sub.2] = 0; (17a)
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]; (17b)
[F.sub.12]: [[bar.u].sub.1] = [[[??].sub.1][[gamma].sub.1] - [[??].sub.2][[gamma].sub.2] - [[??].sub.1][[??].sub.2] + 1 [+ or -] R/2([[??].sub.1] - 1)], (17c)
[[bar.v].sub.1] = [[[bar.u].sub.2]/[[??].sub.1]], [[bar.v].sub.2] = [[bar.u].sub.1]/[[??].sub.2], (17d)
R = [square root of [([[??].sub.1][[gamma].sub.1] - [[??].sub.2][[gamma].sub.2] - [[??].sub.1][[??].sub.2] + 1).sup.2] + 4[[??].sub.2] ([[??].sub.1] - 1) ([[??].sub.1][[gamma].sub.1] - [[gamma].sub.2])]. (18)
Boundary equilibria F, [F.sub.1], and [F.sub.2] always exist. [F.sub.12] exists if
1/[[??].sub.2] < [[gamma].sub.2]/[[gamma].sub.1] < [[??].sub.1] for [[??].sub.1] > 1 [conjunction] [[??].sub.2] > 1 (strong coupling), (19)
[[??].sub.1] < [[gamma].sub.2]/[[gamma].sub.1] < 1/[[??].sub.2] for [[??].sub.1] < 1 [conjunction] [[??].sub.2] < 1 (weak coupling). (20)
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:
[lambda](F): [[gamma].sub.1]/[epsilon], [[gamma].sub.2]/[epsilon], -1, -1. (21)
Correct to O(1) in [epsilon], the eigenvalues for [F.sub.1] and [F.sub.2] are
[lambda]([F.sub.1]): -1, -[1/2] ([[gamma].sub.1] + 1) [+ or -] i[square root of [[gamma].sub.1]/[epsilon]], [[[gamma].sub.2] - [[??].sub.1][[gamma].sub.1]/[epsilon]]; (22a)
[lambda]([F.sub.2]): 1, -[1/2] ([[gamma].sub.2] + 1) [+ or -] i[square root of [[gamma].sub.2]/[epsilon]], [[[gamma].sub.1] - [[??].sub.2][[gamma].sub.2]/[epsilon]], (22b)
whence it follows that (17b) and 17c) are stable when
[[gamma].sub.2]/[[gamma].sub.1] < [[??].sub.1], (23a)
[[gamma].sub.1]/[[gamma].sub.2] < [[??].sub.2], (23b)
The necessary and sufficient conditions for all the eigenvalues of the Jacobian matrix,
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (24)
evaluated at [F.sub.12], to have negative real parts are, from the Routh-Hurwitz criterion,
[c.sub.0] > 0, (25a)
[c.sub.3] > 0, (25b)
[c.sub.2][c.sub.3] - [c.sub.1] > 0, (25c)
[c.sub.1] ([c.sub.2][c.sub.3] - [c.sub.1]) - [c.sub.0][c.sup.2.sub.3] > 0, (25d)
where [c.sub.0], [c.sub.1], [c.sub.2], and [c.sub.3] are the coefficients of the characteristic polynomial [[lambda].sup.4] + [c.sub.3][[lambda].sup.3] + [c.sub.2][[lambda].sup.2] + [c.sub.1][lambda] + [c.sub.0] of (24):
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (26)
To analyze the validity of (25a)-(25d) we put [[??].sub.2] = [mu][[??].sub.1], where [mu] = O(1). Then for (25a) to be true, [[??].sub.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.sub.12] to have negative real part. For [[??].sub.1] = o(1) it boils down to
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (27)
whence it follows that
[[??].sub.1] = o ([[epsilon].sup.1/2]). (28)
With regard to a fairly small value of [epsilon], 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.sub.12]--if it exists--is always unstable by growing oscillations.
Existence and stability conditions of possible nonnegative equilibrium points are summarized in Table 1.
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.sub.1] and [F.sub.2] are stable with an unstable coexistence steady state [F.sub.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, [[gamma].sub.2] = [[gamma].sup.*.sub.2], while the inflow of resource 1, [[gamma].sub.1], steadily increases from a value less than [[gamma].sup.*.sub.2]/[[??].sub.1] along the line [[gamma].sub.2] = [[gamma].sup.*.sub.2] in the ([[gamma].sub.1], [[gamma].sub.2]) parameter plane, as shown in Figure 3(a). Referring to (17a)-(17d) and 23a) and (23b), one sees that initially [F.sub.1] is unstable; [F.sub.2] is stable and [F.sub.12] does not exist with the complete dominance of consumer 2. Within the interval [[gamma].sup.*.sub.2]/[[??].sub.1] < [[gamma].sub.1] < [[??].sub.2][[gamma].sup.*.sub.2] steady state [F.sub.1] becomes stable yet empty and [F.sub.12] becomes existent (according to (19)) yet unstable, so the situation remains unchanged until point ([[??].sub.2][[gamma].sup.*.sub.2],[[gamma].sup.*.sub.2]) in Figure 3(a) has been reached from the left. For a larger [[gamma].sub.1], state [F.sub.2] gives up its stability and the system jumps to [F.sub.1]. Consumer 2 gets washed out, while consumer 1 takes over. If now we start reducing [[gamma].sub.1], the system remains in [F.sub.1] until [[gamma].sub.1] drops to the lower critical value [[gamma].sup.*.sub.2]/[[??].sub.1], beyond which [F.sub.1] is no longer stable and there is a reverse jump to [F.sub.2]. In other words, as [[gamma].sub.1] progresses along [[gamma].sub.2] = [[gamma].sup.*.sub.2] there is a discontinuous switch from consumer 2 to consumer 1 at [[??].sub.2][[gamma].sup.*.sub.2], while as [[gamma].sub.1] retraces its steps, there is a discontinuous switch from consumer 1 to consumer 2 at [[gamma].sup.*.sub.2]/[[??].sub.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 [[gamma].sub.1] [member of] ([[gamma].sup.*.sub.2]/[[??].sub.1], [[??].sub.2][[gamma].sup.*.sub.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.
3.4. Weak Coupling: Antiphase Relaxation Oscillations. As seen from Table 1, the very existence of interior equilibrium [F.sub.12] in a case of weak coupling (condition (20)) implies instability of both boundary fixed points, [F.sub.1] and [F.sub.2]. System (1a)-(1d) happens to possess four nonnegative steady states, none of them being stable. As we have found, [F.sub.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 . For negligible [delta] the problem reduces to weakly coupled quasi-conservative oscillators. Introduce the new time s = t/[[mu].sup.2], where [[mu].sup.2] = [square root of [epsilon]/[[gamma].sub.1]]. The new dynamic variables [[xi].sub.1,2] and [[eta].sub.1,2] are defined by formulas [[xi].sub.1,2] = [u.sub.1,2]/[[mu].sup.2][[gamma].sub.1,2] and [[eta].sub.1,2] = ([v.sub.1,2]-[[gamma].sub.1,2])/[[gamma].sub.1,2]. Thus the variables [[xi].sub.1,2] and [[eta].sub.1,2] are the respective deviations of [u.sub.1,2] and [v.sub.1,2] from their standalone steady-state values. With these new variables equations (1a)-(1d) become
[[??].sub.1] = - [eta]1 - [[mu].sup.2][[xi].sub.1](1 + [[gamma].sub.1] (1 + [[eta].sub.1])), (29a)
[[??].sub.1] = [[xi].sub.1] (1 + [[eta].sub.1]) - ([[gamma].sub.2]/[[gamma].sub.1]) [[rho].sub.2] (1 + [[eta].sub.1]) (1 + [[eta].sub.2]), (29b)
[[??].sub.2] = -[[eta].sub.2] - [[mu].sup.2][[xi].sub.2] (1 + [[gamma].sub.2] (1 + [[eta].sub.2])), (29c)
[[??].sub.2] = ([[gamma].sub.2]/[[gamma].sub.1])[[xi].sub.2](1 + [[eta].sub.2]) - [[rho].sub.1] (1 + [[eta].sub.1])(1 + [[eta].sub.2]), (29d)
where dots mean differentiation with respect to s and [[rho].sub.1,2] = [[??].sub.1,2]/[[mu].sup.2] are the scaled coupling strengths.
For brevity, we restrict our consideration to the case of identical CR pairs and symmetric coupling setting [[gamma].sub.1,2] = [gamma] and [[rho].sub.1,2] = [rho]. 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 [[xi].sub.1] and [[xi].sub.2]. Secondly, (29b) and (29d) are differentiated with respect to time to get [[??].sub.1] and [[??].sub.2]. Thirdly, [[xi].sub.1], [[xi].sub.2], [[??].sub.1], and [[??].sub.2] are plugged into the equations for [[??].sub.1] and [[??].sub.2]. And finally, [[??].sub.1,2] are replaced by the new variables [[zeta].sub.1,2]: [[eta].sub.1,2] = [mu][[zeta].sub.1,2].
As a consequence, we arrive at the following coupled equations:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (30)
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (31)
We seek the solution of (30) in the quasi-harmonic form
[[zeta].sub.1,2] = [a.sub.1,2]cos([omega]s + [[phi].sub.1,2]),
[[??].sub.1,2] = -[a.sub.1,2][omega] sin([omega]s + [[phi].sub.1,2]), (32)
where [a.sub.1,2] and [[phi].sub.1,2] are slowly varying amplitudes and phases, and to is the frequency of the synchronous oscillations. Differentiating the assumed form of [[zeta].sub.1,2] and equating the result to the assumed form of [[??].sub.1,2] yields the first pair of relationships between [a.sub.1,2] and [[phi].sub.1,2]:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (33)
Then, differentiating [[??].sub.1,2] and substituting the resulting expression for [[??].sub.1,2] as well as the assumed forms for [[zeta].sub.1,2] and [[??].sub.1,2] into (30) yields the second pair of equations relating [a.sub.1,2] and [[phi].sub.1,2]. Separating into equations for the rate of change of [a.sub.1,2] and [[phi].sub.1,2] one obtains
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (34)
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (35)
To this point no approximations have been made except of the expansion in powers of [mu] in (30). Averaging equations (34) over the period 2[pi]/to and considering [a.sub.1,2], [[phi].sub.1,2], [[??].sub.1,2] and [[??].sub.1,2] to be constants while performing the averaging, one obtains the following equations describing the slow variations of [a.sub.1,2] and [[phi].sub.1,2]:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (36a)
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (36b)
where [phi] = [[phi].sub.1,2] - [[phi].sub.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 ([phi] = const), and the oscillations are periodic with constant amplitudes [a.sub.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, [phi], we see that stable is only antiphase steady-state regime with [phi] = [pi].
Let R 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,
[??] = R - <R>, (37)
and integrate [??] up to an arbitrary function of the amplitudes chosen for simplicity to be equal to zero:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (38)
Now, upon calculation of the Jacobian matrix 9R/ [partial derivative][([a.sub.1], [a.sub.2], [[phi].sub.1], [[phi].sub.2]).sup.T], we can write down the second approximation symbolically as
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (39)
(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 set [a.sub.1] = [a.sub.2] = a enabling us to present the resulting equations of the second approximation as compact as possible:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (40a)
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (40b)
The stable steady-state solution for the phase difference is [phi] = [pi]. Inserting this value into (40a) and putting [??] = 0 we get three different steady-state solutions for a (including the trivial one), of which only one,
a = 2[square root of 3][square root of [rho][[omega].sup.2] - [[mu].sup.2]([gamma] + 1)([[rho].sup.2] + [[omega].sup.2])/[mu][square root of [rho](4[[omega].sup.2] + 7)]], (41)
being stable. Formula (41) indicates a soft transition to sustained oscillations as the coupling parameter p is progressively increased from a critical value. Thus we expect that the Hopf bifurcation is supercritical.
Figure 4 plots ([u.sub.1], [v.sub.1]) phase planes for coupling strength below and above the Hopf bifurcation. When [[??].sub.1,2] = 0 (consumers are decoupled) the internal equilibrium is a stable focus (Figure 4(a)). For [[??].sub.1,2] = 0.011 (very weak coupling; just below the bifurcation) [F.sub.12] is still a stable focus, though a very gently winding one: the decay is slow (Figure 4(b)). For [[??].sub.1,2] = 0.012 (very weak coupling; just above the bifurcation) there is an unstable focus at [F.sub.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)).
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 < [epsilon] [much less than] 1, meaning that system (1a)-(1d) is singularly perturbed. The slow variables are resources, [u.sub.1] and [u.sub.2], and the fast variables are consumers, [v.sub.1] and [v.sub.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 . In particular, Tikhonov's theorem  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 [tau] = t/[epsilon]. Now rescale (1a)-(1d) by replacing t with [tau][epsilon] and, after taking [epsilon] = 0, it becomes
[u'.sub.1] = [u'.sub.2] = 0, (42a)
[v'.sub.1] = ([u.sub.1] - [delta][v.sub.1] - [[??].sub.2][v.sub.2])[v.sub.1], (42b)
[v'.sub.2] = ([u.sub.2] - [delta][v.sub.2] - [[??].sub.1][v.sub.1])[v.sub.2], (42c)
where prime means differentiation with respect to [tau]. This is the fast subsystem, where [u.sub.1] and [u.sub.2] are replaced by their initial values and treated as parameters. It yields inner solution, valid for t = 0([epsilon]).
Setting [epsilon] = 0 in (1a)-(1d) leads to the slow subsystem
[[??].sub.1] = [[gamma].sub.1] - ([u.sub.1] + 1)[v.sub.1] - [u.sub.1], (43a)
[[??].sub.2] = [[gamma].sub.2] - ([u.sub.2] + 1)[v.sub.2] - [u.sub.2], (43b)
0 = ([u.sub.1] - [delta][v.sub.1] - [[??].sub.2][v.sub.2])[v.sub.1], (43c)
0 = ([u.sub.2] - [delta][v.sub.2] - [[??].sub.1][v.sub.1])[v.sub.2], (43d)
which produces outer solution, valid for t = O(1). In this singular limit as [epsilon] [right arrow] 0, the subsystem defines a slow flow on the surface (slow manifold) given by (43c) and (43d). Outer solution is valid for those [u.sub.1] and [u.sub.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.sub.1], [u.sub.2], [v.sub.1], [v.sub.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.sub.1], [u.sub.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):
Q: [[??].sub.1] = 0, [[??].sub.2] = 0; (44a)
[Q.sub.1]: [[??].sub.1] = [u.sub.1]/[delta], [[??].sub.2] = 0; (44b)
[Q.sub.2]: [[??].sub.1] = 0, [[??].sub.2] = [u.sub.2]/[delta]; (44c)
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (44d)
Existence and stability of these quasi-steady states are determined by ([u.sub.1], [u.sub.2])--position of a representing point in the phase plane of the slow subsystem, shown in Figure 5(a).
Q, [Q.sub.1] and [Q.sub.2] always exist for all [u.sub.1] and [u.sub.2] from the positive quadrant of the slow phase plane. For not-too-weak coupling, such that [[??].sub.1][[??].sub.2] > [[delta].sup.2], [Q.sub.12] exists for all [u.sub.1] and [u.sub.2] satisfying the condition [delta][u.sub.1]/[[??].sub.2] < [u.sub.2] < [u.sub.1][[??].sub.1]/[delta], that is, within the opening of the angle formed by lines [delta][u.sub.1] - [[??].sub.2][u.sub.2] = 0 and x1u1 - Su2 = 0 in Figure 5(a). The opening shrinks as the coupling strengths get weaker.
Jacobian matrix of the fast subsystem,
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (45)
has the following sets of eigenvalues at (44a)-(44d):
[lambda](Q): [u.sub.1]/[epsilon], [u.sub.2]/[epsilon]; (46a)
[lambda]([Q.sub.1]): -[[??].sub.1][u.sub.1]/[epsilon][delta] + [u.sub.2]/[epsilon], -[u.sub.1]/[epsilon]; (46b)
[lambda]([Q.sub.2]): -[[??].sub.2][u.sub.2]/[epsilon][delta] + [u.sub.1]/[epsilon], -[u.sub.2]/[epsilon]; (46c)
[lambda]([Q.sub.12]): [+ or -] [square root of [u.sub.1][u.sub.2]] + O([delta])/[epsilon]. (46d)
Based on (46a)-(46d) one concludes that Q (the origin) is always an unstable node for all [u.sub.1] and [u.sub.2] from the positive quadrant of the slow phase plane. [Q.sub.1] is a stable node for [delta][u.sub.2] < [[??].sub.1][u.sub.1], that is, below the line [[??].sub.1][u.sub.1] - [delta][u.sub.2] = 0 in Figure 5(a), otherwise it is a saddle. Similarly, [Q.sup.2] is a stable node for [delta][u.sub.1] < [[??].sub.2][u.sub.2], that is, above the line [delta][u.sub.1] - [[??].sub.2][u.sub.2] = 0 in the plane of slow variables, otherwise it is a saddle. The interior quasi-steady state, Q12, 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 initially [Q.sub.1] is stable and [Q.sub.2] is not. Consumer 1 completely dominates. This corresponds to slow variables [u.sub.1] and [u.sub.2] being somewhere below the line [delta][u.sub.1] - [[??].sub.2][u.sub.2] = 0 of Figure 5(a). Fast subsystem 42a)-(42c) has phase portrait of a type shown in Figure 5(b). While [u.sub.1] remains much greater than [[gamma].sub.1][delta], the dynamics of the resources (treated as bifurcation parameters in reference to the consumers) is described by a system of two independent equations:
[[??].sub.1] = [[gamma].sub.1] - ([u.sub.1] + 1/[delta] = 1)[u.sub.1], (47a)
[[??].sub.2] = [[gamma].sub.2] - [u.sub.2], (47b)
which is a piecewise version of the slow subsystem (43a)-(43d) for [v.sub.1] = [u.sub.1]/[delta] and [v.sub.2] = 0. System (47a) and (47b) has stable steady state
[[??].sup.(1).sub.1] = 1/2([r.sub.1] - [delta] - 1) = [[gamma].sub.1][delta] + O([[delta].sup.2]),
[[??].sup.(1).sub.2] = [[gamma].sub.2], (48)
where [r.sub.1] = [square root of 1 + [delta](4[[gamma].sub.1] + 2 + [delta])]. This equilibrium lies in the upper left corner of Figure 5(a).
While heading to (48), the trajectory crosses the line [delta][u.sub.1] - [[??].sub.2][u.sub.2] = 0 and enters the domain of bistability of both [Q.sub.1] and [Q.sub.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 [[xi].sub.1] and [[xi].sub.2] from the respective steady states (48), the system (47a)-(47b) become
[[??].sub.1] = -[[xi].sub.1]([r.sub.1] + [[xi].sub.1])/[delta], (49a)
[[??].sub.2]= -[[eta].sub.2], (49b)
whence one finds
[[xi].sub.1] = [r.sub.1]/([r.sub.1]/[[xi].sub.1](0) + 1)exp([r.sub.1]t/[delta]) - 1,
[[eta].sub.2](t) = [[eta].sub.2](0)exp(-t). (50)
It follows from (50), that the dynamics of variable [u.sub.1] is faster than that of [u.sub.2] due to small [delta]. Clearly, the representing point must have relaxed to the vertical line [u.sub.1] = [[gamma].sub.1][delta] well before approaching the horizontal line [u.sub.2] = [[gamma].sub.2]. Further developments depend on whether the involved individual CR systems are overdamped or underdamped.
If [delta] is great enough to damp intrinsic oscillations in the constituent CR pairs, the representing point will slide along the nullcline [[??].sub.1] = 0 (slow manifold of the system (47a)-(47b)), steadily tending to point (48) Figure 6(a)).
If S is small in terms of the condition (6), then the individual CR pairs are underdamped. In the immediate vicinity of [[??].sup.(1).sub.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:
[[??].sub.1] = [[gamma].sub.1] - ([u.sub.1] + 1)[v.sub.1] - [u.sub.1], (51a)
[epsilon][[??].sub.1] = ([u.sub.1] - [delta][v.sub.1])[v.sub.1]. (51b)
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.sub.1], [u.sub.2]) these oscillations manifest themselves in damped transverse fluctuations superimposed on the independent vertical motion along the nullcline [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] toward point (18) Figure 6(b)).
By virtue of the condition (20), the trajectory has to cross the line [[??].sub.1][u.sub.2] - [delta][u.sub.2] = 0 on its way toward the neighborhood of steady state (48). As soon as this has happened, node [Q.sub.1] in the plane ([v.sub.1], [v.sub.2]) will be absorbed by saddle Q12. 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.sub.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:
[[??].sub.1] = [[gamma].sub.1] - [u.sub.1], (52a)
[[??].sub.2] = [[gamma].sub.2] - ([u.sub.2] + 1/[delta] + 1)[u.sub.2], (52b)
with the initial conditions [u.sub.1](0) = [[??].sup.(1).sub.1] = [[gamma].sub.1][delta] and [u.sub.2](0) = [[gamma].sub.1][[??].sub.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:
[[??].sup.(2).sub.1] = [[gamma].sub.1], [[??].sup.(2).sub.2] = 1/2([r.sub.2] - [delta] - 1) = [[gamma].sub.2][delta] + O([[delta].sup.2]), (53)
where [r.sub.2] = [square root of 1 + [delta](4[[gamma].sub.2] + 2 + [delta])]. This equilibrium lies in the lower right corner of Figure 5(a).
Variable [u.sub.2], being more rapid in comparison to [u.sub.1], quickly enters the neighborhood of the nullcline [[??].sub.2] = 0 given by [u.sub.2] = [[??].sup.(2).sub.2] [approximately equal to] [[gamma].sup.2][delta] and then--depending on the value of [delta]--finally approaches the nullcline either monotonically (Figure 6(a)) or via damped oscillations according to equations
[[??].sub.2] = [[gamma].sub.2] - ([u.sub.2] + 1)[v.sub.2] - [u.sub.2], (54a)
[epsilon][[??].sub.2] = ([u.sub.2] - [delta][v.sub.2])[v.sub.2]. (54b)
(Figure 6(b)). System (54a)-(54b) describes underdamped intrinsic oscillations of uncoupled consumer 2 for small [delta].
At the same time, [u.sub.1] steadily and independently tends to [[??].sup.(2).sub.1] = [[gamma].sub.1]. Again, because point (53) is located below the line [delta][u.sub.1] - [[??].sub.2][u.sub.2] = 0 (on the strengths of the condition (20)), the trajectory would certainly cross that line at a point ([[gamma].sub.2][[??].sub.2], [[gamma].sub.2][delta]), whereupon node [Q.sub.2] would be absorbed by saddle [Q.sub.12]. The system returns to the first branch of the slow manifold, and thereby the oscillatory cycle gets closed.
4. Results and Discussion
Figure 7 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.
The resources [u.sub.1] and [u.sub.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 [delta] (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.sub.1] and [u.sub.2] and, to a zeroth approximation in e and [delta], can be found as solutions to the equations of motion (52a) and (47b) with respective boundary conditions (0, [[gamma].sub.2][[??].sub.2]) and (0, [[gamma].sub.1][[??].sub.1]). In this way one obtains a quite simple estimate for the period:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (55)
It is interesting that, according to (55), the period depends on the ratio of the two resource inflows, [[gamma].sub.1] and [[gamma].sub.2], rather than on each of them individually, and does not depend completely on concrete value of [delta].
The consumers [v.sub.1] and [v.sub.2] change periodically between extinction and respective constant levels [[gamma].sub.1] and [[gamma].sub.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 [delta] [right arrow] 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.sub.1], for example, rises above a threshold value (determined by the amount of losses experienced by [v.sub.1]), then [v.sub.1] comes into dominance causing [v.sub.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.sub.1] is essentially zero, while [v.sub.2] is approximately equal to its uncoupled steady-state value, [y.sub.2]. [u.sub.1] increases due to its inflow until it overcomes losses for consumer 1;
(2) with a sufficient resource stock, [v.sub.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.sub.1] drops. Cross-losses cause [v.sub.2] to wash out;
(3) [v.sub.1] and [u.sub.1] relax to quasi-steady-state values, as if there were only one isolated CR pair. [v.sub.2] is essentially zero. [u.sub.2] is increasing, like [u.sub.1] did in part 1;
(4) [u.sub.2] surpasses the losses, [v.sub.2] emerges and the subsequent cross-losses cause [v.sub.1] to wash out. The spiking [v.sub.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" . 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 equilibrium [F.sub.1], where only consumer 1 is present, to the neighborhood of [F.sub.2], where consumer 2 completely dominates, back to [F.sub.1], and so on in cyclic alternation. This kind of trajectory was termed "heteroclinic cycle" by Kirlinger . 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 (Figure 8). The only difference is strict periodicity, although the period of the limit cycle in a neighborhood of the heteroclinic cycle maybe long. Besides, at the threshold of homo-/heteroclinic bifurcation the period is susceptible to external noise.
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 (Figures 4(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 [[gamma].sub.2]/[[gamma].sub.1] = [[??].sub.1] and [[gamma].sub.1]/[[gamma].sub.2] = [[??].sub.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  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  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 e from (A.6)).
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 maybe 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 a stable 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.
A. Derivation of the Coupled CR Equations
A.1. Ecological Perspective. Of all types of interactions between individuals of the same population (intraspecific interactions) or individuals of different populations (interspecific interactions) of the same trophic level competition is most commonly encountered. In a broad sense, competition takes place when each species (individual) has an inhibiting effect on the growth of the other species (individual). An inhibiting effect should be understood to mean either an increase in the death rate or a decrease in the birth rate.
Consider the famous CR equations proposed by MacArthur [17, 26]:
[[??].sub.j]([r.sub.j](1 - [x.sub.j]/[K.sub.j]) - [n.summation over (i=1)][c.sub.ij][y.sub.i])[x.sub.j], j = 1, ..., m, (A.1a)
[[??].sub.i] = ([m.summation over (j=1)][c.sub.ij][w.sub.j][x.sub.j] - [b.sub.i])[y.sub.i], i = 1, ..., n. (A.1b)
Here, dots indicate differentiation with respect to time t, [x.sub.j] represents the total biomass of jth resource (prey), [y.sub.i] stands for the total biomass of ith consumer (predator) species, the constant [r.sub.j] defines the growth rate of jth resource, [K.sub.j] is the carrying capacity of jth resource, is the rate of uptake of a unit of jth resource by each individual of ith consumer population, [w.sup.-1.sub.j] is the conversion efficiency parameter representing an amount of jth resource an individual of ith consumer population must consume in order to produce a single new individual of that species, and [b.sub.i] is the loss rate of ith consumer due to either natural death or emigration. All parameters in (A.1a) and (A.1b) are nonnegative.
MacArthur assumed population dynamics of the resources to be much faster than that of the consumers which enabled him to approximate [x.sub.j] in (A.1b) by its quasi-steady-state value derived by setting the right-hand side of (A.1a) to zero. As a result, he succeeded to reduce slow-scale equation (A.1b) to the well-known LVG model :
[[??].sub.i] = ([k.sub.i] - [n.summation over (s=1)][a.sub.is][y.sub.s])[y.sub.i], i = 1, ..., n, (A.2)
[a.sub.is] = [m.summation over (j=1)][c.sub.ij][c.sub.sj]([w.sub.j][K.sub.j][c.sub.ij][c.sub.sj]([w.sub.j][K.sub.j]/[r.sub.j]), i = 1, ..., n; s = 1, ..., n, (A.3a)
[k.sub.i] = [m.summation over (j=1)][c.sub.ij][w.sub.j][K.sub.j] - [b.sub.i], i = 1, ..., n. (A.3b)
More recently, such an asymptotic reduction has also been carried out for a model of competition where species (with continuous trait) consume the common resource that is constantly supplied, under the assumption of very fast dynamics for the supply of the resource and fast dynamics for death and uptake rates .
CR model (A.1a)-(A.1b) assume that competition within and between consumer species is purely exploitative: individuals and populations interact through utilizing (or occupying) common resource that is in short supply. Quite on the contrary, LVG model (A.2), (A.3a)-(A.3b) describes competition strictly phenomenologically, as direct interference where consumers experience harm attributed to their mutual presence in a habitat (e.g., through aggressive behavior). However we have to stress that MacArthur's reduction neither claims that interference competition entirely results from "more fundamental" trophic competition, nor urges us to hastily consider direct competition as some derived concept. What it actually states is that at slow-time scale associated with dynamics of the consumers, the effects of exploitation competition are indistinguishable from those of interference competition. And at slow-time scale, coefficients [a.sub.is] of (A.3a) merely add to interference coefficients [a'.sub.is] which are to be present primordially in (A.1b).
Most mathematical models dealing with coupled CR pairs or multilevel trophic chains ignore contributions of intra-specific and interspecific interference. Indeed, the empirical data like  do indicate that [a'.sub.ij] may be negligible in comparison with [a.sub.is]. Still, works advocating the explicit account for direct interference show that incorporation of self-limitation and cross-limitation terms in the equations at the consumers' level can provide for the stable coexistence of many species on few resources [29, p. 31], .
Moreover, if we are to assume dynamics of the resources to be much slower than that of the consumers, it is likely that we have to retain interference competition terms in all equations (A.1b).
Consider the following modification of (A.1a) and (A.1b) representing coupled two-consumer, two-resource equations:
[[??].sub.1] = [p.sub.1] - ([c.sub.1][y.sub.1] + [q.sub.1])[x.sub.1], (A.4a)
[[??].sub.2] = [p.sub.1] - ([c.sub.2][y.sub.2] + [q.sub.2])[x.sub.2], (A.4b)
[[??].sub.1] = ([c.sub.1][w.sub.1][x.sub.1] - [b.sub.1] - [d.sub.1][y.sub.1] - [h.sub.2][y.sub.2]) [y.sub.1], (A.4c)
[[??].sub.i] = ([c.sub.2][w.sub.2][x.sub.2] - [b.sub.2] - [d.sub.2][y.sub.2] - [h.sub.1][y.sub.1])[y.sub.2]. (A.4d)
Instead of the logistic mode of resource supply, as is the case in MacArthur's model, our model is based on so-called "equable" mode of resource exploitation , by which the quantities of available resources are held constant by a continuous-flow system. According to (A.4a) and (A.4b), a constant concentration of jth resource (j = 1,2) flows into a defined volume with the rate [p.sub.j] while unused resource flows out with the per capita rate [q.sub.j], in much the same manner as in a chemostat .
In more exact terms, the true chemostat model for one substrate and one species looks as follows:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (A.5)
where the rate of substrate uptake is expressed by the Monod formula [[mu].sub.xy]/([K.sub.x] + x), [K.sub.x] is a saturation constant numerically equal to the substrate concentration at which the uptake rate is half the maximum, D is the dilution rate defined as the rate of flow of medium over the volume of the bioreactor, and [x.sub.0] is an input concentration of the substrate.
Model (A.5) turns into an uncoupled version of (A.4a)-(A.4d) if we put p = D[x.sub.0] and q = b = D and assume [K.sub.x] [much greater than] x, so that [mu]x/([K.sub.x] + x) [approximately equal to] cx, where c = [mu]/[K.sub.x].
In natural conditions, the equable modes of feeding, for instance, can be found on the first trophic level of ecosystem, among autotrophs.
Besides, in (A.4c) and (A.4d) intraspecific competition strength [d.sub.i] (i = 1,2) measures direct interference of individuals within ith consumer population with each other resulting in an additional per capita loss rate [d.sub.i][y.sub.i]; interspecific competition strength [h.sub.s] (s = 1,2; s [not equal to] i) quantifies direct interference effect from sth consumer on ith consumer resulting in an additional per capita loss rate, [h.sub.s][y.sub.s], of the latter.
Equations (A.4a)-(A.4d) contain two important assumptions. First, they assume that the resources are noninteractive. On higher trophic levels, however, resources may interact and the possibility of competition among the resources was originally pointed out by Lynch . Since then, a whole series of theoretical papers has been published on two-predator, two-prey systems with interference competition between two self-reproducing prey species based on the Lotka-Volterra equations. Specifically, Kirlinger  describes the model in which each predator specializes on one prey only, while Xiang and Song  treat a similar model in which each predator is allowed to feed on both prey.
As seen from (A.4a) and (A.4b), there is no intraspecific interference competition within the resource populations either. Yet the resource level would remain finite even in the absence of the consumer.
A second assumption of our equations is that the consumers interact only directly, through interference competition and cannot compete through their use of resources, as each consumer specializes in one resource only. The theory of pure trophic competition in equable models has been developed in the works [31, 35].
Intraspecific interference competition is allowed within the consumers as well. Even though the available amount of any resource happened to be of a constant level, the population size of the associated consumer would remain finite due to self-limitation caused by direct intraspecific interference.
The novelty of model (A.4a)-(A.4d) is that it considers time hierarchy of MacArthur's CR equations to be reversed by assuming dynamics of the consumers to be much faster than that of the involved resources and articulates the importance of direct competition mechanisms within the framework of this assumption.
Upon the scaling
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (A-6)
equations (A.4a)-(A.4d) take the following nondimensional form:
[[??].sub.1] = [[gamma].sub.1] - [u.sub.1][v.sub.1] - [u.sub.1] - [v.sub.1], (A.7a)
[beta][[??].sub.2] = [[gamma].sub.2] - [u.sub.2][v.sub.2] - [u.sub.2] - [v.sub.2], (A.7b)
[[epsilon].sub.1][[??].sub.1] = [(.sub.u1] - [[delta].sub.1][v.sub.1] - [x.sub.2][v.sub.2])[v.sub.1], (A.7c)
[[epsilon].sub.2][[??].sub.2] = ([u.sub.2] - [[delta].sub.2][v.sub.2] - [x.sub.1][v.sub.1])[v.sub.2]. (A.7d)
Note that in (A.7a)-(A.7d) dots mean differentiation with respect to nondimensional "slow" timescale variable t', as defined by (A.6).
The parameters [[beta].sup.-1], [[epsilon].sup.-1.sub.1], and [[epsilon].sup.-1.sub.2] reflect the rapidity of the dynamics of [u.sub.2], [v.sub.1], and [v.sub.2] with reference to that of [u.sub.1]. It is assumed that [beta] = O(l) and [[epsilon].sub.1,2] [much less than] 1.
In studying the effect of coupling, the parameters of interest are obviously the coupling strengths, [x.sub.1] and [x.sub.2]. The parameters of interest are also those which characterize the difference between the states of the uncoupled systems. The resource income rates [[gamma].sub.1] and [[gamma].sub.2] are used as the control parameters that distinguish the relative base states of the two systems.
For the sake of simplicity but without any loss of generality, we set [beta] = 1, [[epsilon].sub.1] = [[epsilon].sub.2] = [epsilon], and [[delta].sub.1] = [[delta].sub.2] = [delta] and also drop the prime at Qtoobtain 1a)-(1d).
A.2. Laser Dynamics Perspective. Laser rate equations originally proposed by Statz and deMars  are differential equations that relate two quantities: injected carrier density (n) and photon density (p). For a single-mode semiconductor laser, these equations take the form [12, ch. 6]:
[??] = J/qd = [GAMMA]va(n - [n.sub.0])p - n/[[tau].sub.e], (A.8a)
[??] = [GAMMA][v.sub.g]a(n - [n.sub.0])p - p/[[tau].sub.p], (A.8b)
where dots mean differentiation with respect to time t, J is the injection current density (pump parameter), q is the magnitude of the electron charge, d is the active-layer thickness, r is the confinement factor accounting for the fraction of the light power contained in the active region, [v.sub.g] is the group velocity of light that can be expressed through the speed of light in vacuum (c) and the group refractive index of the dispersive semiconductor material ([[mu].sub.g]) as vg = c/[[mu].sub.g], a is the gain coefficient, [n.sub.0] is the carrier density at transparency corresponding to the onset of population inversion, [[tau].sub.e] is the lifetime of the electrons in the conduction band before being lost by escape from the active region, and [[tau].sub.p] is the lifetime of photons inside the cavity before going out of the cavity or being absorbed inside the cavity. In (A.8b) the contribution of spontaneous emission is neglected.
Typical parameter values for a semiconductor laser (mostly borrowed from [12, p. 238]) are given in Table 2. These numerical values are used in the calculations throughout the present paper, unless otherwise noted.
Consider two (not necessarily identical) lasers of type (A.8a)-(A.8b) and introduce additional intensity-dependent losses such that each laser, i, of the two has a total loss represented by the sum of the constant loss, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], plus the loss proportional to its own intensity, [D.sub.i][p.sub.i], plus the loss proportional to the intensity of the other laser, [h.sub.j][p.sub.j];
[[??].sub.1] = [J.sub.1]/q[d.sub.1] [[GAMMA].sub.1]v[g.sub.1][a.sub.1]([n.sub.1] - [n.sub.01]) [p.sub.1] - [n.sub.1]/[[tau].sub.e1], (A.9a)
[[??].sub.2] = [J.sub.2]/q[d.sub.2] - [[GAMMA].sub.2][v.sub.g2][a.sub.2]([n.sub.2] - [n.sub.02]) [p.sub.2] - [n.sub.2]/[[tau].sub.e2], (A.9b)
[[??].sub.2] = [J.sub.2]/q[d.sub.2] - [[GAMMA].sub.2][v.sub.g][a.sub.2]([n.sub.2] - [n.sub.02]) [p.sub.2] - [n.sub.2]/[[tau].sub.e2], (A.9c)
[[??].sub.1] = ([[GAMMA].sub.1][v.sub.g1][a.sub.1][a.sub.1]([n.sub.1] - [n.sub.01]) - 1/ [[tau].sub.p1] - [D.sub.1][p.sub.1] - [h.sub.2][p.sub.2])[p.sub.1], (A.9d)
[[??].sub.2] = ([[GAMMA].sub.2][v.sub.g2][a.sub.2]([n.sub.2] - [n.sub.02]) - 1/ [[tau].sub.p2] - [h.sub.1][p.sub.1])[p.sub.2], (A.9d)
where [D.sub.i](i = 1,2) is the second order loss constant of the isolated ith laser and ht is the coupling strength measuring the cross-loss effect of ith laser on jth laser.
Thus according to (A.9a)-(A.9d), two lasers happen to be cross-coupled through their resonators, so that each of them can modulate the cavity loss of the other. Technically, intensity-dependent intrinsic and cross-losses may be implemented, for example, using an intracavity electro-optic modulator fed by a current proportional to the output power [37, 38].
We will consider the pump currents [J.sub.1] and [J.sub.2], and the coupling strengths [h.sub.1] and [h.sub.2], as free parameters of the model. For the present, we cannot judge with any confidence the numerical value of [D.sup.1,2], however, as it is demonstrated in the main body of the paper, the exact value of the intrinsic second order loss is not all that critical and does not affect principal results of our analysis, providing that this parameter is small in a sense. For the purposes of model calculations, we adopt [D.sub.1,2] to be somewhat less than 4 x [10.sup.-5] [cm.sup.3]/s.
Presumably, Hofelich-Abate and Hofelich  were first to introduce (in general form) an intensity-dependent self-limitation term in the photon-density rate equation. Later, that has been done in an explicit form of the second order loss [40, 41]. Those and subsequent studies [42, 43] showed efficient damping of relaxation oscillations in the presence of intensity-dependent losses.
As to the formal analysis of laser coupling via intensity-loss cross-modulation, very few attempts have been done so far to consider such a mechanism--for example, [44,45]. The former work, however, does not treat intrinsic second order losses.
The works which address loss-coupled modes of a single laser rather than loss-coupled single-mode lasers are much more plentiful and varied, and we refer the reader for the reviews in [46, ch. 8] and [47, ch. 12]. In his seminal paper , Baer studied multimode regime of Nd:YAG (neodymium-doped yttrium aluminum garnet) laser with the intracavity-doubling KTP (potassium titynal phosphate) crystal both experimentally and numerically, and proposed the following rate equations for two coupled longitudinal modes:
[[tau].sub.f][[??].sub.1] = [G.sup.0.sub.1] - ([beta][I.sup.1] + [[beta].sub.12][I.sub.2] + 1)[G.sub.1], (A.10a)
[[tau].sub.f][[??].sub.2] = [G.sup.0.sub.2] - ([beta][I.sub.2] + [[beta].sub.21][I.sub.1] + 1)[G.sub.2], (A.10b)
[[tau].sub.e][[??].sub.1] = ([G.sub.1] - [[alpha].sub.1] - [epsilon][l.sub.1] - 2[epsilon][I.sub.2])[I.sub.1], (A.10c)
[[tau].sub.e][[??].sub.2] = ([G.sub.2] - [[alpha].sub.2] - [epsilon][I.sub.2] - 2[epsilon][I.sub.1]) [I.sub.2], (A.10d)
where [G.sub.i] and [I.sub.i] (i = 1,2) are the respective gain and intensity of ith mode, [[tau].sub.f] is the fluorescence lifetime, [G.sup.0.sub.i] is the small-signal gain (pump parameter) for ith mode, [beta] is the saturation parameter which determines how strongly the intensity depletes the available gain, [[beta].sub.12] = [[beta].sub.21] is the cross-saturation parameter for modes 1 and 2, [[tau].sub.e] is the cavity roundtrip time, [[alpha].sub.i] is the loss of ith mode, and e is the nonlinear coupling coefficient, which models the intracavity-doubling crystal as an intensity-dependent loss in the laser resonator.
The last two terms in (A.10c) and (A.10d) represent second order losses that are due to intracavity second-harmonic generation and sum-frequency generation, respectively. Numerical calculations revealed antiphase synchronous oscillations in (A.10a)-(A.10d). In antiphase state, either mode has precisely the same time profile being shifted by 1/2 of a period from its counterpart. This type of dynamics was later observed in a multimode [Nd.sup.3+]:YAG laser with intracavity doubling crystal . For n antiphase oscillators, the phase shift between any nearest neighbors is 1/n of a period.
At [[beta].sub.12] = [[beta].sub.21] = 0, system (A.10a)-(A.10d) is essentially identical to (A.9a)-(A.9d), except that intensity-dependent intrinsic and cross-losses in (A.9c) and (A.9d) are allowed to be independent.
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (A.11)
turns (A.9a)-(A.9d) into nondimensional form (A.7a)-(A.7d). Since the lasers are made of the same material, we can put [d.sub.1] = [d.sub.2], [[GAMMA].sub.1] = [[GAMMA].sub.2], [v.sub.g1] = [v.sub.g2], [a.sub.1] = [a.sub.2], [n.sub.01] = [n.sub.02], [[tau].sub.e1] = [[tau].sub.e2], [[tau].sub.p1] = [[tau].sub.p2], and [D.sub.1] = [D.sub.2], whence [beta] = 1, [[epsilon].sub.1] = [[epsilon].sub.2] = [epsilon], and [[delta].sub.1] = [[delta].sub.2] = [delta], and we arrive, dropping the prime at t, at the set of (1a)-(1d).
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.
 A. Balanov, N. Janson, D. Postnov, and O. Sosnovtseva, Synchronization: From Simple to Complex, Springer, Berlin, Germany, 2009.
 F. C. Hoppensteadt and E. M. Izhikevich, Weakly Connected Neural Networks, vol. 126 of Applied Mathematical Sciences, Springer, New York, NY, USA, 1997
 A. Pikovsky, M. Rosenblum, and J. Kurths, Synchronization: A Universal Concept in Nonlinear Sciences, Cambridge University Press, Cambridge, Mass, USA, 2001.
 S. Strogatz, Sync: The Emerging Science of Spontaneous Order, Hyperion, New York, NY, USA, 2003.
 J. Vandermeer, "Oscillating populations and biodiversity maintenance," BioScience, vol. 56, no. 12, pp. 967-975, 20 06.
 S. Smale, "A mathematical model of two cells via Turing's equation," in Some Mathematical Questions in Biology V, J. D. Cowan, Ed., vol. 6 of Lectures on Mathematics in the Life Sciences, pp. 15-26, American Mathematical Society, Providence, RI, USA, 1974.
 Y. Loewenstein, Y. Yarom, and H. Sompolinsky, "The generation of oscillations in networks of electrically coupled cells," Proceedings of the National Academy of Sciences of the United States of America, vol. 98, no. 14, pp. 8095-8100, 2001.
 A. Gomez-Marin, J. Garcia-Ojalvo, and J. M. Sancho, "Self-sustained
spatiotemporal oscillations induced by membranebulk coupling," Physical Review Letters, vol. 98, no. 16, Article ID 168303, 2007.
 I. Szatmari and L. O. Chua, "Awakening dynamics via passive coupling and synchronization mechanism in oscillatory cellular neural/nonlinear networks," International Journal of Circuit Theory and Applications, vol. 36, no. 5-6, pp. 525-553, 2008.
 W. W. Murdoch, C. J. Briggs, and R. M. Nisbet, Consumer-Resource Dynamics, Princeton University Press, Princeton, NJ, USA, 2003.
 J. D. Murray, Mathematical Biology: I. An Introduction, Springer, New York, NY, USA, 3rd edition, 2002.
 G. P. Agrawal and N. K. Dutta, Semiconductor Lasers, Van Nostrand Reinhold, New York, NY, USA, 2nd edition, 1993.
 W.-B. Zhang, Synergetic Economics: Time and Change in Nonlinear Economics, vol. 53 of Springer Series in Synergetics, Springer, Berlin, Germany, 1991.
 M. V. Volkenstein, General Biophysics, Academic Press, New York, NY, USA, 1983.
 D. S. Chernavskii, E. K. Palamarchuk, A. A. Polezhaev, G. I. Solyanik, and E. B. Burlakova, "A mathematical model of periodic processes in membranes (with application to cell cycle regulation)," BioSystems, vol. 9, no. 4, pp. 187-193, 1977
 G. F. Gause and A. A. Witt, "Behavior of mixed populations and the problem of natural selection," The American Naturalist, vol. 69, no. 725, pp. 596-609, 1935.
 R. MacArthur, "Species packing and competitive equilibrium for many species," Theoretical Population Biology, vol. 1, no. 1, pp. 1-11, 1970.
 T. Baer, "Large-amplitude fluctuations due to longitudinal mode coupling in diode-pumped intracavity-doubled Nd :YAG lasers," Journal of the Optical Society of America B: Optical Physics, vol. 3, no. 9, pp. 1175-1180, 1986.
 T. Erneux and P. Mandel, "Minimal equations for antiphase dynamics in multimode lasers," Physical Review A, vol. 52, no. 5, pp. 4137-4144, 1995.
 A. N. Tikhonov, "Systems of differential equations containing small parameters in the derivatives," Matematicheskii Sbornik, vol. 31(73), no. 3, pp. 575-586, 1952.
 J. Sanders, F. Verhulst, and J. Murdock, Averaging Methods in Nonlinear Dynamical Systems, vol. 59 of Applied Mathematical Sciences, Springer, New York, NY, USA, 2nd edition, 2007.
 F. Verhulst, Methods and Applications of Singular Perturbations: Boundary Layers and Multiple Timescale Dynamics, vol. 50 of Texts in Applied Mathematics, Springer, New York, NY, USA, 2005.
 G. Kirlinger, "Permanence in Lotka-Volterra equations: linked prey-predator systems," Mathematical Biosciences, vol. 82, no. 2, pp. 165-191, 1986.
 R. M. May and W. J. Leonard, "Nonlinear aspects of competition between three species," SIAM Journal on Applied Mathematics, vol. 29, no. 2, pp. 243-253, 1975.
 J. Vandermeer, "Intransitive loops in ecosystem models: from stable foci to heteroclinic cycles," Ecological Complexity, vol. 8, no. 1, pp. 92-97, 2011.
 P. Chesson, "MacArthur's consumer-resource model," Theoretical Population Biology, vol. 37, no. 1, pp. 26-38, 1990.
 S. Mirrahimi, B. Perthame, and J. Y. Wakano, "Direct competition results from strong competition for limited resource," Journal of Mathematical Biology, vol. 68, no. 4, pp. 931-949, 2014.
 M. Devetter and J. Seda, "The relative role of interference competition in regulation of a Rotifer community during spring development in a eutrophic reservoir," International Review of Hydrobiology, vol. 93, no. 1, pp. 31-43, 2008.
 A. D. Bazykin, Nonlinear Dynamics of Interacting Populations, vol. 11 of World Scientific Series on Nonlinear Science, Series A, World Scientific, River Edge, NJ, USA, 1998.
 Y. Kuang, W. F. Fagan, and I. Loladze, "Biodiversity, habitat area, resource growth rate and interference competition," Bulletin of Mathematical Biology, vol. 65, no. 3, pp. 497-518, 2003.
 F. M. Stewart and B. R. Levin, "Partitioning of resources and the outcome of interspecific competition: a model and some general considerations," American Naturalist, vol. 107, no. 954, pp. 171-198, 1973.
 D. Herbert, R. Elsworth, and R. C. Telling, "The continuous culture of bacteria; a theoretical and experimental study," Journal of General Microbiology, vol. 14, no. 3, pp. 601-622, 1956.
 M. Lynch, "Complex interactions between natural coexploiters--Daphnia and Ceriodaphnia," Ecology, vol. 59, no. 3, pp. 552-564, 1978.
 Z. Xiang and X. Song, "Extinction and permanence of a two-prey two-predator system with impulsive on the predator," Chaos, Solitons & Fractals, vol. 29, no. 5, pp. 1121-1136, 2006.
 D. Tilman, Resource Competition and Community Structure, Princeton University Press, Princeton, NJ, USA, 1982.
 H. Statz and G. A. deMars, "Transients and oscillation pulses in masers," in Quantum Electronics, C. H. Townes, Ed., pp. 530-537, Columbia University Press, New York, NY, USA, 1960.
 A. Ghosh, B. K. Goswami, and R. Vijaya, "Nonlinear resonance phenomena of a doped fibre laser under cavity-loss modulation: experimental demonstrations," Pramana-Journal of Physics, vol. 75, no. 5, pp. 915-921, 2010.
 J. Tian, Y. Yao, J. J. Xiao, X. Yu, and D. Chen, "Tunable multiwavelength erbium-doped fiber laser based on intensity-dependent loss and intra-cavity loss modulation," Optics Communications, vol. 285, no. 9, pp. 2426-2429, 2012.
 E. Hofelich-Abate and F. Hofelich, "Phenomenological theory of damped and undamped intensity oscillations in lasers," Zeitschriftfur Physik, vol. 221, no. 4, pp. 362-378, 1969.
 C. J. Kennedy and J. D. Barry, "Stability of an intracavity frequency-doubled Nd:YAG laser," IEEE Journal of Quantum Electronics, vol. 10, no. 8, pp. 596-599, 1974.
 A. Kellou, K. Ait-Ameur, and F. Sanchez, "Suppression of undamped oscillations in erbium-doped fibre lasers with intensity-dependent losses," Journal of Modern Optics, vol. 45, no. 9, pp. 1951-1956, 1998.
 F. Li, J. N. Kutz, and P. K. A. Wai, "Characterizing bifurcations and chaos in multiwavelength lasers with intensity-dependent loss and saturable homogeneous gain," Optics Communications, vol. 285, no. 8, pp. 2144-2153, 2012.
 A. El Amili, G. Kervella, and M. Alouini, "Experimental evidence and theoretical modeling of two-photon absorption dynamics in the reduction of intensity noise of solid-state Er:Yb lasers," Optics Express, vol. 21, no. 7, pp. 8773-8780, 2013.
 B. A. Nguyen and P. Mandel, "Stability of two loss-coupled lasers," Journal of Optics B: Quantum and Semiclassical Optics, vol. 1, no. 3, pp. 320-324, 1999.
 A. Mustafin, "Two mutually loss-coupled lasers featuring astable multivibrator," Physica D: Nonlinear Phenomena, vol. 218, no. 2, pp. 167-176, 2006.
 P. Mandel, Theoretical Problems in Cavity Nonlinear Optics, Cambridge University Press, Cambridge, Mass, USA, 1997.
 T. Erneux and P. Glorieux, Laser Dynamics, Cambridge University Press, Cambridge, Mass, USA, 2010.
 K. Wiesenfeld, C. Bracikowski, G. James, and R. Roy, "Observation of antiphase states in a multimode laser," Physical Review Letters, vol. 65, no. 14, pp. 1749-1752, 1990.
Department of General and Theoretical Physics, Kazakh National Technical University, 22 Satpayev Street, Almaty 050013, Kazakhstan
Correspondence should be addressed to Almaz Mustafin; firstname.lastname@example.org
Received 9 June 2014; Revised 17 July 2014; Accepted 28 July 2014; Published 10 September 2014
Academic Editor: Qingdu Li
TABLE 1: Existence and stability conditions of nonnegative equilibria of system (1a)-(1d). Equilibrium Existence Stability F Always Unstable [F.sub.1] Always [[gamma].sub.2]/ [[gamma].sub.1] < [[ALEPH].sub.1] [F.sub.2] Always [[gamma].sub.2]/ [[gamma].sub.1] > 1/[[ALEPH].sub.2] [F.sub.12] 1/[[ALEPH].sub.2] < Unstable [[gamma].sub.2]/ [[gamma].sub.1] < [[ALEPH].sub.1] for [[ALEPH].sub.1,2] > 1 (strong coupling) [[ALEPH].sub.1] < [[ALEPH].sub.1,2] = [[gamma].sub.2]/ o([[epsilon].sup.1/2]) [[gamma].sub.1] < 1/[[ALEPH].sub.2] for [[ALEPH].sub.1,2] < 1 (weak coupling) TABLE 2: Laser parameters used for numerical simulations. The values in parentheses stand for associated nondimensional quantities. Notation Quantity Meaning Value J ([gamma]) Pump current density 5 x [10.sup.3] A/[cm.sup.2] (1.19375) d Active layer thickness 2 x [10.sup.-5] cm [GAMMA] Confinement factor 0.3 [[mu].sub.g] Group refraction index 4 a Differential gain 2.5 x [10.sup.-16] coefficient [cm.sup.2] [n.sub.0] Carrier density at 1 x [10.sup.18] transparency [cm.sup.-3] [[tau].sub.e] Carrier lifetime 2.2 x [10.sup.-9] s [[tau].sub.p] Photon loss time 1.6 x [10.sup.-12] s [epsilon] [[tau].sub.p]/ 0.727273 x [10.sup.-3] [[tau].sub.e] D ([delta]) Intrinsic second 0.773 x [10.sup.-5] order loss [cm.sup.3]/s (0.01)
|Printer friendly Cite/link Email Feedback|
|Title Annotation:||Research Article|
|Publication:||Journal of Applied Mathematics|
|Date:||Jan 1, 2014|
|Previous Article:||A Crank-Nicolson difference scheme for solving a type of variable coefficient delay partial differential equations.|
|Next Article:||Fair optimization of video streaming quality of experience in LTE networks using distributed antenna systems and radio resource management.|