A modeling framework for biological pest control

1.   Introduction

We develop a model describing the evolution, in time and space, of a predator prey system, where predators can be "used" to eliminate prey. Whenever prey are any sort of harmful pest or parasite and predators constitute a biological counter measure to prey diffusion, the present model provides a prototype for some sort of biological pest control strategy. Within this framework, various possible strategies can be tested and compared. While the analytic framework provides the necessary existence and stability results, it is a biological or economical question to choose which criterion has to be used to prefer a strategy to another.

More precisely, prey are assumed to diffuse all over the available space, while predators move following the prey density gradient towards regions with a higher prey density. In this setting, we introduce a term in the predators' equation able to describe the artificial insertion of predators in the environment.

The mathematical modeling of pest distribution, evolution and control has a long tradition, see for instance the review ^[1]. Biological pest control, in particular, is receiving a growing attention in the recent literature, see e.g., ^[2,3,4]. We recall in particular ^[5], where the same subject is considered, but in the entirely different setting provided by ordinary differential equations. The present approach relies on a thorough use of the techniques of non–local partial differential equations and is amenable to rigorously state and tackle a variety of optimal control problems.

From the analytic point of view, we stress that in the present model the propagation speed of predators is finite. While prey's movement is described through the usual diffusion operator (i.e., the Laplacian), here predators move according to a Hyperbolic Balance Law, which notoriously display finite speed of propagation ^[6]. Moreover, as it is reasonable from the biological point of view, each predator moves according to the distribution of prey in a region all around it, not only at its exact location, leading to a so called non–local interaction and, hence, non–local equation. Thus, the dynamics of predators essentially depends also on their horizon, i.e., on the distance at which they can "feel" the presence of prey.

As it is natural to expect, and as it is confirmed by the model's equations (2.1), various functions and parameters contribute to the description of the present predator-prey dynamics. Here, we target the role of the control parameter, namely $q$ in (2.1), which represents an external time and space dependent intervention aimed at increasing the predators' population to limit the prey harmful presence.

The variety of strategies that can be adopted is vast and many questions naturally arise: are a limited number of intense global spraying more effective than less intense but periodic treatments? Is it relevant to introduce predators uniformly all over the space domain? Are there effective and, presumably, cheaper strategies where predators are introduced only in suitable carefully selected regions? Present day analytic techniques are, to our knowledge, currently not able to single out the best strategies: the use of numerical integrations is essential and unavoidable.

The next section is devoted to the description of the model and to the rigorous statement of its well posedness and stability properties. Then, in the subsequent section, various numerical integrations allow to compare the results of various pest control strategies.

2.   Model description

The pest population is described through its density $w = w (t, x)$ while its antagonist's density is $u = u (t, x)$, where as usual $t \in {\mathbb{R}}_+$ is the time variable and $x \in {\mathbb{R}}^2$ is the space coordinate. The pest population $w$ is assumed

(w.1) to diffuse isotropically, with a movement governed by the usual Laplace operator $\Delta w$;

(w.2) to naturally increase through reproduction with a net birth rate $\gamma$;

(w.3) to be predated by the $u$ population, so that the $w$ birth rate is reduced by the quantity $\delta \, u$.

On the other hand, the predators in the $u$ population

(u.1) move where the pest density is higher, their velocity $v$ being a non–local function of the prey gradient $\nabla w$;

(u.2) have a birth rate $\alpha \, w$ proportional to the prey density;

(u.3) suffer a fixed mortality rate $\beta$.

In both cases, the latter terms (w.2)–(w.3) and (u.2)–(u.3) reflect the usual Lotka–Volterra interaction, see for instance [7,§ 3.1].

We thus get the following mixed hyperbolic-parabolic system,

which extends ^[8] thanks to the control term $q$ and fits into the more general setting of ^[9], where the various coefficients are time and space dependent.

We allow the velocity vector field of the predators $v = v (w)$ to depend on $\nabla w$ through a space average. Indeed, according to (u.1), it is proportional to the average density gradient $\nabla w$ seen by the predators within their horizon, so that, in general,

The role of adimensional function $\eta$, chosen as specified above, is to provide an average of the gradient $\nabla w$. Indeed, as soon as $w$ is sufficiently regular, $\nabla (w * \eta) = (\nabla w) *\eta$. Hence, $v (w)$ in (2.2) is $\kappa$ times the average gradient of $w$, smoothly normalized. Thus, $\kappa$ is the maximal predators' speed. The constant $\ell$ is the predators' horizon, i.e., the maximal distance at which predators feel prey.

