Linear Augmentation for Stabilizing Stationary Solutions: Potential Pitfalls and Their Application

Linear augmentation has recently been shown to be effective in targeting desired stationary solutions, suppressing bistablity, in regulating the dynamics of drive response systems and in controlling the dynamics of hidden attractors. The simplicity of the procedure is the main highlight of this scheme but questions related to its general applicability still need to be addressed. Focusing on the issue of targeting stationary solutions, this work demonstrates instances where the scheme fails to stabilize the required solutions and leads to other complicated dynamical scenarios. Examples from conservative as well as dissipative systems are presented in this regard and important applications in dissipative predator—prey systems are discussed, which include preventative measures to avoid potentially catastrophic dynamical transitions in these systems.


Introduction
Studies on coupled nonlinear systems have explored a wide variety of emergent dynamical phenomena, namely synchronization [1], oscillator suppression [2], multistability [3], hysteresis [4], extreme-events [5,6] etc. which can be exploited in applications: to model natural phenomena or in regulating the system behavior for instance. Controlling dynamical systems towards a desired behavior is an important research topic in nonlinear sciences [7]. Starting with chaos control [8][9][10][11], research in this domain now also extends towards control of multistability [12], patterns and spatio-temporal chaos [13,14], noisy systems [15,16], methods of stabilizing unstable stationary solutions [17,18] etc. A better understanding of these different regulatory aspects has greatly contributed towards development of related novel and highly efficient procedures. Considering noninvasive (without changing the intrinsic system parameters) mechanisms leading to stabilization of stationary solutions, oscillator suppression via coupling nonlinear systems has been discussed extensively in literature (see Refs. [17,18] for detailed reviews). This suppression is majorly observed as a consequence of parameter heterogeneity between coupled units [19][20][21], presence of time-delayed [22,23]/conjugate variables [24] in the coupling function or through dynamic coupling [25], etc.
Recently, linear augmentation has been suggested as another practical alternative leading to oscillator suppression, which is achieved by coupling systems to a linear feedback which simply consists of an exponentially decaying function [26] in the uncoupled state. Interestingly, the coupling structure of linear augmentation is quite reminiscent of indirect or environmental coupling procedures [27][28][29] which are motivated by observations of collective behaviors in several real world systems, namely, behavior of chemical relaxation oscillators globally coupled through the concentration of chemicals in a common solution [30], dynamics of multi-cell systems where the cells interact through common complex proteins [31], and collective behavior of cold atoms in the presence of a coherent electromagnetic field and atomic recoil [32,33] for instance. These instances therefore also serve as good examples of systems where linear augmentation can exist naturally. Lately, studies have also effectively used linear augmentation in controlling bistability [34], in controlling the dynamics of drive response systems [35] and in controlling hidden attractors [36].
With respect to stabilizing stationary solutions, Ref. [26] discusses results for an augmented Lorenz oscillator where either the stationary solutions of the Lorenz system (desired) or those of the augmented system could be stabilized by picking an appropriate feedback function. The paper also presents some parameter space plots highlighting the regimes where linear augmentation works in stabilizing the solutions of the original system and where it does not. These results although instructive, are also extremely system specific. At this point one must question the ability of linear augmentation towards stabilizing the desired stationary solutions in a more general sense, namely, for which systems, parameter values and coupling configurations will the scheme work? Since the stability of stationary solutions is determined by the eigenvalues of the augmented system's Jacobian in the linear approximation, if the eigenvalues are independent of augmentation, or if they remain positive for all values of the augmentation strength, then the procedure will never stabilize these stationary solutions. Furthermore, linear augmentation introduces new stationary solutions in the dynamics which might get stabilized instead of the intended solutions based on similar arguments. In this paper, we will look at some simple examples of linearly augmented systems to highlight that we need to be quite careful before choosing linear augmentation in such applications. These examples illustrate that there could be situations where even picking an appropriate feedback function does not guarantee that the required stationary solutions will be necessarily stabilized; since the mechanism is highly dependent on the intrinsic properties of the oscillators in consideration, the stationary solutions to stabilize, and also on how these systems are coupled to the feedback/ augmented. More importantly, we will also discuss instances where the failures of this procedure can be exploited in meaningful applications, especially in predator-prey systems where we can potentially avoid catastrophic transitions using linear augmentation.
The manuscript is arranged as follows: Sec. Methods briefly introduces the linear augmentation scheme. In Sec. Results, we will look at augmented Harmonic oscillator, and Duffing oscillator as examples of augmented conservative systems followed by augmented Dissipative predator-prey models where the scheme fails to stabilize the desired stationary solutions, and also discuss possible applications for these observations in the latter. Details regarding certain dynamical aspects of harmonic oscillator and Duffing system are provided in Sec. Additional details for completeness. The results and outlook of this work are summarized in Conclusions.

