Abstract
The paper is concerned with the development and numerical analysis of mathematical models used to describe complex biological systems in the framework of Integrated Pest Management (IPM). Established in the late 1950s, IPM is a pest management paradigm that involves the combination of different pest control methods in ways that complement one another, so as to reduce excessive use of pesticides and minimize environmental impact. Since the introduction of the IPM concept, a rich set of mathematical models has emerged, and the present work discusses the development in this area in recent years. Furthermore, a comprehensive parametric study of an IPM-based impulsive control scheme is carried out via path-following techniques. The analysis addresses practical questions, such as how to determine the parameter values of the system yielding an optimal pest control, in terms of operation costs and environmental damage. The numerical study concludes with an exploration of the dynamical features of the impulsive model, which reveals the presence of codimension-1 bifurcations of limit cycles, hysteretic effects, and period-doubling cascades, which is a precursor to the onset of chaos.
1. Introduction
Food losses due to pests and plant diseases are nowadays one of the major threats to food security, particularly in large parts of the developing world. As reported by the United Nations [1], the world’s population in 2014 was estimated at 7.2 billion, with an approximate yearly growth of 82 million, a quarter of which occurs in the least developed countries. This unprecedented amount of people in the world poses serious challenges for food producers and policy-makers, especially regarding the minimization of crop losses due to pests and plant diseases, which have been estimated to be as high as 40% of the world production [2]. This issue has been a matter of active research for many decades, where the main challenge lies in the unavoidable trade-off between pest reduction, financial costs, effects on human health, and environmental impact. Therefore, the problem of pest control has necessarily to be addressed in an integrated manner, which has motivated the development of various integrated approaches, such as Integrated Pest Management (IPM) [2, 3]. IPM’s basic principle consists in the judicious and coordinated use of multiple pest control mechanisms (e.g., biological control, cultural practices, and selected chemical methods) in ways that complement one another, maintaining pest damage below acceptable economic levels, while minimizing hazards to humans, animals, plants, and the environment. The literature on Integrated Pest Management is vast, and Section 2 will present a discussion of the historical development of the IPM concept over the past decades, which will serve as motivation for the mathematical description of the underlying pest control methods (see below).
One of the critical factors of success in the implementation of IPM programmes is the fundamental understanding of the interplay between the different elements of the associated agricultural ecosystems, such as crops, pests, natural enemies, and biopesticides, which quite often can only be achieved via a certain degree of mathematical abstraction. The mathematical description of ecological processes is one of the main subjects of study in the field of Ecological Modeling, which can be considered to a great extent as a nonlinear science, since almost all ecological interactions, both trophic and competitive, are nonlinear. In this respect, one of the main challenges in this area is the construction of mathematical models that provide reliable predictions and understanding of field observations in real ecosystems, in such a way that these models can be used as decision support tools. The possible approaches that can be employed to tackle this issue have been a matter of extensive debate among scientists for at least four decades [4, 5], and discussions on this topic are still going on [6–11]. In the particular case of agricultural pest control, the challenge consists in developing robust mathematical models able to at least qualitatively describe different pest control methods, so as to make the development of pest control strategies and policies less intuitive or empirical.
Since the origins of the IPM concept in the late 1950s [12, 13], a rich set of mathematical models has emerged in the literature focusing on different features of IPM-based applications, and in the present work we will briefly discuss several aspects of the recent development in this area (Section 3). One of the first systematic reviews of IPM-related mathematical models was published by Shoemaker [5], which was based upon her doctoral dissertation submitted in 1971. She considered suitable combinations of chemical and biological pest control applied to simplified ecosystems consisting of crops, pests, and parasites. The parasites are able to naturally control the pest population to a certain extent, but this effect is compromised by the pesticides applications, as Shoemaker assumed that the parasites are also killed by the chemical. As a practical example, she developed charts that the growers can use to determine whether pesticides should be sprayed, given in terms of the time until harvest and pest and parasite densities. Subsequent surveys of mathematical models for controlling pest populations have been carried out by Jaquette [4], Wickwire [14], Beddington et al. [15], and Barclay [16], and Section 3 will discuss some representative models that have been further introduced since then.
Despite the vast literature on modeling pest control strategies in agricultural ecosystems, comprehensive parametric studies of the mathematical models are rather scarce, with numerical investigations conducted mainly at the simulation level. This fact motivates the main part of the present work (Section 4), which will be devoted to the utilization of specialized numerical techniques in order to address practical questions relevant to IPM applications. Specifically, an impulsive pest control model will be chosen and reformulated as a hybrid dynamical system [17], thus allowing the parametric study of the periodic response of the system by means of numerical continuation (path-following) methods [18]. In this way, we will be able to tackle questions as to the optimal implementation of the impulsive pest control, in terms of minimizing operation costs and environmental damage. The paper will end with further numerical investigations of the dynamic response of the impulsive model, which will reveal the presence of codimension-1 bifurcations of limit cycles, hysteretic phenomena, and period-doubling cascades, which is a precursor to the onset of chaos.
2. A Brief Overview of Integrated Pest Management
An apparently promising approach to reduce crop losses due to pests appeared during the 1940s, with the discovery of synthetic pesticides such as DDT, which marked a new era in pest control. Pesticides became soon popular in the agricultural industry, as they were easy to apply and effectively killed a significant amount of the targeted pest, due to which their use spread rapidly worldwide. However, the overreliance on synthetic pesticides was proven to be unsustainable already by the end of the 1950s. Just a few years after the first use of DDT, resistance to the chemical was observed in a variety of insect pests [2], which meant more frequent applications and higher dosages of pesticides in order to keep acceptable pest population levels. Another negative effect was the reduction of beneficial species (such as natural enemies), which intensified the problem of pest resurgence and allowed nonpest species to increase in number and become pests themselves. In addition to these on-site crop problems produced by the overuse of pesticides, their negative impact extended beyond the agricultural framework, causing damage to water sources and further ecosystems, as well as posing serious human health hazards due to, for example, pesticide residuals in food and pesticide exposition [19, 20].
Recognition of the problems associated with the indiscriminate use of synthetic pesticides encouraged the development of alternative pest control paradigms, such as the concept of Integrated Pest Management (IPM) [2, 3, 12, 13]. Having its origins in the seminal work by Stern et al. [21], IPM is an interdisciplinary pest control approach that relies heavily upon natural mortality factors, such as natural predators and environmental conditions, combined with further control mechanisms. These include biological control, selected chemical methods, and cultural practices. The basic idea is that, instead of employing a single control method, efforts are directed to the judicious and coordinated use of multiple tactics in ways that complement one another, maintaining pest damage below acceptable levels, while minimizing hazards to humans, animals, plants, and the environment.
Integrated Pest Management (IPM) has been recognized as one of the most robust constructs in agricultural sciences to deal with the challenges related to the excessive use of pesticides [12], already outlined above (see also [22–25]). The key concept for the implementation of a pest control programme in an IPM framework is that of economic injury level. This term was introduced for the first time by Stern et al. [21] (see also [3]) and means the lowest pest population density that will cause economic damage. The latter term can be defined as the amount of injury that justifies the application of control measures, and it can vary depending on the area, season, and other economic or ecological factors.
In general, it is assumed that a number of pest control mechanisms are available, for instance, biological methods, cultural practices, natural enemies, habitat management, and synthetic pesticides. The basic decision rules rely on a predefined economic injury level and an economic threshold, which gives the pest population density above which control actions must be taken so as to prevent the pest population from reaching the economic injury level. An IPM-based pest control scheme, in its simplest form, will then require that whenever the amount of pests is less than the economic threshold only ecologically benign control measures are applied, that is, those that enhance natural control. If natural control is not capable of preventing the pest population from reaching the economic injury level, then synthetic pesticides come into play, nevertheless, in adequate combination with environmentally friendly control measures so as to minimize the amount of pesticides released into the underlying ecosystem. In practice, however, the task of developing and implementing an IPM-based pest control programme sustainable in both ecological and economic terms is by no means a trivial one.
As has been recognized in the past [2, 12, 16], the implementation of IPM-based control strategies requires a profound knowledge of the interactions between the different components of the underlying agricultural ecosystem, for instance, crops, pests, natural enemies, and habitat conditions. Given the remarkable complexity of such ecosystems, the required fundamental understanding of the involved biological interactions can quite often only be obtained via a certain degree of mathematical abstraction. Ideally, the final result should be a mathematical model tailored for a specific application, which provides reliable predictions and understanding of field observations, in such a way that the model can be used as a decision support tool for devising effective pest control schemes [26]. Nevertheless, the main challenge lies in the sheer complexity of the involved biological processes, which often hinders the search for appropriate mathematical representations of the laws governing the ecosystem. In the next section we will discuss in detail some representative mathematical models used for describing various pest control methods, including synthetic (pesticides), biological (natural enemy predation, biopesticides), and cultural (roguing, replanting), which are precisely some of the most common control mechanisms used in combination in IPM-based control programmes.
3. Mathematical Models for Pest Control
In this section we will present a short overview on the available mathematical models used to describe the ecological interactions in pest control. The discussion will be mainly guided by two criteria: model type, in a mathematical sense, and the underlying ecological phenomena in the framework of agricultural pest control. Special attention will be paid to those models related to Integrated Pest Management (IPM), where the main purpose is to minimize damage to nontarget organisms and harmful environmental effects, by combining classical chemical strategies with alternative control methods, such as biological control, host-plant resistance breeding, crop rotation, harvest management, and cultural techniques (see Section 2).
Before starting our discussion, let us give some remarks regarding notation. As is well known, biological pest control is characterized by the reduction of pest population as a result of the introduction of a natural enemy [15]. The interaction between a pest and its natural enemies can be understood in terms of the dynamics observed in prey-predator models. Therefore, most of the models presented in our discussion will have in their core a certain type of prey-predator system. In this context, state variables related to pest (prey) and natural enemies (predator) populations will be denoted by scalar, nonnegative variables and , respectively. Furthermore, all system parameters are assumed to be positive numbers, unless otherwise stated, and they will be used throughout the manuscript in a consistent manner, in such a way that, whenever possible, they will have the same meaning in different models. Finally, the prime symbol will denote time differentiation.
3.1. Pest Control as a Time-Continuous Process
Pest Diseases and Natural Enemies as Control Measures. In the context of Integrated Pest Management, one of the most representative control strategies is that of an artificial spread of an infection among a pest population combined with a different control method, such as continuous pesticide spraying or natural enemy predation. Typically, the pest population is divided into two classes: susceptible and infective. The infective population is used to spread a certain disease or virus created in a laboratory. At the beginning, a small amount of infected pest is introduced into the ecosystem with the purpose of generating an epidemic. The susceptible population becomes infected through direct contact with the infective pest, thereby causing a significant reduction of the pest population as a direct consequence of the disease, or due to a decrease in its reproductive ability. A model describing this control mechanism has been recently proposed by Jana and Kar [27], based on the classical susceptible-infective (SI) paradigm [28]:The subscripts and denote the susceptible and infective pest population, respectively, while represents the natural enemy (predator) population. From the two classes of pest, only the susceptible population is able to reproduce, according to a logistic law with intrinsic growth rate and environmental carrying capacity . The infective pest is assumed to also contribute towards the carrying capacity of the ecosystem; however, the significance of this environmental impact can differ from that of the susceptible population, which is reflected by the parameter . The disease spreads within the susceptible pest through direct contact with the infected population, whose effect can be seen in the second term of the first equation of system (1), where the coefficient represents the force of infection. A second control action in this ecosystem is given by the presence of a natural enemy (predator). The predator population feeds on the pest according to a Holling type II trophic function [29–31], where is the maximum predator’s capturing rate and stands for the susceptible pest population density at which the capturing rate is half the maximum value.
As already mentioned above, the infective population is not able to reproduce; hence its survival relies heavily upon the ability of the disease to turn the susceptible population into infected pest, which depends on the force of infection . On the other hand, the infective pest is assumed to have a natural death rate . This produces an exponential decay of the infective population in the absence of natural enemies and susceptible pest. The infective population is attacked by the natural enemy following a Holling type I functional response with capturing rate ; see the second equation of the model (1).
The third equation in (1) describes the evolution of the natural enemy population. Here, and represent the conversion factors from susceptible and infective pest, respectively, into predators. A negative effect of the infection on the predator population is accounted for by the parameter . In addition, denotes the mortality rate and measures the intensity of competition among the natural enemies for space, food, and so on. At low pest population densities, the natural enemies are able to find alternative food sources, leading to an additional growth rate . As the pest population grows, the predators make less use of the alternative food, and when the pest approaches the environmental carrying capacity, the natural enemy feeds almost on the pest only.
A systematic study of the dynamics of system (1) has been carried out by Jana and Kar [27]. They analyze in great detail the effect of the relevant parameters on the proposed ecosystem, and the theoretical results are illustrated by representative numerical simulations. As has been already explained above, the pest population is controlled by two methods: infection release and natural enemy predation, which fall into the class of biological pest control. From a practical point of view, the role of the alternative food source for the natural enemy has been shown to be crucial. If the predator growth rate due to the alternative food is greater than its mortality rate (i.e., ), a nontrivial pest-free equilibrium is feasible, where the natural enemy is the only species present in the ecosystem. Furthermore, if exceeds a certain threshold, the equilibrium is asymptotically stable, which means that small (pest) perturbations in the ecosystem will be controlled entirely by the natural enemies. Jana and Kar conclude their investigation by considering a third control measure in the system consisting in chemical methods (pesticides), and they discuss how to combine these three control techniques so as to minimize the deadly effects on the natural enemies and environmental damage, while effectively reducing the pest population.
Plant Disease Control via Cultural Methods. Another nonchemical mechanism to deal with plant diseases can be implemented via cultural practices, which fall into the category of biological control. They are carried out by means of human actions only, such as roguing (removing) infected plants, replantation of disease-free plants, crop rotation, intercropping, and strip farming [20]. The main purpose is to reduce the negative consequences of a disease to levels that are acceptable in economic terms, while causing minimal damage to the environment. However, there are certain limitations for the application of such control methods, since they typically involve high labor costs and a complete eradication of plant diseases through cultural practices is generally not possible.
A low-dimensional model considering cultural practices as the only disease control action is proposed by van den Bosch et al. [32]:with . Here, and denote the population of susceptible and infected plants, respectively, in an agricultural crop field contaminated with a certain viral disease. New plants are introduced in the ecosystem at a constant rate , from which corresponds to in vitro-germinated healthy plants and to cuttings taken from the crop (susceptible and infected). By visual inspection or other diagnostic methods, infected plants from the cuttings are discarded with probability . Furthermore, it is assumed that the infected plants are able to recover due to reversion with probability . Consequently, the introduction of new plants contribute to both the susceptible and infected populations (see, e.g., the first term in the second equation of (2), which gives the proportion of infected plants that enter the crop as a result of the introduction of cuttings). In addition to the usual death rate , the infected plants are removed from the system (roguing) at a rate . As can be seen, only cultural actions are considered as a measure of disease control. On the other hand, another parameter considered in the ecosystem is the within-plant virus titre (denoted hereafter by ). As pointed out in [32], this parameter deeply influences both the viral transmission in the crops and the corresponding disease symptoms, which in turn affect the roguing rates and recovery and detection probabilities. Therefore, van den Bosch et al. assume the parameters , , , and to be -dependent, and they suggest certain functional relations to quantify these interactions, motivated by some previous field studies.
A remarkable feature of model (2) is the fact that it illustrates two well-known disease transmission modes: horizontal and vertical. The first one typically occurs through herbivorous insects that ingest spores on host-plant leaves and carry the disease from one plant to another. The intensity of this transmission mode is quantified by the transmission coefficient . On the other hand, the vertical transmission is caused by the use of cuttings to establish new crops. In this respect, van den Bosch et al. focused special attention on the effects of a trade-off between horizontal and vertical transmission modes operating in plant pathogens. Specifically, the authors analyze in detail the underlying cultural practices and how they are reflected in the proposed ecosystem (2). Their results then suggest that the considered control mechanisms should be carefully adapted and combined in order to minimize the risks of failure, which to a great extent depends on how effectively the vertical and horizontal transmission modes are adequately dealt with. Analogous models considering disease control based on cultural practices have been proposed and investigated in the past and can be found in, for example, [33–38].
3.2. Pest Control Strategies via Impulsive Perturbations
As already explained before, pest control strategies based on an Integrated Pest Management approach involve a suitable combination of biological, cultural, and chemical control techniques, including pesticide spraying, pest harvesting or trapping, plant roguing, and acceleration of the pest mortality through the introduction of natural enemies or pest diseases. Some of these control methods were mathematically described in the models presented in Section 3.1, in terms of smooth ordinary differential equations. In doing so, the pest control strategies were conceived as a time-continuous process of autonomous type, where the system parameters have to be chosen in such a way that the underlying control mechanisms effectively reduce the negative pest effects on agricultural crop fields, without any external perturbation. However, this approach neglects the fact that pest control methods are in reality implemented in a discontinuous manner, which is a consequence of the discrete nature of human activities. Furthermore, there can be exogenous factors in the ecosystems (e.g., temperature, air composition, and further human actions) that may lead to pest population densities changing very rapidly in a short period of time, which can be modeled via impulsive perturbations.
In our case, the above-outlined sudden changes in the ecosystems will be described in the framework of impulsive differential equations. This class of models is particularly suited for the representation of dynamical phenomena subject to short-term perturbations whose duration is negligible in comparison to the duration of the system evolution. Therefore, these perturbations can be assumed to act instantaneously in the form of impulses, which generally leads to jumps in the state space (discontinuous evolution). Processes of this nature can be found in numerous applications, for instance, in mechanics, population dynamics, ecology, biology, and economy, and the theoretical foundations have been developed to a great extent [39–42].
Susceptible-Infective Control Scheme with Impulsive Roguing and Replanting. A straightforward mathematical description of a pest control strategy based on impulsive perturbations is given by Tang et al. [43]:
Here, cultural practices (roguing and replanting) are considered for pest control, and the state variables and parameters are of the same nature as those of system (2), with . The actions of planting and roguing are assumed to be carried out in a periodic and impulsive manner, with period , which differs from the approach employed in the model (2), where those actions are continuously executed. As can be seen, the key parameters in this control scheme are the roguing coefficient and the impulse period , that is, the parameters that are directly influenced by human decisions. Consequently, the main question here is how to choose those parameters so as to keep the number of infected plants below a certain predefined critical level. In this respect, Tang et al. identify a lower bound for the impulse period that guarantees the extinction of the infected plants, for the special case (i.e., replanted plants contribute only to the healthy class). Specifically, they prove the existence of a periodic solution of system (3) of the form , , which is asymptotically stable provided is large enough.
Periodic Impulsive Control at Different Fixed Times. When different control measures are combined, for example, natural enemy release and pesticide spraying, it can be beneficial to carry out those actions at different moments within a period of control. This is particularly convenient when the chemicals used to control the pest population also have the collateral effect of killing the natural enemies. Therefore, a suitable amount of natural enemies has to be regularly introduced into the ecosystem at an appropriate time, in order to compensate for the undesired pesticide-induced predator mortality. Another reason for which control measures applied at different times can be advantageous is that some natural enemies are effective only at certain life stages of the pest population. For instance, some predators are able to attack the adult population only, hence leaving larval or other previous stages of the pests unaffected. In such cases, the natural enemy releases and pesticide spraying have to be carried out according to the life cycle of the pest population at different times, so as to effectively cover all possible life stages of the pest.
Zhang and coworkers have proposed several pest control models based on the scheme described above, for example [44],where and stand for the densities of the pest and natural enemy population, respectively, and the parameters have the usual meaning as in the previous models, with and . As can be seen, two control actions applied at different times are considered in model (4). The first one consists in spraying pesticides into the ecosystem, which has not only the effect of killing pest individuals but also of reducing the predator population, by proportions and , respectively. In order to compensate for the undesired effect of the applied pesticide, a fixed amount of natural enemies is periodically introduced into the ecosystem, which corresponds to the second control action in system (4).
The trophic interaction between the pest and its natural enemy is characterized by a Beddington–DeAngelis functional [31, 45, 46] (see the second term in the first equation of model (4)), and it has some similarities with the Holling type II trophic function described before. The coefficient is a weighting factor that determines how fast the predator’s capturing rate approaches its saturation value as the pest population increases. In addition, the Beddington–DeAngelis functional considers mutual interference between the natural enemies, and the intensity of this effect is regulated by the parameter . Furthermore, the constant gives a measure of the abundance of pests and natural enemies relative to the ecosystem in which they interact, and it can be interpreted as a protection provided to the pest by the environment.
In the study presented in [44], the authors established a critical value of the impulse period that separates the system behavior into two cases. In the first one, model (4) admits a stable periodic solution of the form , ; that is, the pest population can be completely eradicated by means of the control measures. To prove this, the authors give an explicit construction of the fundamental solution matrix of the linear variational equation around the pest-free periodic orbit, along with the corresponding transition (also referred to as saltation [47]) matrix. Based on the resulting explicit forms, an integral condition is derived so as to guarantee that the Floquet multipliers of the periodic orbit lie within the unit circle, thus ensuring the local stability of the pest-free solution. If the so obtained integral stability condition is not satisfied, then a second type of system response takes place in which the pest-free trajectory is not stable anymore. Instead, the system possesses a periodic solution corresponding to the case when the pest and its natural enemies coexist, and the stability of this system response is shown to strongly depend on the value of the impulse period and the proportion of pest population killed by the pesticide .
Apart from the preliminary results mentioned above, there is little information on how the remaining system parameters actually affect the behavior of model (4), although various studies related to this class of pest control are available [48–54]. For instance, a crucial feature of this model is that the control actions of natural enemy release and pesticide spraying are carried out at different times, and this effect is controlled by the parameter , whose influence on the underlying pest control strategy has not been systematically discussed in the past. Therefore, in Section 4 we will use system (4) as a toy model to show how specialized numerical methods (based on path-following algorithms for nonsmooth dynamical systems) can be used to study in detail the behavior of such models under parameter variations.
4. Numerical Analysis of a Pest Control Scheme with Impulsive Effects
In the previous section, various existing mathematical models were presented and discussed in detail, with special attention given to those describing agricultural pest control methods in the framework of Integrated Pest Management. The model types ranged from smooth ordinary differential equations to differential equations with impulsive perturbations. As was pointed out, several authors have contributed to the theoretical analysis of those systems, with particular emphasis on the existence and stability of pest-free solutions, as well as finding explicit thresholds for control parameters at which stability is lost (bifurcations). While this is generally a straightforward task for models belonging to the class of smooth ordinary differential equations, the situation can be more involved for other types of systems, for example, impulsive differential equations; see [43, 44, 48, 55, 56]. Although impulsive systems describing pest control methods have received a good deal of attention in the past, numerical investigations of such models are rather scarce in the literature and are mostly carried out at the simulation level. This is the main motivation of the present section, in which we will present a comprehensive numerical analysis of one of the impulsive models introduced before, namely, system (4), with particular emphasis on how specialized numerical techniques can be employed to study the model behavior and gain insight into the underlying pest control methods.
In order to investigate the dynamics of the impulsive system (4), we will employ two different kinds of numerical approaches, namely, direct numerical integration and path-following methods. As is well known [57–59], impulsive models of the type of (4) can be formulated in the framework of hybrid dynamical systems, which are characterized by a continuous-time behavior interrupted by discrete-time events [17]. In our case, these interruptions are defined by the impulse times at which the pest control actions are carried out. In order to get reliable numerical simulations of the model behavior, the impulse times need to be accurately detected, which can be achieved by means of the standard MATLAB ODE solvers together with their built-in event location routines [60, 61], as suggested in [62]. In this way, direct numerical integration will be implemented in the present work.
As will be seen later, our investigation will primarily focus on the periodic behavior of system (4), with special attention given to the pest population and its response to the control actions. Since the impulsive model is parameter-dependent, a family of periodic solutions can generically be tracked by varying one (control) parameter, which can be numerically realized via path-following (continuation) methods. These are well-established techniques in applied mathematics [18] that enable a systematic study of system solutions subject to parameter variations, without having to make recourse to direct numerical integration, which sometimes can be time-demanding and inefficient. For the analysis of periodic solutions of hybrid dynamical systems, various specialized computational tools are available, such as COCO [63], SlideCont [64], and TC-HAT [17], and the latter will be employed in the current work for the numerical study of the impulsive model. In the next section we will formulate system (4) as a hybrid dynamical system, which will allow the implementation of the model in TC-HAT.
4.1. Impulsive Pest Control Model as a Hybrid Dynamical System
Before starting the numerical investigation of model (4), it is convenient to introduce a rescaling of the system as follows: According to these transformations, we obtain the following scaled version of the impulsive system (4):where the tildes have been dropped for the sake of simplicity.
As mentioned earlier, this impulsive model can be formulated as a hybrid dynamical system. For this purpose, the trajectories are divided into smooth segments consisting of the following components: a smooth vector field that governs the system behavior during the segment; a smooth event function whose zeroes define the terminal point of the segment; and a smooth jump function, which maps the terminal point of the current segment to the initial point of the next one. Each segment is labeled with an index , , so that any solution of the hybrid dynamical system is fully characterized by its solution signature , where defines the length of the signature. This mathematical framework enables the application of path-following algorithms by means of the software package TC-HAT [17], a driver of AUTO 97 [65] for numerical continuation and bifurcation detection of periodic solutions of hybrid dynamical systems. Recent applications of TC-HAT can be found in [66–70], where the continuation package is employed to study the bifurcation scenario of various engineering applications.
Denote by and the parameters and state variables of the system, respectively, with being the set of nonnegative numbers. The auxiliary variable will be used to embed the time into the state space, in such a way that each impulsive period , , will be mapped to the interval . In this setting, a solution of the impulsive model (6) will be divided into smooth segments, as defined as follows.
Pesticide Spraying (, P-Spr). This segment occurs for , and the dynamics of the system during this regime is governed by (cf. (6))This segment terminates when a crossing with the discontinuity boundary is detected. According to the pest control scheme, pesticide is sprayed at this terminal point, and this is implemented in the hybrid dynamical system via the jump function which gives the initial point for the next segment.
Predator Release (, Pr-Re). In this segment we have that , and the system behavior is determined by the ODE (7). The segment ends when the solution hits the discontinuity boundary with the initial point for the next segment given by the jump function The second component of this function represents the control action of introducing natural enemies into the ecosystem (cf. (6)), while the third one has the purpose of resetting the variable , so that it is always kept within the interval .
Moreover, throughout our numerical investigations the following solution measures will be used:where is assumed to be a -periodic solution of the impulsive model (6). The quantities and give the maximum amount of pest and natural enemy populations attained in a period, respectively, and can be used to investigate the impact of the system parameters on the pest control scheme from a practical point of view. For instance, if we consider the auxiliary boundary conditionit is possible to trace a curve in a two-parameter space (see Section 4.2.2) for which the pest population achieves a maximum, fixed critical value corresponding to, for example, a predefined economic threshold (see Section 2).
According to the mathematical framework presented above, the pest control model (6) can be written in compact form as follows:In Figure 1 we present a periodic orbit of this system illustrating the solution segmentation introduced above. With this mathematical framework we are now ready to carry out specialized numerical investigations of this type of periodic response, via the numerical package TC-HAT.
4.2. Numerical Results
This section will be devoted to a detailed numerical study of the impulsive pest control scheme introduced in Section 3.2 and modeled by system (14). The focus will be on the effect of the system parameters on (see (12)), which can be used to monitor the amount of pest population in the ecosystem. One of the main questions here will be how the control parameters should be chosen so as to keep the pest population below certain admissible levels. In addition, various dynamical phenomena will be investigated, such as fold and period-doubling bifurcations, as well as chaotic responses. Unless otherwise indicated, the parameter values used in the numerical results reported here are given in Table 1.
4.2.1. Behavior of the Pest Control Method under One-Parameter Perturbations
As was already mentioned in Section 3.2, Zhang et al. [44] studied the stability of pest-free (also referred to as pest-eradication) periodic solutions of system (14), as the impulse period is varied. They determined a threshold , depending on some of the remaining system parameters, so that for the pest-free solution is asymptotically stable, while for the system response is dominated by periodic solutions for which the pest and its natural enemies coexist. Specifically, Zhang et al. showed that the pest-free solution undergoes a change of stability (bifurcation) at ; however, they did not determine the actual type of bifurcation that produces this qualitative change, and this is precisely the first question that will be addressed numerically in this section.
Let us then begin our study with the numerical continuation of a pest-free solution with respect to , as shown in Figure 2(a). In this picture, the solid black line denotes stable pest-free solutions as displayed in Figure 2(d). If the impulse period is increased, a critical value is found, at which the pest-free response loses stability. As was confirmed numerically, this bifurcation corresponds to a branching point [71], wherefrom two branches of periodic solutions emanate (black dashed and solid green lines). The dashed curve represents unstable pest-free trajectories, while the green one corresponds to stable periodic solutions with pests and predators coexisting in the ecosystem. The green branch has a critical point where the curve loses smoothness. This singularity is produced by a change in the position of the peak value of the pest population. For , the maximum amount of pest is attained exactly at the end of the segment (Figure 2(b)), while for the maximum value occurs at the end of the segment (Figure 2(c)). Another feature of this curve is that the amount of pest population, measured by , increases as grows. From a practical point of view, this means that, in order to keep low levels of pest population in the ecosystem, the impulse period should be chosen as small as possible, ideally below . Nevertheless, having small impulse periods may amount to high operation costs, since the control actions are carried out more frequently if is reduced, and therefore a compromise should be made.
The next step in our numerical investigation is to study the model response when further system parameters are varied, one at a time. The result can be seen in Figure 3. Figures 3(a) and 3(b) correspond to the numerical continuation of the periodic response of the impulsive model (14) with respect to and , respectively. In both diagrams, it can be observed that the predator population is not significantly affected by those parameters. On the other hand, the pest population shows a decreasing tendency when is increased, whereas the effect of on the pest population is exactly the opposite. This is consistent with the biological meaning of those parameters: an increment in means that the natural enemy enhances its ability to catch pest individuals, while a larger implies more competition between predators, which has a detrimental effect on their capturing rate. Although the numerical results indicate that both and can be effectively used to reduce the pest population, this may be in practice difficult or too expensive to implement, as it would require a certain mechanism to modify biological attributes of the natural enemy (via, e.g., genetic engineering or selective breeding).
From a practical perspective, variations in the control parameters , , , and may be more accessible for the users. As can be seen in Figure 3(d), the influence of (proportion of natural enemies killed by pesticides) on the system response is rather marginal. Even if one compares the extreme cases and , no significant difference can be observed. This is due to the fact that the pest control scheme introduces periodically a certain amount of natural enemies into the ecosystem, which compensates for the mortality of the predators due to the pesticide. On the other hand, the proportion of pest individuals killed by pesticides affects significantly the presence of pests in the ecosystem; see Figure 3(c). As this parameter approaches the upper boundary , the pest population suffers a significant reduction, which suggests that the pesticide mortality should be maximized in order to eradicate the pest. However, this can have well-known negative consequences for the environment; hence the effectiveness of the pest control method cannot be based upon such a strategy. Alternatively, one can try to find an environmentally acceptable, yet optimal, pest mortality . For instance, the local minimum shown in Figure 3(c).
The next parameter to be discussed is the number of predators introduced periodically into the ecosystem , whose impact on the model response is presented in Figure 3(e). The result is biologically consistent in that the larger , the larger the amount of predators in the ecosystem and the smaller the size of the pest population. This provides us with another effective way to control the pest; however, an increment in may imply higher operation costs or negative consequences for other species present in the ecosystem, and therefore this parameter must also be chosen with certain precaution.
From the system parameters considered in Figure 3, is probably the one that can be varied without causing major ecological damage. This parameter determines the instant at which pesticide is sprayed in each impulsive period, and its effect on the pest and predator populations is shown in Figure 3(f). The picture reveals an optimal point where the peak size of pest population in the ecosystem achieves a minimum value . Practically speaking, this indicates that the performance of the pest control method can be improved by just changing the time at which pesticide is sprayed, keeping all remaining parameters fixed. This represents another avenue that can be explored in order to apply the impulsive control scheme in the most effective way, with acceptable levels of environmental impact.
4.2.2. Optimal Implementation of the Impulsive Control Method
In the previous section, we presented a systematic parametric study of the pest control scheme modeled by system (14). Special emphasis was put on trying to find optimal operation points when one parameter is changed at a time. This process can be further refined by defining objective functions and constraints and letting two or more parameters vary simultaneously. Let us assume that, for a certain application, we have to minimize the utilization of pesticides and introduction of natural enemies in the framework of the pest control scheme studied in the preceding section. This can be motivated by the associated operation costs or environmental damage, as discussed earlier. Nevertheless, a reduction in the usage of pesticides and natural enemies evidently increases the risk of high levels of pest population in the ecosystem. Hence, we need to introduce a restriction for the optimization problem, which can be defined in terms of the size of the pest population not exceeding a predefined economic threshold . Mathematically speaking, we will consider the following optimization problem:Here, the objective function is nothing but the Euclidean norm of the vector , which gives a measure of the amount of pesticides and natural enemies used in the control scheme (control effort). Depending on the specific application, this functional can be further refined by, for example, introducing weighting coefficients, but for the sake of clarity we will carry out the numerical implementation with the simple objective function defined above.
The optimization problem (15) will be solved numerically via path-following techniques, as shown in Figure 4. Figure 4(a) presents a curve in the - plane corresponding to the combination of pesticide and introduced predators yielding a constant peak pest population . A trajectory along this curve is presented in Figure 4(c), where it can be verified that the pest population indeed does not exceed the predefined economic threshold . Figures 4(e) and 4(f) display solutions corresponding to parameter values above and below the point considered in Figure 4(c), respectively. These three panels demonstrate how moves away from the imposed economic threshold when the operation point is perturbed around the computed curve.
As was confirmed numerically, the curve shown in Figure 4(a) divides the - plane locally into two parts. The one to the left corresponds to parameter values for which the peak pest population exceeds the economic threshold , while in the region to the right we have that . Therefore, the optimal operation point must lie on the right part of the - plane, including the computed boundary curve. This point can be located numerically by monitoring the values of the objective function on the boundary and is found to be ; see Figure 4(b). The system response corresponding to this operation point is displayed in Figure 4(d). In comparison to the solution shown in Figure 4(c), it can be observed that the reduction of the pest population due to pesticide spraying is significantly smaller, meaning that less pesticide is being used in the optimal case. This is nonetheless compensated with a slight increment of the amount of predators introduced periodically into the ecosystem (from to ). In both cases the condition is satisfied; however, the environmental damage, measured in terms of the objective function defined above, is minimized for the optimal parameter values found.
4.2.3. Further Dynamical Analysis of the Pest Control Scheme
So far the main tool in our numerical study has been path-following algorithms, which enabled a detailed parametric study of the periodic response of the pest control scheme described by the impulsive model (14). In this way, we were able to tackle practical questions such as how to find the most suitable operation conditions in terms of pest reduction and minimization of harmful environmental effects. Nevertheless, the dynamic behavior of impulsive systems is a subject of scientific interest in its own right, and the present section will be devoted to this matter. Specifically, we will employ both path-following methods and direct numerical integration in order to gain a deeper understanding of the dynamics of the impulsive pest control scheme considered in our investigation.
The theoretical foundations for a numerical study of the impulsive system (14) have been established by Zhang et al. [44]. They addressed the question of existence and uniqueness of nontrivial periodic solutions in terms of a fixed point problem of a suitably defined operator. Following this approach, they also determined a threshold for the impulse period after which pest-free solutions lose stability. Moreover, in the conclusion part of [44] the authors raised the question whether chaotic behavior may be present in the system, and this will be precisely one of the main motivations for the numerical study presented in this section. Unless otherwise indicated, the parameter values used in the discussion below are those given in Table 1.
We begin our numerical study with the continuation of the periodic solution shown in Figure 1 with respect to the predator’s mortality rate , see Figure 5(a). As the parameter is varied from larger to lower values, a first qualitative change is observed at (PD1). Here, one real Floquet multiplier of the periodic orbit crosses the complex unit circle from the inside, passing through . This phenomenon is referred to as a period-doubling (flip) bifurcation of limit cycles [71] and is characterized, in the supercritical case, by the birth of a stable periodic solution with twice the period of the original limit cycle, which in turn loses stability (schematically represented by the dashed line in Figure 5(a)). This unstable solution regains stability via another flip bifurcation at (PD2), where the critical Floquet multiplier comes back inside the unit circle and the stable periodic solution with double period disappears. If is further decreased, a turning point (also known as fold bifurcation) is found at (F1), in which case a pair of stable and unstable periodic orbits collide and then disappear for lower parameter values. From this point a branch (dashed segment) of unstable solutions is born, which finishes at F2 (), where the system undergoes another fold bifurcation of limit cycles, and hence stability is regained. The last stability change is found at (PD3), corresponding to a supercritical flip bifurcation, from which a branch of unstable periodic solutions emanates. Another feature of the bifurcation diagram shown in Figure 5(a) is the presence of a parameter window for which the impulsive system possesses two stable periodic solutions that coexist (coexisting attractors [72]). This is produced by the interplay between the fold bifurcations F1 and F2, typically giving rise to hysteretic effects, as will be discussed below.
In order to gain further insight into the long-time dynamics of the pest control model, we will carry out a parametric study of the impulsive system (14) via direct numerical integration, when the predator’s mortality rate is varied. For this purpose, we fix a starting value for and then integrate the system over 300 impulse periods to allow for the decay of transients. After this, we extend the numerical integration for another interval of 100 periods and store samples of the extended solution at the times , , with being a fixed shift coefficient. This parameter is introduced so as to avoid sampling the solution at the impulse times and , . Next, is increased (or decreased, depending on the sweep direction) by a small amount and then the same procedure is repeated, now using the last sample of the previous step as initial value. This process terminates when a predefined final value for is reached.
The result of the numerical procedure described above is shown in Figure 5(b). The picture confirms the qualitative changes (bifurcations) predicted in Figure 5(a), labeled PD and F. In particular, the bifurcation diagram shown in Figure 5(b) allows us to visualize the creation and disappearance of orbits with double period, observed at the flip bifurcations and detected before. The increasing (red) and decreasing (blue) parameter sweeps reveal the presence of parameter hysteresis in the system, produced by the coexistence of periodic attractors, as predicted above. The blow-up of the red part of the bifurcation diagram, shown in Figure 5(c), contains the bifurcation points PD3 and F2 found in Figure 5(a). These two points define a parameter window in which a stable orbit of period survives. When decreases through PD3, the period- solution becomes unstable and a stable periodic orbit of twice the period appears. If is further reduced, the period- trajectory loses stability via another flip bifurcation (), giving rise to a stable periodic solution of period . This phenomenon repeats again and again as the parameter continues to decrease, until a critical value is reached where the flip bifurcations accumulate, leading to chaotic behavior. This infinite sequence of period-doubling bifurcations caused by the variation of a parameter over a finite interval is referred to as a period-doubling cascade and is one of the classical routes to chaos in dynamical systems [72]. Figure 5(d) shows the intersection of a chaotic attractor of the impulsive system (14) with a Poincaré section, for . This numerical study gives a positive answer to one of the open questions outlined by Zhang et al. [44] related to chaotic behavior, since we have identified a parameter set and a mechanism through which chaos can appear in the pest control model.
5. Concluding Remarks
The intrinsic premise in the present work is that pest control is a dynamical process. As such, mathematical models are essential for understanding and providing useful abstractions of the underlying biological phenomena and ecological interactions taking place in pest control applications. Since the introduction of the notion of Integrated Pest Management (IPM) in the late 1950s [20, 21], a rich set of mathematical models has emerged focusing on various aspects of IPM-based applications, and Section 3 was devoted to providing the reader with an overview of the development in this area. The discussion presented here was guided by two criteria: model class (in mathematical sense) and type of pest control method. As a result, the models considered in this section ranged from classical smooth differential equations to differential equations with impulses. Moreover, the types of pest control methods considered in our discussion can be briefly grouped in the following categories: chemical (pesticides), biological (natural enemy predation, biopesticides), and cultural (roguing, replanting), which are suitably combined in the framework of IPM in ways that complement one another.
Although the literature on modeling pest control strategies in agricultural ecosystems is vast, comprehensive parametric studies of the underlying mathematical models have received rather little attention in the past, with numerical investigations carried out primarily at the simulation level. This fact motivated the main contribution of the present work (Section 4), which concerned the application of specialized numerical techniques in order to address practical questions relevant to IPM. For this purpose, an impulsive pest control model was chosen (namely system (4)) and reformulated in the framework of hybrid dynamical systems, thus enabling the employment of path-following (continuation) methods for the numerical study of the impulsive system under parameter variations via the software package TC-HAT [17]. With the introduction of appropriate numerical indicators, a comprehensive parametric study was carried out in Sections 4.2.1 and 4.2.2, where the main question was how to determine the parameter values of the system so as to control the pest population in an optimal way, in terms of minimizing operation costs and environmental damage. Further numerical investigations of the model (4) were conducted in Section 4.2.3, with special attention focused on the dynamical features of the system. This study revealed the presence of fold and flip bifurcations of limit cycles, period-doubling cascades leading to chaotic behavior, and hysteretic effects.
As has been pointed out in the past [2, 12, 16], a decisive factor of success in the application of IPM programmes is the fundamental understanding of the interactions between the different components of the underlying agricultural ecosystem, such as crops, pests, natural enemies, and biopesticides. Much of the challenge lies in the sheer complexity of the involved biological processes, which often hinders the search for appropriate mathematical representations of the laws governing the ecosystem. The resulting models are usually the product of a trade-off between model simplicity, predictive capability and accuracy, biological consistency, and amenability for experimental validation and calibration. It is then reasonable to assume that the validity of a model correlates with how well it satisfies the aforementioned criteria, which makes model development by no means a trivial task. Once a valid model has been obtained, it is crucial to tackle formal questions regarding the well-posedness of the model in a mathematical sense, for instance, those related to existence and uniqueness of solutions, degree of smoothness, and dependence on initial values, which enable a confident application of numerical methods in order to explore the model behavior. Future progress in this area will therefore require a multidisciplinary collaborative work between agronomists, ecological modelers, mathematical analysts, and numerical specialists, aimed at constructing mathematical models that provide reliable predictions and understanding of field observations in real ecosystems, in such a way that these models can be used as support tools for the development of new pest control strategies and paradigms.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The first author wishes to express his gratitude to the Alexander von Humboldt Foundation (Georg Forster Research Fellowship Programme) for the financial support to carry out the present research.