A key role is played by the quantity $q$ in the first equation: it is this term that is meant to describe an external action aiming at increasing the $u$ population, thus allowing predators to naturally reduce the pest population $w$. In general, $q$ may depend on space and time, when describing an open–loop control, or also on the $w$ population, resulting in a feedback strategy.

The case of the various coefficients $\alpha, \beta, \gamma, \delta, \kappa$ or of $v$ to be explicitly dependent also on time $t$ and on the space coordinate $x$ can also be included in the well posedness and stability results related to (2.1)–(2.2), as well as the case of speed laws more general than (2.2), see ^[9] for more details.

Remark that the above choice (2.2) of the speed law is grounded on our specific modeling setting. For instance, with the different choice $v (w) = \nabla w$, which is quite usual in the different context of chemotaxis, the notion of predators horizon would be lost, as also the fact that the propagation of predators is a priori bounded. The requirement that $\eta$ is of class $\mathbf{C}^{3}$ is of a merely technical nature.

3.   The analytic framework

The next result ensures the well posedness of the model (2.1)–(2.2) and provides the key stability estimates that describe the dependence of the solution on the initial datum $(u_o, w_o)$ and on the open-loop control $q$.

We refer to ^[8,9] for the rigorous description of the analytic framework. Here, we only recall that by solution to (2.1)–(2.2) we mean a pair $(u_*, w_*)$, where $u_*$ is a weak entropy solution to the balance law $\partial_t u + \mathinner{\mathop{{{\rm{div}}}}} \left(u \, v (w_*) \right) = (\alpha \, w_* - \beta) \, u + q$, while $w_*$ is a solution to the parabolic equation $\partial_t w - \mu \, \Delta w = (\gamma - \delta \, u_*) \, w$.

Theorem 3.1. Fix positive $\alpha, \beta, \gamma, \delta, \kappa, \mu$ and $\ell$. Let $v$ and $\eta$ be as in (2.2). Then, for any initial datum $(u_o, w_o)$ and for any control $q$, with

problem (2.1)–(2.2) admits a unique global solution $(u_*, w_*)$ with

Moreover, with obvious notation, we have the following estimates on the dependence of the solutions $(u_*, w_*)$ and $(u_*', w_*')$ on the initial data $(u_o, w_o)$ and $(u_o', w_o')$:

On the other hand, if $(u_*, w_*)$ and $(u_*', w_*')$ are the solutions corresponding to the control $q$ and $q'$,

The quantity $\mathcal{O} (t)$ depends on the parameters in (2.1)–(2.2), on $q$ and on upper bounds on norms of the initial data in ${\mathbf{L}^1}$, $\mathbf{L}^\infty$ and on the predators' initial total variation; typically, $t \to \mathcal{O} (t)$ grows exponentially.

We refer for the proof of Theorem 3.1 to ^[9], where the more general case of time and space dependent coefficients is considered and further stability estimates are provided. In particular, [8,Lemma 4.1] shows that the assumptions above on $\eta$ yield a non–local operator $v$ in (2.2) meeting the requirements necessary for Theorem 3.1 to hold.

Here, we stress the qualitative information provided by Theorem 3.1. First, it ensures the existence, uniqueness and continuous dependence on the initial data of the evolution prescribed by (2.1)–(2.2). Then, it opens the way to the search for optimal strategies, ensuring the continuous dependence of the solution on the control. Mostly, Theorem 3.1 tells that the relevant norm, that is, the one that allows to consider data or controls as near or not, is the ${\mathbf{L}^1}$ norm. In particular, we do not claim that the solution, in particular $w$, can not be more regular than merely ${\mathbf{L}^1}$.

For completeness, we also note that the dependence of the ${\mathbf{L}^1}$ difference between solutions corresponding to different controls or initial data resembles that typically appearing in ordinary differential equations. As time increases, this difference may well grow unboundedly. Concerning the dependence of the solution on $\eta$, ${\mathbf{L}^1}$ estimates can be found in ^[8,9], while the case of $\eta$ approaching a Dirac delta is treated in ^[10].

4.   Comparing strategies

4.1.   General setting

The next paragraphs show various evolutions as described by the model (2.1)–(2.2). Our aim is to investigate the different evolutions resulting from different control choices. Therefore, we let only $q$ in (2.1)–(2.2) vary, keeping all other parameters (both physical and numerical) fixed.