Linear augmentation
General representation of a linearly augmented dynamical system is, where the column vector .] T corresponds to the transpose) contains the systems variables, and u is the augmentation variable. ε ¼ ½ε 1 ; ε 2 ; . . . ; ε N T 2 R N is the column vector with information regarding the coupling strength of the interaction between the dynamical variables and u; augmentation/coupling term corresponding to the i th component x i is ε i u 8 i = 1, 2, . . ., N, and ε i = 0 if x i is not coupled to u. b 2 R N is an arbitrary vector and k is the decay constant [37] which apparently can be used to control the transient time leading to stationary solutions [26].
gives the dot product of the corresponding column vectors.
In the following, we will look at examples of augmented conservative and dissipative dynamical systems which highlight the limitations of this procedure. The terms augmentation/ augmented and coupling/coupled are used synonymously in the following text.

Results
Here we will discuss some instances of systems controlled via linear augmentation. We will first look at two examples of conservative systems, namely the harmonic oscillator and conservative Duffing oscillator where linear augmentation is employed to stabilize their stationary solutions.

Harmonic oscillator
Equations describing a partially linearly augmented harmonic oscillator are: where ω is the frequency of the oscillator, k is an augmentation parameter, and ε is the coupling strength. The first two equations governing the evolution of x and y correspond to the original harmonic oscillator dynamics of position and momentum respectively. In the absence of augmentation, harmonic oscillator conserves total energy, which stays constant on the ellipses shown in Fig 1b. Each of these ellipses correspond to the systems' evolution following different initial values of position and momentum, and hence, different conserved total energies. Harmonic oscillator has x Ã = 0, y Ã = 0 as the only stationary solution and note that the augmentation term only appears in the rate equation of the position variable x at this point. In case of a successful stabilization, the required stationary solution of the full system should be (x Ã , y Ã , u Ã ) = (0, 0, 0) (origin) where the system effectively decouples from the controller. The characteristic eigenvalue equation at the origin for this system is, For ε = 0, the eigenvalues for the full system are λ 1,2 = ±iω corresponding to the non-hyperbolic stationary solution at the origin, and λ 3 = −k corresponds to the decay of the control variable u: u(t) / exp(−kt) in this case. The other parameter values for the following calculations are fixed at ω = 2, and k = 2. For the evolution of the augmented system (ε > 0), the bifurcation diagram of the system with increasing ε values is shown in Fig 1a (black dots). For this calculation, the system is evolved for 10 different initial conditions for each value of ε sufficiently for the transients to be discarded. System is then further evolved to capture the possible dynamical regimes by recording all the extrema during the evolution. It is seen that with an increasing ε, the system which was conservative for ε = 0 becomes dissipative and gets into a stable origin regime even for arbitrarily small values of ε. Rewriting Eq 3 as, and applying the Routh-Hurwitz criteria (RHC) [38], we can deduce that the roots of this equation are all either negative or have negative real parts (in case of complex roots) 8 ε > 0. Largest eigenvalues of the Jacobian (red symbols) are also plotted along with the bifurcation diagram in Fig 1a which demonstrate the transition from oscillatory to stationary state for any ε > 0. Considering the behavior of this system for large ε, we can see that the largest eigenvalue λ = λ 1 ! 0 from Eq 4 in this limit. Since the discriminant (for a general cubic equation of the cubic characteristic Eq 4 is negative 8 ε ! 0, this implies that the system has one real eigenvalue and a pair of complex conjugate eigenvalues in this range. The negative real part of these complex eigenvalues for large ε can therefore be estimated by equating the sum of all eigenvalues to the trace of the Jacobian tr(J), giving Re(λ 2,3 ) = −k/2. This further implies that as ε ! 1, convergence to the origin gets slower although origin is stable in the entire ε > 0 range and any change in stability will only occur as ε ! 1 when λ 1 = 0. Now let us consider a more general case of an augmented harmonic oscillator given by, where the augmentation now appears in the rate equations of both position x and momentum y with coupling strengths ε 1 , ε 2 respectively. The eigenvalue equation in this case is, Substituting ε 2 = 0 and ε 1 = ε in Eq 5 yields the dynamics of Eq 2. Similarly, for ε 1 = 0 and ε 2 = ε, we obtain a case where the system is only coupled in the y variable for which the characteristic Eq 6 is exactly identical to Eq 4, and therefore the stability characteristics of the origin are identical and independent of whether the system is augmented in x or y. In the previous example, we saw a situation where linear augmentation successfully stabilized the origin for the entire range of ε > 0. Now considering ε 1 = ε 2 = ε (the system is similarly augmented in both variables), in which case Eq 6 gives, Using the RHC, it can be checked that this equation will have all negative eigenvalues iff kω 2 and for higher values of ε, RHC suggests appearance of positive eigenvalue/eigenvalues. For large ε, we can get an estimate of largest eigenvalue l 1 ! ðo 2 À 1Þ 2 > 0 8 ω > 1. Since the discriminant is negative 8 ε > 0, the remaining complex conjugate eigenvalue pair have a negative real part given by Reðl 2;3 Þ ¼ ðtrðJÞ À . Therefore, unlike in the previous example, we can see that origin here is unstable for large ε. The expression for ε Ã also shows that a higher value of k can extend the coupling range for a stable origin. This result is the first instance of unexpected behavior as we would normally expect a higher coupling value to keep the origin stable. Furthermore, in the ε > ε Ã regime it is numerically observed that the trajectories escape to infinity which is also quite unexpected.
One of the primary reasons behind considering augmented harmonic oscillator in this study is the fact that it is highly solvable and therefore can provide necessary insights into the physical mechanisms behind the desirable as well as undesirable behaviors. It turns out that in case of a successful stabilization, this system represents a forced harmonic oscillator where the steady state solution (which is completely determined by the forcing) decays to the origin along with an exponentially decaying force. Similarly, the case where the trajectories escape to infinity corresponds again to a forced system but this time being driven by an exponentially diverging force which is analogous to a situation where energy is being pumped into the system. Therefore the steady state solution in this case diverges along with the diverging force explaining the unexpected behavior of escaping trajectories observed for the fully augmented system. Details of the calculations leading to these deductions are available in Sec. Harmonic oscillator: diverging trajectories.

