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.
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, , such that the reaction equations for the set of metabolites, , 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 , where is a diagonal matrix with the elements , 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:Here dots indicate differentiation with respect to the nondimensional time variable , and () are quantities measuring the respective population sizes of th resource and th consumer, () is the inflow rate of th resource, is a parameter representing consumer self-limitation, ( and ) quantifies the inhibitory effect of th consumer on the growth of th consumer due to coupling, and 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, directly represents population density of consumer species and is always nonnegative. Quantity , 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 , where and 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) corresponds to zero population size in reality. Nevertheless, from here on we shall apply the term “resource” to 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 , the communities are uncoupled and completely independent. An isolated CR pair is governed by equations:
There exist two nonnegative steady states:The linearization of (2a) and (2b) takes the form where (, ) is one of the above steady states (3a) and (3b).
At (3a), one eigenvalue is negative and one is positive: , . Thus (3a) is a saddle point. At (3b), and , so the steady state is a stable node/focus. Specifically, focus is the case for whence one obtains an asymptotic estimate 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.
(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, , , and , involved in the CR system (2a) and (2b) when it is overdamped, and two— and —when underdamped.
System (2a) and (2b) is singularly perturbed because the derivative of one of its state variables, , 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.
Replacing in (2a) and (2b) with a fast time variable and setting we obtain the fast subsystem:where prime means differentiation with respect to . In the stretched timescale the slow resource variable , 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 valid for .
After a lapse of considerable time (in the fast scale) converges to either of two fixed points depending on a sign of : It means that every single trajectory starting within the positive quadrant of plane far enough from the stable steady state (3b) will run almost parallel to the vertical axis and hit the line practically in a finite time of order . This is shown in Figures 1(a) and 1(b).
Now set in (2a) and (2b) to get the slow subsystem:Equation (10b) describes a slow manifold consisting of two lines in the plane: and . 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), , is an isolated root of (10b) and is a stable solution of (10b) for any ; the assumptions of Tikhonov’s theorem [20] are satisfied, and one may proceed to approximate and 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 , 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 .
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 (), and one slow (): Since , trajectories starting off the associated eigenvector 2 (which is tangent to the slow manifold ) converge to that vector along lines almost parallel to eigenvector 1 (which is parallel to the vertical axis ). 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 .
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 (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 : 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 .
To for small there is a way to recast (2a) and (2b) near the focal equilibrium in a convenient form where the perturbation parameter does not multiply any right hand side. Namely, performing the scaling , , and , where , one obtains Here dots stand for differentiation with respect to the time variable .
Equations (13) represent a weakly perturbed Hamiltonian system. For the system is pure Hamiltonian and admits the first integral which is a conserved quantity (). The periodic solutions of the Hamiltonian system form a one-parameter family with the equilibrium as center point. The condition 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 only: This is the equation for a nonlinear quasi-conservative oscillator. As is seen from (15), for and the oscillator is not only quasi-conservative, but also quasi-linear with the frequency . These conditions motivate introducing the new variable . As a result one obtains 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, , and weak, .
Physically feasible equilibria, or fixed points, , of (1a)–(1d) are those for which . We denote the interior fixed point by , 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 and designate either of one-consumer equilibria corresponding to dominance, while means both consumers having been washed out.
Model (1a)–(1d) has four feasible steady states. To for small where
Boundary equilibria , , and always exist. exists if or In case (19) the expression with “+” in (17d) is realized, while in case (20) the one with “−” is feasible.
is always unstable, because two of the four associated eigenvalues are positive:
Correct to in , the eigenvalues for and arewhence it follows that (17b) and (17c) are stable whenrespectively.
The necessary and sufficient conditions for all the eigenvalues of the Jacobian matrix, evaluated at , to have negative real parts are, from the Routh-Hurwitz criterion, where , , , and are the coefficients of the characteristic polynomial of (24):
To analyze the validity of (25a)–(25d) we put , where . Then for (25a) to be true, 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 to have negative real part. For it boils down to whence it follows that 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 —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 and are stable with an unstable coexistence steady state , 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, , while the inflow of resource 1, , steadily increases from a value less than along the line in the parameter plane, as shown in Figure 3(a). Referring to (17a)–(17d) and (23a) and (23b), one sees that initially is unstable; is stable and does not exist with the complete dominance of consumer 2. Within the interval steady state becomes stable yet empty and becomes existent (according to (19)) yet unstable, so the situation remains unchanged until point in Figure 3(a) has been reached from the left. For a larger , state gives up its stability and the system jumps to . Consumer 2 gets washed out, while consumer 1 takes over. If now we start reducing , the system remains in until drops to the lower critical value , beyond which is no longer stable and there is a reverse jump to . In other words, as progresses along there is a discontinuous switch from consumer 2 to consumer 1 at , while as retraces its steps, there is a discontinuous switch from consumer 1 to consumer 2 at . 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 . 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.
(a)
(b)
3.4. Weak Coupling: Antiphase Relaxation Oscillations
As seen from Table 1, the very existence of interior equilibrium in a case of weak coupling (condition (20)) implies instability of both boundary fixed points, and . System (1a)–(1d) happens to possess four nonnegative steady states, none of them being stable. As we have found, 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 , where . The new dynamic variables and are defined by formulas and . Thus the variables and are the respective deviations of and from their standalone steady-state values. With these new variables equations (1a)–(1d) becomewhere dots mean differentiation with respect to and are the scaled coupling strengths.
For brevity, we restrict our consideration to the case of identical CR pairs and symmetric coupling setting and . 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 and . Secondly, (29b) and (29d) are differentiated with respect to time to get and . Thirdly, , , , and are plugged into the equations for and . And finally, are replaced by the new variables : .
As a consequence, we arrive at the following coupled equations: where
We seek the solution of (30) in the quasi-harmonic form where and are slowly varying amplitudes and phases, and is the frequency of the synchronous oscillations. Differentiating the assumed form of and equating the result to the assumed form of yields the first pair of relationships between and :
Then, differentiating and substituting the resulting expression for as well as the assumed forms for and into (30) yields the second pair of equations relating and . Separating into equations for the rate of change of and one obtains where
To this point no approximations have been made except of the expansion in powers of in (30). Averaging equations (34) over the period and considering , , and to be constants while performing the averaging, one obtains the following equations describing the slow variations of and : where .
Any stable fixed points of (36a) and (36b) could mean that the phase difference between the coupled oscillators do not change in time (), and the oscillations are periodic with constant amplitudes . 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 .
Let be the vector of the right-hand sides of (34). Also denote the operation of averaging by the angle brackets . Then will mean the right-hand sides of (36a) and (36b). Following the conventional procedure we retain only vibrational terms in , and integrate up to an arbitrary function of the amplitudes chosen for simplicity to be equal to zero: Now, upon calculation of the Jacobian matrix , we can write down the second approximation symbolically as (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 enabling us to present the resulting equations of the second approximation as compact as possible:
The stable steady-state solution for the phase difference is . Inserting this value into (40a) and putting we get three different steady-state solutions for (including the trivial one), of which only one, 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.
Figure 4 plots phase planes for coupling strength below and above the Hopf bifurcation. When (consumers are decoupled) the internal equilibrium is a stable focus (Figure 4(a)). For (very weak coupling; just below the bifurcation) is still a stable focus, though a very gently winding one: the decay is slow (Figure 4(b)). For (very weak coupling; just above the bifurcation) there is an unstable focus at 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)).
(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, , meaning that system (1a)–(1d) is singularly perturbed. The slow variables are resources, and , and the fast variables are consumers, and . 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 . Now rescale (1a)–(1d) by replacing with and, after taking , it becomeswhere prime means differentiation with respect to . This is the fast subsystem, where and are replaced by their initial values and treated as parameters. It yields inner solution, valid for .
Setting in (1a)–(1d) leads to the slow subsystemwhich produces outer solution, valid for . In this singular limit as , the subsystem defines a slow flow on the surface (slow manifold) given by (43c) and (43d). Outer solution is valid for those and , 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 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 , 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 (the slow variables are deemed to be frozen):Existence and stability of these quasi-steady states are determined by —position of a representing point in the phase plane of the slow subsystem, shown in Figure 5(a).
(a)
(b)
(c)
(d)
, and always exist for all and from the positive quadrant of the slow phase plane. For not-too-weak coupling, such that , exists for all and satisfying the condition , that is, within the opening of the angle formed by lines and in Figure 5(a). The opening shrinks as the coupling strengths get weaker.
Jacobian matrix of the fast subsystem, has the following sets of eigenvalues at (44a)–(44d):Based on (46a)–(46d) one concludes that (the origin) is always an unstable node for all and from the positive quadrant of the slow phase plane. is a stable node for , that is, below the line in Figure 5(a), otherwise it is a saddle. Similarly, is a stable node for , that is, above the line in the plane of slow variables, otherwise it is a saddle. The interior quasi-steady state, , 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 is stable and is not. Consumer 1 completely dominates. This corresponds to slow variables and being somewhere below the line of Figure 5(a). Fast subsystem (42a)–(42c) has phase portrait of a type shown in Figure 5(b). While remains much greater than , the dynamics of the resources (treated as bifurcation parameters in reference to the consumers) is described by a system of two independent equations:which is a piecewise version of the slow subsystem (43a)–(43d) for and . System (47a) and (47b) has stable steady state where . This equilibrium lies in the upper left corner of Figure 5(a).
While heading to (48), the trajectory crosses the line and enters the domain of bistability of both and . 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 and from the respective steady states (48), the system (47a)-(47b) becomewhence one finds
It follows from (50), that the dynamics of variable is faster than that of due to small . Clearly, the representing point must have relaxed to the vertical line well before approaching the horizontal line . 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 (slow manifold of the system (47a)-(47b)), steadily tending to point (48) (Figure 6(a)).
(a)
(b)
If is small in terms of the condition (6), then the individual CR pairs are underdamped. In the immediate vicinity of 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: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 these oscillations manifest themselves in damped transverse fluctuations superimposed on the independent vertical motion along the nullcline () toward point (48) (Figure 6(b)).
By virtue of the condition (20), the trajectory has to cross the line on its way toward the neighborhood of steady state (48). As soon as this has happened, node in the plane will be absorbed by saddle . 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 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:with the initial conditions and .
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: where . This equilibrium lies in the lower right corner of Figure 5(a).
Variable , being more rapid in comparison to , quickly enters the neighborhood of the nullcline given by and then—depending on the value of —finally approaches the nullcline either monotonically (Figure 6(a)) or via damped oscillations according to equations(Figure 6(b)). System (54a)-(54b) describes underdamped intrinsic oscillations of uncoupled consumer 2 for small .
At the same time, steadily and independently tends to . Again, because point (53) is located below the line (on the strengths of the condition (20)), the trajectory would certainly cross that line at a point , whereupon node would be absorbed by saddle . 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.
(a)
(b)
(c)
(d)
The resources and 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, . These times are determined mainly by the dynamics of the slow resource variables and and, to a zeroth approximation in and , can be found as solutions to the equations of motion (52a) and (47b) with respective boundary conditions and . In this way one obtains a quite simple estimate for the period: It is interesting that, according to (55), the period depends on the ratio of the two resource inflows, and , rather than on each of them individually, and does not depend completely on concrete value of .
The consumers and change periodically between extinction and respective constant levels and . 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 , 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 , for example, rises above a threshold value (determined by the amount of losses experienced by ), then comes into dominance causing 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) is essentially zero, while is approximately equal to its uncoupled steady-state value, . increases due to its inflow until it overcomes losses for consumer 1;(2)with a sufficient resource stock, 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 drops. Cross-losses cause to wash out;(3) and relax to quasi-steady-state values, as if there were only one isolated CR pair. is essentially zero. is increasing, like did in part 1;(4) surpasses the losses, emerges and the subsequent cross-losses cause to wash out. The spiking 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 equilibrium , where only consumer 1 is present, to the neighborhood of , where consumer 2 completely dominates, back to , 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 (Figure 8). 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.
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 and . 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.
Appendix
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]: Here, dots indicate differentiation with respect to time , represents the total biomass of th resource (prey), stands for the total biomass of th consumer (predator) species, the constant defines the growth rate of th resource, is the carrying capacity of th resource, is the rate of uptake of a unit of th resource by each individual of th consumer population, is the conversion efficiency parameter representing an amount of th resource an individual of th consumer population must consume in order to produce a single new individual of that species, and is the loss rate of th 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 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 [16]: where
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 [27].
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 of (A.3a) merely add to interference coefficients which are to be present primordially in (A.1b).
Most mathematical models dealing with coupled CR pairs or multilevel trophic chains ignore contributions of intraspecific and interspecific interference. Indeed, the empirical data like [28] do indicate that may be negligible in comparison with . 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], [30].
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:
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 [31], 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 th resource () flows into a defined volume with the rate while unused resource flows out with the per capita rate , in much the same manner as in a chemostat [32].
In more exact terms, the true chemostat model for one substrate and one species looks as follows: where the rate of substrate uptake is expressed by the Monod formula , is a saturation constant numerically equal to the substrate concentration at which the uptake rate is half the maximum, is the dilution rate defined as the rate of flow of medium over the volume of the bioreactor, and is an input concentration of the substrate.
Model (A.5) turns into an uncoupled version of (A.4a)–(A.4d) if we put and and assume , so that , where .
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 () measures direct interference of individuals within th consumer population with each other resulting in an additional per capita loss rate ; interspecific competition strength (; ) quantifies direct interference effect from th consumer on th consumer resulting in an additional per capita loss rate, , 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 [33]. 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 [23] describes the model in which each predator specializes on one prey only, while Xiang and Song [34] 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 equations (A.4a)–(A.4d) take the following nondimensional form:Note that in (A.7a)–(A.7d) dots mean differentiation with respect to nondimensional “slow” timescale variable , as defined by (A.6).
The parameters , , and reflect the rapidity of the dynamics of , , and with reference to that of . It is assumed that and .
In studying the effect of coupling, the parameters of interest are obviously the coupling strengths, and . The parameters of interest are also those which characterize the difference between the states of the uncoupled systems. The resource income rates and 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 , , and and also drop the prime at , to obtain (1a)–(1d).
A.2. Laser Dynamics Perspective
Laser rate equations originally proposed by Statz and deMars [36] are differential equations that relate two quantities: injected carrier density () and photon density (). For a single-mode semiconductor laser, these equations take the form [12, ch. 6]: where dots mean differentiation with respect to time , is the injection current density (pump parameter), is the magnitude of the electron charge, is the active-layer thickness, is the confinement factor accounting for the fraction of the light power contained in the active region, is the group velocity of light that can be expressed through the speed of light in vacuum () and the group refractive index of the dispersive semiconductor material () as , is the gain coefficient, is the carrier density at transparency corresponding to the onset of population inversion, is the lifetime of the electrons in the conduction band before being lost by escape from the active region, and 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, , of the two has a total loss represented by the sum of the constant loss, , plus the loss proportional to its own intensity, , plus the loss proportional to the intensity of the other laser, : where () is the second order loss constant of the isolated th laser and is the coupling strength measuring the cross-loss effect of th laser on th 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 and , and the coupling strengths and , as free parameters of the model. For the present, we cannot judge with any confidence the numerical value of , 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 to be somewhat less than .
Presumably, Hofelich-Abate and Hofelich [39] 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 [18], 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: where and () are the respective gain and intensity of th mode, is the fluorescence lifetime, is the small-signal gain (pump parameter) for th mode, is the saturation parameter which determines how strongly the intensity depletes the available gain, is the cross-saturation parameter for modes 1 and 2, is the cavity round-trip time, is the loss of th mode, and 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 of a period from its counterpart. This type of dynamics was later observed in a multimode :YAG laser with intracavity doubling crystal [48]. For antiphase oscillators, the phase shift between any nearest neighbors is of a period.
At , 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.
The scaling turns (A.9a)–(A.9d) into nondimensional form (A.7a)–(A.7d). Since the lasers are made of the same material, we can put , , , , , , , and , whence , , and , and we arrive, dropping the prime at , 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.