We use throughout the Lax–Friedrichs scheme [11,§ 12.5] with dimensional splitting [11,§ 19.5] and a further splitting to take care of the source term [11,§ 17.1]. Further details on the numerical algorithm exploited below can be found in ^[12], while related numerical studies are in ^[13,14].

Solutions are displayed on the domain $\Omega = [-1, 1] \times [-1, 1]$ and along its boundary $\partial\Omega$ we fix $u$ and $w$ to $0$, which roughly corresponds to set all parameters in (2.1) and the initial data to $0$ outside $\Omega$. The presence of the non–local speed (2.2) imposes to use a numerical domain larger than $\Omega$, that is $[-1-\ell, 1+\ell] \times [-1-\ell, 1+\ell]$, where $\ell$ quantifies the distance at which predators "feel" prey, see (2.2).

The integration time is $T = 12$, the space mesh consists of $2^9\times2^9$ squares. Moreover, with reference to the notation in (2.1)–(2.2), we set:

The initial datum used in this section is detailed in (4.2), see also Figure 1.

Remark that the diffusion term in the second equation causes the initial position of the prey to be essentially irrelevant at large times. To evaluate the efficiency of each strategy, we choose the total number of prey, averaged over time between time $t_o = 4$ and the final time $T = 12$, that is we compute the quantity

Our dropping the initial interval $t \in [0, 4]$ is aimed at reducing the initial oscillations due to the particular initial datum chosen. As it is clear in the examples below, it may well happen that different strategies lead to similar values of the average $\overline{W}$, but with very different ranges in the values attained by the total number $\int_\Omega w$ of prey at different times. Therefore, we also measure the variation of $\int_\Omega w$ in time through the ${\mathbf{L}^1}$–oscillation defined as^*

^*If $A$ is a measurable subset of ${\mathbb{R}}^m$, by ${\left|{A}\right|}$ we denotes its Lebesgue measure.

Our preferring the ${\mathbf{L}^1}$ based quantity above is justified by system (2.1)–(2.2) being well posed in ${\mathbf{L}^1}$, see ^[9]. Note however that in a realistic environment several other aspects of the pest distribution can be relevant, such as the size of their peaks or their location, or also their reaching very low values, essentially suggesting extinction.

First, as a reference case, we integrate (2.1)–(2.2) with no external input of predators, i.e., with $q \equiv 0$, which corresponds to the case considered in ^[8,12]. The total number of predators and prey display a typical Lotka–Volterra like behavior, with synchronized oscillations in the total number of individuals in the two populations, see Figure 2, left. The geometric distribution of prey is strongly affected by the very fast diffusion driven by the Laplacian, see Figure 2, right. In turn, the predators' geometric distribution follows, since predators move towards high prey density regions, due to (2.2).

The resulting key final total values of predators and prey, as well as the values of $\overline{W}$ and of $\mathcal{O}_w$ are the following:

Below, while comparing different strategies, an obvious key role is played by the quantity $\int_0^T \int_{ {\mathbb{R}}^2} q (t, x) \mathinner{\mathrm{d}{{x}}} \mathinner{\mathrm{d}{{t}}}$, which represents the total amount of predators introduced in the environment during the time interval $[0, T]$. Clearly, in general, higher values of this integral lead to higher effects on the prey population. Therefore, whenever $q$ is present in the integrations below, we keep

4.2.   Periodic treatments

The first strategy we investigate consists in releasing predators periodically into the environment all along two strips symmetric with respect to the origin:

We now let the frequency $\nu$ in (4.7) vary. The resulting average presence $\overline{W}$, as defined in (4.3), of the $w$ population is affected up to about $10\%$ by changes of $\nu$ in (4.7), while the oscillation $\mathcal{O}_w$, defined in (4.4), spans a wider range:

The best performance, measured through the value of $\overline{W}$, is obtained with $\nu = 2/3$. The total populations resulting from the integrations with the smallest and largest choices of $\nu$ are shown in Figure 3. These results open the way to a variety of further investigations. Indeed, while the mean presence of pest $\overline{W}$ is similar in the two situations (with variations of about $10\%$), the two evolutions are entirely different. A single criterion can be too little to evaluate the success of a strategy. When $\nu = 1/4$ pests almost vanish and then rise to rather high values. On the contrary, when $\nu = 2$, the total pest population displays small oscillations around the mean value. Clearly, on the basis of considerations specific to the real situation at hand, a measure of the dispersion of $t \to \int_\Omega w (t, x)\mathinner{\mathrm{d}{{x}}}$ should better be taken into account. A possible choice can be for instance $\mathcal{O}_w$, as defined in (4.4), so that the optimal strategy $q$ might be sought as the one minimizing a function of the type $q \to \overline{W} (q) + \lambda \; \mathcal{O}_w (q)$, for a carefully chosen positive weight $\lambda$.