Duffing oscillator
General equations for a linearly augmented Duffing oscillator with no damping or forcing can be written as: Uncoupled Duffing system has an invariant of motion (also the Hamiltonian) given by H(x, y) = y 2 /2 − x 2 /2 + x 4 /4 and stationary solutions: (x Ã , y Ã ) = (±1, 0), (0, 0). The trajectories of this system evolve on the double well potential surface of H(x, y) on starting with different initial conditions for ε 1 = ε 2 = 0. Similar to the previous example, for a successful stabilization, the required stationary solutions of the full system should be (x Ã , y Ã , u Ã ) = (0, 0, 0) or (±1, 0, 0) where the system effectively decouples from the augmentation. These solutions will be referred to as (x Ã , y Ã ) = (0, 0) (origin) or (±1, 0) in the following.
For the system in Eq 8, the characteristic eigenvalue equation can be expressed as, Now similar to the harmonic oscillator example, considering partial augmentation with ε 1 (2) = ε, and ε 2(1) = 0 first, Eq 9 suggests that the stability characteristics for the stationary solutions are again independent of whether the system is being augmented in x or y. For this partial augmentation, Eq 9 gives, Substituting x Ã = 0,±1, and rearranging the terms, we can obtain the characteristic eigenvalue equations for these stationary solutions as, for (x Ã , y Ã ) = (0, 0) (hyperbolic for ε = 0), and for (x Ã , y Ã ) = (±1, 0) (non hyperbolic for ε = 0) respectively. It is straightforward to check that the largest eigenvalue λ 1 ! 0 for larger ε values in both these cases which implies that a stable/ unstable stationary solution will stay the same until a stability change (zero crossing of the eigenvalue/s) occurs in the ε ! 1 limit. Furthermore using the RHC, it is easily verifiable that Eq 11 will always have positive root/roots, whereas Eq 12 will have all negative roots 8 ε > 0; which implies that the x Ã = 0 is always unstable and x Ã = ±1 is always stable. Therefore, we see that partial augmentation works for stabilizing (x Ã , y Ã ) = (±1, 0) but fails completely to stabilize the origin (x Ã , y Ã ) = (0, 0). Fig 2 shows the largest eigenvalue calculations which verify these deductions. This brings us to an important observation that there might exist situations where it is not possible to target the desired stationary solution even on using an appropriate feedback function with any combination of k and ε values. Now considering identical augmentation with ε 1 = ε 2 = ε and we will see that this system has some interesting properties. . It is observed that even for very small coupling values, the system quickly gets into a stable stationary state regime, although, for smaller values of ε, it exhibits bistability. The transient trajectories in this parameter regime are shown in Fig 3(a.1). We observe that the augmentation is stabilizing our desired stationary solution at (x Ã , y Ã ) = (1, 0), but along with it, other stationary solutions which are ε dependent are also getting stabilized on starting with different initial conditions. These other stationary  solutions for the augmented system here are given by, and solutions (x Ã − , y Ã − , z Ã − ) are observed to coexist along with (x Ã , y Ã ) = (1, 0). For higher coupling values, bistability terminates via a saddle node bifurcation when the stable branch of In absence of augmentation, the eigenvalues for (x Ã , y Ã ) = (1, 0) are complex: l 1;2 ¼ AEi ffiffi ffi 2 p . For the augmented system, the characteristic equation can therefore be written as: The RHC shows that this equation will have all negative roots for (2k − ε 2 )>0 and positive root/roots appear for ε > ffiffiffiffiffi 2k p . This gives us the transition threshold for the destabilization of the stationary solution as ε Ã ¼ ffiffiffiffiffi 2k p , at which the eigenvalue/s cross the zero axis. Since the discriminant is negative, the characteristic equation has one real and two complex conjugate roots. Considering large ε behavior, it is seen that the largest eigenvalue λ 1 ! 1/2 which implies that (x Ã , y Ã ) = (1, 0) is unstable in this range. The real part of the remaining complex conjugate eigenvalue pair is Re(λ 2,3 ) = −(2k + 1)/4. At ε Ã , Eq 14 can be rewritten as, which consequently gives the eigenvalues as λ 1 = 0 and l 2; We can see that λ 2,3 will be a complex conjugate pair for k 2 ð8 À 6 ffiffi ffi 2 p ; 8 þ 6 ffiffi ffi 2 p Þ. For our calculations, we have considered k = 2 which gives that at ε Ã ¼ ffiffiffiffiffi 2k p ¼ 2, λ 1 crosses the zero line as can be seen in Fig 3a (red symbols). For higher values of ε > ffiffiffiffiffi 2k p , stationary states (x Ã + , y Ã + , u Ã + ) (circles) in Fig 3(a) for ε > ε Ã (= ε TC ) get stabilized via a transcritical bifurcation where (x Ã , y Ã ) = (1, 0) and (x Ã + , y Ã + , u Ã + ) exchange their stability. This again is quite unexpected since the controller is designed to stabilize (x Ã , y Ã ) = (1, 0) for higher ε values. A brief discussion regarding the behavior of this system in the (ε, k) plane is available in Sec. Duffing system: (ε, k) plane behavior.
For the origin at (x Ã , y Ã ) = (0, 0), the bifurcation diagram (black dots) is shown in Fig 3b. Appropriate transient trajectories corresponding to different augmentation regimes are also shown in related Fig 3(b.1) and 3(b.2). We observe bistability for a range of lower ε values before the origin gets stabilized. The stationary solutions obtained in the bistable regime are given by Transient trajectories in this regime demonstrating the two observed stationary solutions are shown in Fig 3(b.1). These solutions approach and collapse at the origin in a pitchfork bifurcation for ε PF ¼ ffiffi ffi k 2 r (= 1 for k = 2 in this case) beyond which the solutions (x o ± , y o ± , u o ± ) become imaginary and the origin is the only stable real stationary solution. A transient trajectory in this parameter regime is shown in Fig 3(b.2).
The characteristic equation for the origin is, which has all negative eigenvalues for 2ε 2 −k > 0 giving us a stability regime of ε > ffiffiffiffiffiffiffi ffi k=2 p and a transition value of ε PF ¼ ε Ã ¼ ffiffiffiffiffiffiffi ffi k=2 p ¼ 1 (since k = 2) when the eigenvalue/s cross the zero axis. From Eq 17, we get λ 1 = −1 in the large ε limit. It is numerically observed here that the discriminant Δ < 0 in this range and therefore Re(λ 2,3 ) = (1−k)/2 = −0.5 and consequently, the origin is stable in the large ε limit. The largest eigenvalue for the origin is plotted as red symbols in Fig 3(b) which shows the changes in the stability of the origin from unstable in ε 2 (0, 1) to stable 8 ε > 1. For a discussion regarding the system behavior in the (ε, k) plane, please again see Sec. Duffing system: (ε, k) plane behavior for details.
For Since this solution is a symmetric counterpart of (x Ã , y Ã ) = (1, 0), the corresponding analysis similarly carries over in this case.
These simple examples demonstrate the fact that targeting the required stationary solutions using linear augmentation is not quite straightforward and the procedure is quite sensitive to how the systems are augmented, the stationary solutions being targeted and to the properties of systems. In the following, results for a specific class of dissipative dynamical systems are presented to further highlight these limitations.

Dissipative predator-prey models
Considering predator-prey population models, general evolution equations for these systems with logistic prey growth can be written as, where x and y correspond to prey and predator populations respectively and the parameters r, K, ρ, and γ are positive. Considering the evolution equation for preys, the first term rx(1 − x/K) represents the logistic growth rate of the prey species with the maximum growth rate of r and carrying capacity K which is the maximum population size that the environment can sustain indefinitely. The second term f(x)y corresponds to the prey mortality via predation. f(x) is the functional response governing the rate of per capita prey consumption by the predators [39][40][41]. The parameter ρ governs the biomass conversion efficiency for the predators in the sense of how many predators are added to the population via predation, and γ is the intrinsic predator mortality parameter. One of the stationary solutions of this system corresponds to vanishing predator-prey populations, i.e. (x Ã , y Ã ) = (0, 0). The other stationary solutions are dependent on the type of functional response considered. Most commonly employed f(x) forms in such models are the Holling type with the following general expressions:  [45] which is used in modeling phytoplankton and zooplankton interactions leading to harmful algal blooms, and consequently, corresponding stationary solutions can be obtained. The parameter a in expressions above corresponds to the maximum per capita predation rate, and b is the half saturation constant governing how quickly the predators attain their maximum consumption rate. In the following, we will have a closer look at the stability properties of the trivial stationary solution (x Ã , y Ã ) = (0, 0): origin. Considering a general augmented population model, it turns out that the Jacobian for this system is identical for all three functional responses at the origin. The identical characteristic equation therefore is, ðr À lÞðg þ lÞðk þ lÞ þ ε 2 2 ðr À lÞ À ε 1 2 ðg þ lÞ ¼ 0: For ε 1 = ε 2 = 0, we obtain the eigenvalues as λ 1 = r, λ 2 = −γ, and λ 3 = −k where λ 1 and λ 2 are the eigenvalues for the original system in Eq 18 implying that the origin is unstable, and λ 3 corresponds to the exponentially decaying control variable u. Since the Holling type I case with insatiable predators is quite unrealistic, we will focus here on systems with Holling type II (H II) and III (H III) behaviors. In the following analysis, the parameter values are fixed at: r = 0.5, K = 0.5, a = 1/3, b = 1/15, ρ = 0.5, γ = 0.1 for the H II [46] system, and r = 0.43, K = 1, a = 1, b = 0.053, ρ = 0.05, γ = 0.028 for H III [45,47]. Let us now look at different augmentation situations.
For ε 1 = ε and ε 2 = 0, i.e. only prey populations are augmented, substituting these values in Eq 20 gives, ðg þ lÞ½ðr À lÞðk þ lÞ À ε 2 ¼ 0: Since one of the roots λ = −γ is independent of ε therefore the remaining roots of this equation determine the stability of the origin. The remaining two roots are l AE ¼ Àðk À rÞ=2 AE ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðk þ rÞ 2 À 4ε 2 q =2 out of which, λ − < 0, 8 ε. It is easily verifiable that the eigenvalue λ + (which also is the largest) is positive 8 ε < ffiffiffiffi ffi kr p and crosses the zero axis at ε Ã ¼ ffiffiffiffi ffi kr p leading to all negative eigenvalues and hence a stable origin. This is quite similar to the harmonic oscillator case where increasing/decreasing the value of the decay parameter k could increase/ decrease the threshold value of stable ! unstable transition (unstable ! stable in this case). Furthermore, in the large ε limit, we obtain the largest eigenvalue λ 1 = −γ and therefore the origin is stable in this regime. Fig 4: top row shows the bifurcation diagram and the largest eigenvalue behavior for H II (left) and III (right). The unstable ! stable transition in both these systems with increasing coupling values can be seen in the figure. Although for the H III system in the regime ε > ffiffiffiffi ffi rk p ð¼ 0:927; for r = 0.43, k = 2), certain initial conditions lead to the trajectories escaping to infinity (not shown) which accounts for the missing dots in the bifurcation figure. Since x and y are population variables by definition, population models are constrained to work for non-negative values of x and y respectively. What we observe here is that the augmentation forces the prey populations into negative values which leads to a breakdown in the model constraints and the logistic function in the rate equation of x leads to diverging solutions as time increases. For H II system this appears not to be the case and all considered initial conditions lead to a stable origin 8 ε > ffiffiffiffi ffi rk p ð¼ 1; for r = 0.5, k = 2). For ε 1 = 0 and ε 2 = ε, i.e. only the predator populations are augmented, substituting these values in Eq 20 gives, and we see that an eigenvalue λ = r is always positive since r > 0, and therefore this setup will never stabilize the origin. The remaining eigenvalues are l AE ¼ Àðg þ kÞ=2 AE ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ðk À gÞ þ 4ε 2 p =2. In Fig 4: second row, for low ε values, the systems exhibit periodic behavior similar to the one shown for the augmented prey case. For higher values of ε, both H II and H III settle on a stationary solution of the original system (x Ã , y Ã ) = (K, 0) which we did not intent to stabilize. For this solution, the preys exist at their carrying capacity and the predators vanish. For H II and H III, the carrying capacities considered for simulations are K = 0.5 and K = 1 respectively, and hence the observations in Fig 4 (middle row).
For ε 1 = ε 2 = ε, i.e. both prey and predator populations are augmented, substituting these values in Eq 20 and rearranging terms gives, l 3 þ ðg þ k À rÞl 2 þ ð2ε 2 þ ðg À rÞk À rgÞl þ ε 2 ðg À rÞ À rgk ¼ 0: Using the RHC, one of the conditions for this equation to have all negative roots is ε > ffiffiffiffiffiffiffiffiffiffiffiffiffiffi rgk ðg À rÞ s which is impossible to achieve since r > γ. Therefore, this setup will not stabilize the origin either. Fig 4: bottom row shows the behavior of H II and H III. For smaller ε values, these systems exhibit periodic behavior similar to the augmented prey. On increasing ε further, systems enter a regime where all considered initial conditions lead to escaping trajectories. The reason behind this behavior here again is due to a breakdown in modeling constraints. Examination of transient trajectories reveals that augmentation in this case is forcing the predator populations into y < 0 axis which leads to a breakdown in the model, thereby initiating a positive feedback loop in the prey populations leading to the diverging behaviors observed in simulations. Beyond this regime for higher values of ε, H II system exhibits bistability between different stationary solutions where in one case, preys exceed their carrying capacity (x Ã > K) and the predator populations are negative (y Ã < 0), and in the other case, the prey populations are negative (x Ã < 0) and predators assume a small positive value. It is important to note yet again that these solutions are impractical because the populations cannot exist above their carrying capacities nor can they take negative values. For H III system, we only observe the Augmented predator-prey model dynamics. Bifurcation diagram (black dots), largest eigenvalues (red circles) and time series (insets for specific ε values) for predator-prey systems with H II (left column) and H III (right column) functional responses. Top row: For augmented prey, insets show the time series of x, y and u for two different ε values before (with oscillatory u) and after the stabilization of (x*, y*) = (0, 0) (with u* = 0). The Hopf bifurcation is marked as H, and the transcritical bifurcation points have been highlighted by T 1 , and T 2 respectively. Middle row: For augmented predator, the systems exhibit oscillatory behavior similar to augmented preys (not shown) for low ε before they settle on the stationary solution (x*, y*) = (K, 0) with u* = 0; where the preys reach their carrying capacity in the absence of predators for higher ε values. Bottom row: For augmented predator and prey, for low ε, the systems are oscillatory (not shown). With increasing ε, both the systems lose the oscillatory behavior and all trajectories escape to infinity in the regime marked by E. For higher ε values, unrealistic stationary solutions where either the preys exceed their carrying capacity (x* > K) with negative predator populations (y* < 0) (H II and H III), or where the prey populations are negative with small positive predator population (for H II) and u* 6 ¼ 0 get stabilized. equilibrium solutions with x Ã > K and y Ã < 0 (see inset). In both the cases, we have a non vanishing u Ã > 0 and therefore these solutions exist due to augmentation and cannot be observed otherwise. Following this analysis, we can conclude that augmenting the prey is the correct strategy to stabilize of the origin and the other coupling schemes can lead to complicated dynamics. Even though the analysis here is limited to the origin, we can expect these behaviors to be quite general with respect to other stationary solutions as well. Now considering applications, as already mentioned, origin corresponds to an equilibrium for which the predators and the preys vanish. Persistence of populations for a proper ecosystem function is very imperative and has been studied extensively from several perspectives, contributing towards a better understanding of the processes leading to species extinction [48][49][50][51][52]. Knowledge regarding these processes can help in devising procedures which can contribute towards better species conservation efforts. For the simple models considered in the previous analysis, it is clear that either by coupling the system appropriately or by using specific parameter values for k and ε, we can avoid stabilizing the origin. For instance, considering the prey augmented case, for low ε, the systems exhibit periodic oscillations. On increasing the coupling strength, stationary solutions of the augmented system (x Ã > 0, y Ã > 0, u Ã < 0), satisfying ð23Þ get stabilized through a reverse Hopf bifurcation (marked as H in Fig 4 (top row)). For these stationary solutions, the value of x Ã stays constant while y Ã and u Ã = −εx Ã /k show a variation for a range of ε values (plateau between H and T 1 in Fig 4 (top row)). It is also important to note that some initial conditions in this regime can lead to trajectories escaping to infinity. This branch of solutions undergoes a transcritical bifurcation (marked T 1 in Fig 4 (top row)) where it exchanges stability with another branch of solutions with u Ã ! 0 for increasing ε. At ε Ã ¼ ffiffiffiffi rk p , u Ã = 0 and the predator-prey system effectively decouples from augmentation which is accompanied by another transcritical bifurcation (T 2 in Fig 4 (top row)) between the continuing branch of stationary solutions (x Ã > 0, y Ã > 0, u Ã ! 0) and the origin. In ε > ε Ã regime, origin is the only dynamical attractor. Fig 5 shows the parameter scans for H II and H III systems highlighting these different dynamical regimes. In region A these systems exhibit periodic behavior and the boundary between A and B is the locus of the Hopf bifurcation in the ε, k plane, which leads to the stabilization of stationary solutions (x Ã > 0, y Ã > 0, u Ã < 0). B corresponds to the regime where stable stationary solutions (x Ã > 0, y Ã > 0, u Ã < 0) and (x Ã > 0, y Ã > 0, u Ã ! 0) are observed and the boundary between B and C is the locus of the second transcritical bifurcation T 2 which leads to the stabilization of the origin. Therefore by using appropriate values of ε and k, we can keep the system in either a periodic state, or a stationary state with non vanishing populations and can expect this control to work in experiments and be robust with respect to demographic noise; Ref. [26] experimentally stabilized a stationary solution in an electronic Lorenz system at permitted noise level. Furthermore, in the other instances of augmented predators, or augmented predators and preys we already observe a complete lack of origin stabilization. Therefore, we can employ these schemes as well to avoid stabilizing the origin but one needs to be careful since these cases can lead to other complications as discussed. Another useful application for these observations could be in cases where maximization of prey yield is required. Augmenting the predator populations is seen to stabilize the equilibrium where the prey populations exist at their carrying capacity and the predators vanish. This can find applications in fisheries [53,54], algae fuel generation [55,56]; where maximal sustainable yields are crucial, and also in biomedical research, for e.g. in HIV-1 infection models [57] where a portion of human immune system i.e. activated CD4 + T cells are the primary target of the HIV-1 infection [58,59] which can be modeled via predator-prey dynamics.