It is now natural to investigate the effects of strategies similar to that in (4.7), but alternating the insertion of predators between the two strips, that is:

Clearly, the case $\nu = 0$ is equivalent to the reference case (4.5). The values of $\overline{W}$ resulting from different frequencies, defined in (4.3), are now:^†

^†In some of the integrations detailed here, the chosen values of $\nu$ lead to a small error in the constraint (4.6).

The values of the cost $\overline{W}$ are similar in the two cases (4.10) and (4.8), now the minimal cost being obtained with $\nu = 2/3$ in (4.7). However, the behavior of the total amount of prey is significantly different, as it clearly results comparing Figure 4, related to the present case (4.9), to Figure 3. Indeed, alternating the introduction of prey in the two strips as in (4.9) appears significantly more effective in reducing the peaks in the pest density than the strategy (4.7). Coherently, the range of values attained by the oscillation $\mathcal{O}_w$, as defined in (4.4), is narrower than that in (4.8), as shown in the table in (4.10).

An example of the dynamics resulting from (2.1)–(2.2)–(4.1) with initial datum (4.2) and control (4.9), where $\nu = 1/4$, is in Figure 5: the deployment of predators reaches its maximum for $x_2 \in [0.5, \; 0.75]$ at $t = 1$, leftmost diagrams, and for $x_2 \in [-0.75, \, -0.5]$ at $t = 3$, rightmost diagrams. At time $t = 2$, middle diagrams, the same amount of predators is released in the two regions $x_2 \in [0.5, \; 0.75]$ and $x_2 \in [-0.75, \, -0.5]$. The effects of the increased presence of predators is a temporary reduction in the pest density in the same region.

For completeness, we remark that the actual deploying of alternated treatments might reasonably have a higher cost.

4.3.   Global Intense Treatments

A strategy entirely different from the one considered above consists in deploying predators uniformly over all of $\Omega$. Below we investigate in detail the effects of treatments applied for relatively short, but carefully chosen, time intervals.

First, as a reference situation, we present the effects of a uniform insertion of predators all during the interval $[0, T]$ and all over the region $\Omega$, always fulfilling the constraint (4.6). The plots of the resulting total amount of predators and of prey are in Figure 6. Remarkably, this strategy leads to the rather low value $\overline{W} = 11.932$, but pests remain significantly present throughout all the time interval $[0, T]$ with relatively small oscillations, as testified by the value $\mathcal{O}_w = 2.185$.

We now consider the (probably to some extent rather realistic) strategy consisting in a single treatment deployed to the whole domain for a relatively short period, so that

Different choices of $\bar q$ and $I$ lead to quite different evolutions of the model (2.1)–(2.2). When effective, such a treatment causes a sharp increase in the $u$ population which, in turn, may bring the pests $w$ near to extinction. Note however that, as it is usual in Lotka–Volterra type systems, as soon as $w$ is not identically and exactly $0$, the term $\gamma \, w$ in the source of the second equation in (2.1) eventually ensures an exponential growth to the $w$ population. Therefore, the value of $\bar q$ and the position of the interval $I$ have to be carefully tuned to the interval of interest $[t_o, T]$.

We choose $I$ of length $2$, so that $\bar q = 150$ according to (4.6). The integration of (2.1)–(2.2)–(4.1)–(4.11) thus yields the following results:

Table (4.12) clearly shows that the cost $\overline{W}$ is not a monotone function of the position of the interval $I$ (or, more precisely, say of its minimum). Indeed, the fact that the cost $\overline{W}$ depends on the solution computed on the interval $[4,12]$ imposes that applying a treatment too early may cause the pest first to essentially disappear, but then to be reborn before time $T = 12$. On the contrary, a late treatment allows the growth of the pests for a too long time and is then not able any more to significantly reduce them. Indeed, the first treatment (see Figure 7, left), is carried on the time interval $[3,5]$ and brings prey almost to extinction. But it is too early and pests have the chance to re-develop. Similarly to a typical Lotka–Volterra dynamics far from the stable equilibrium, the total masses of predators and prey display very wide oscillations.