Additional details
This section presents certain dynamical aspects related to augmented harmonic and Duffing oscillators for completeness. These have been compiled in a separate section to preserve the manuscript flow and for those who might be interested in these details.

Limitations of Linear Augmentation
From this expression we see that the trajectories will exponentially decay to the origin 8 k m < 0 and diverge 8 k m > 0. Fig 7 shows the numerical estimation of k m along with the largest eigenvalue of the Jacobian λ 1 at the origin as a function of ε; for partially (Fig 7(a)) and fully augmented cases (Fig 7(b)). k m here was calculated as the average rate of convergence/divergence in the Euclidean distance of the current systems' state from its previous state, for every time step along the trajectory. These results suggest that k m = λ 1 and this observation has some interesting consequences. For k m = λ 1 < 0, we have a case of a harmonic oscillator being driven by an exponentially decaying force and consequently the system settles on the origin as t ! 1. Similarly the other case of an unstable origin (k m = λ 1 > 0) is equivalent to the oscillator under the influence of an exponentially diverging force (energy being pumped into the system) which leads to diverging trajectories as time increases.
Duffing system: hysteresis. For the bistable, completely and identically augmented Duffing system in Eq 8 with x Ã = 1, y Ã = 0, the largest Lyapunov exponent was calculated for increasing and decreasing values of augmentation strength ε. Starting with initial conditions leading to the solutions (x Ã − , y Ã − , z Ã − ) at ε = 0.1, the initial conditions for the next calculation at ε = 0.1 + δε were considered as the final values of x, y, u from the previous calculation for ε = 0.1, with δε = 0.001 and so on for the entire range in the forward direction. Similarly for backwards calculation, the process was repeated starting from ε = 0.7 where (1, 0) is the only stable solution with δε = −0.001. The results of the calculation are shown in Fig 8 and as one would expect, this system exhibits hysteresis in the interval of bistablity.
Duffing system: (ε, k) plane behavior. Different dynamical regimes for the Duffing system in Eq 8 with (x Ã , y Ã ) = (1, 0) and (0, 0) are shown in (again from Eq 13) and the boundary of the blue dot regime gives the locus of this saddle node bifurcation; parabolic function F SN (ε, k)(= 5ε 2 −k) = 0, estimated from the expression for x Ã ± as the limiting value of k and ε to get real solutions, which is obtained by equating the discriminant in the expression of x Ã ± to zero. The boundary between regimes A and B corresponds to the locus of the transcritical bifurcation between (1, 0) and (x Ã + , y Ã + , z Ã + ), and is given by the zero crossing of the largest eigenvalue for (1, 0) which satisfies F TC (ε, k)(= 2ε 2 −k) = 0. In region B either the now stable coupling dependent stationary solutions (x Ã + , y Ã + , z Ã + ), or escaping trajectories are observed.
Similarly for x Ã = 0, y Ã = 0 in Fig 9(b): B highlights the regime where the system settles on the solutions (x o + , y o + , u o + ) from Eq 16. The blue dots correspond to the initial conditions leading to stationary solutions (x o − , y o − , u o − ) also from Eq 16. Since the system is bistable in this regime, we should expect the entire region B to be filled with these blue dots but that is not the case. The reason behind this behavior is a difference in the relative basin size of these two solutions; number of initial conditions leading to (x o + , y o + , u o + ) is more than the ones which lead to (x o − , y o − , u o − ). This difference is even more pronounced for higher values of k. Both these solutions vanish via a pitchfork bifurcation and the locus of this bifurcation which separates regimes B and A is defined by the function F PF (ε, k)(= 2ε 2 −k) = 0 which is obtained by equating the discriminant in the expression of x o ± to zero.

Conclusions
In this work we studied the ability of linear augmentation towards stabilizing desired stationary solutions of oscillatory systems in a more general sense. Through some very simple examples discussed in this paper, it is clear that the effectiveness of this scheme is quite sensitive to the augmentation parameters, the class of oscillatory systems considered, the stationary solutions to stabilize and also on the way the systems are augmented. Therefore, although the simplicity of linear augmentation makes it a very compelling choice for applications, a careful analysis is required to test the system for potential pitfalls associated with the scheme. As highlighted by the examples, apart from failing to target the appropriate stationary solutions, linear augmentation can also lead to other complicated dynamical situations which include escaping trajectories, stabilization of unintended stationary solutions or the stabilization of stationary solutions which are not permitted under the modeling constraints; preys existing above their carrying capacities and negative predator populations in Sec. Dissipative predator-prey models for instance. Nevertheless, one can find ways to exploit the failures of the scheme in applications.
Although we can expect to see these results in experiments, an in-depth study of the control procedure in presence of noise, and also for larger systems is required. Extending on the results in the ecological context, one needs to check the control behavior in presence of multiple preys and predators, for a food chain, and also for other functional responses [60]. Furthermore, linear augmentation has been proposed as a mechanism to control bistability [34] but how it fares in controlling more general instances of multistability including extreme multistability [47,61] is still an open question and will be addressed in subsequent studies [60].