The opposite situation is that of a too late treatment, see Figure 7, right. Here, predators are able to wipe out pests, but it happens too late, so that the contribution to $\overline{W}$ from the times before predators were introduced is too high to make this strategy convenient.

The middle diagrams in Figure 7 display the most convenient choice among the ones considered. Pests essentially disappear thanks to the insertion of predators and are not able to recover before the final time $T = 12$.

The present global intense treatments introduce a peak in the predators' density and, as a consequence, pests almost disappear. Thus, the ${\mathbf{L}^1}$–oscillation (4.4) only gives a rough measure of the dispersion of the prey density around the value $\overline{W}$.

A possible refinement of the strategy (4.11) might consist in repeating the insertion of predators twice, so that

In the present setting, the only constraint on (4.13) is (4.6), which now reads

A thorough optimization of $\overline{W}$ as a function of the $6$ variables $\bar q_1, \bar q_2, \inf I_1, \sup I_1, \inf I_2, \sup I_2$ subject to the constraints (4.13) and (4.14) is now a question of a conceptually simple but numerically quite heavy procedure, the Lipschitz continuous dependence of the solution to (2.1)–(2.2)–(4.1)–(4.13) from these variables being ensured in ^[9], see also Theorem 3.1. Here, we only consider the two sample cases detailed below:

Again, it stems also from the integrations detailed in Figure 8 that the choice of the strategy $q$ may lead to a variety of evolutions. Depending on the relevance for the controller of the possible appearance of peaks in the pest population and of their geometric distribution, the fine tuning of the various parameters can lead to quite different results. For instance, a key role is played by the positions of these peaks resulting from the Lotka–Volterra type dynamics relative to the time interval chosen for the optimization.

Note that, in spite of the high peak appearing in the pest distribution after time $t = 10$, the latter integration results in the lowest values of $\overline{W}$ obtained so far: $\overline{W} = 10.44$.

4.4.   The Role of Geometry

The present modeling framework allows to investigate the role of the regions where predators are deployed. We consider below two somewhat extreme situations, with $q$ supported either in the corners of $\Omega$, so that

or supported in the center of $\Omega$, that is

Once again, the fact that pests' movements are driven by the Laplacian, with predators hunting for pests, characterizes the overall behavior of both populations. Prey accumulate at the center of $\Omega$, attracting predators. In the former case, the "hunting" attitude of the predators described by the speed law (2.2) makes them move towards the center of $\Omega$, where prey accumulate, see Figure 9, left. In the latter case, predators are released directly where the pests density is higher. Hence, treatments that privilege deployments of predators at the center of $\Omega$ can be more effective, as confirmed by the comparison of the values of $\overline{W}$ in the two cases:

In both cases (4.17) and (4.18), in the long run the pest population is stabilized and oscillations are relatively weak, see Figure 10.

We note, as it has to be expected, that the dynamic equilibrium in the case (4.17) of $q$ supported near to corners also leads to a higher average presence of predators, coherently with the higher average presence of prey.

5.   Conclusions

The above paragraphs show that by means of model (2.1)–(2.2) we are able to test a variety of pest control strategies. Above, we describe some of the possible scenarios that fall within the proposed approach. Many different criteria can be used in testing, comparing and choosing among the different strategies. Here, particular relevance is given to the cost $\overline{W}$ and to the quantity $\mathcal{O}_w$ as a measure of the variability of the solution, but a variety of alternatives are conceivable and can be addressed through the present techniques.

A key aspect in the dynamics described by (2.1)–(2.2) is the presence of the Laplace operator in the second equation to describe pest diffusion. While it is a very usual choice, it allows an extremely fast propagation of the pest population, which is in contrast with the predator speed, bounded by $\kappa$, see (2.2). Whenever modeling reasons justify it, a possible improvement of the present model might consist in modifying the second equation in (2.1) so that also prey move with bounded speed.

Acknowledgment

The first author was partly supported by the PRIN 2015 project Hyperbolic Systems of Conservation Laws and Fluid Dynamics: Analysis and Applications and by the GNAMPA 2018 project Conservation Laws: Hyperbolic Games, Vehicular Traffic and Fluid dynamics. The second author acknowledges the support of the Lorentz Center. The IBM Power Systems Academic Initiative substantially contributed to the numerical integrations.

Conflict of interest

The authors declare there is no conflict of interest.