Optimization in microbial pest control: An integrated approach

Applied Mathematical Modelling 34 (2010) 1382–1395

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

Optimization in microbial pest control: An integrated approach S. Ghosh *, D.K. Bhattacharya Department of Mathematics and Statistics, York University, 4700 Keele Street, Toronto, ON, Canada M3J 1P3 School of Bioscience and Engineering, Jadavpur University, Kolkata, India

a r t i c l e

i n f o

Article history: Received 19 April 2009 Received in revised form 8 August 2009 Accepted 20 August 2009 Available online 31 August 2009 Keywords: Insect pest control Viral infection Optimization

a b s t r a c t The paper deals with optimal management of agricultural pest population under integrated control arising out of viral infection and spraying of pesticide. The costs of the control measures and the profits or projected profits of the biomass of species give rise to a control theoretic optimization problem. We take a four dimensional mathematical model of pest control under viral infection and pesticide, and apply Pontryagin’s maximum principle (PMP) to find out the necessary conditions on economic as well as on ecological parameters to make the control process maximum profitable. Ó 2009 Elsevier Inc. All rights reserved.

1. Introduction In recent years several entomologists have been saying that effective pest control could be obtained using far less insecticide than is currently in use. In terms of an integrated program, on one hand we have chemical controls which are less expensive to use but with high environmental loss; on the other hand, there are biological controls which are more expensive to use, but with little environmental loss. If we are to include the environmental loss with the economic costs of controls, it may be that some combination of chemical and biological controls would give better result. Therefore an optimal pest management program must address itself to the following question: given an agricultural crop infested with pests, what is the best way to distribute controls currently available so that the cost associated with crop damage, controls, and environmental damage is minimized or the net profit under a suitable control process is maximized? This question may be answered via optimal control theory provided that the dynamics of the system involving the pests and the virus can be described analytically. Biological control, by definition can cover a broad spectrum of approaches ranging from the use of obligate parasites and pathogens, to facultative parasites and pathogens, to toxin-producing pathogens, to predators [1]. The agents of pest control which are in greatest usage at the present time are bacteria, fungi and lastly viruses [2]. Insect baculoviruses are particularly attractive as bio-pesticides because of certain few factors: they are safe for non target species and they are generally highly pathogenic, host death being the most likely outcome of an infection [3,4]. Virus production at present involves in vivo methods [5], as insect larvae are efficient producers of baculoviruses, and it is possible that easily cultured insect host could be employed to produce several different viruses [6,7]. In order to be competitive with traditional pest control methods, baculovirus production must be cost effective. Several factors influences virus production, including host insect life-stage, sex, rearing environment and nutritional quality of the diet [8–11]. In general any factor that affects the growth of the host after virus inoculation will influence the virus yield. Optimal production of the viruses is desired for greatest yield of biologically

* Corresponding author. Address: Department of Mathematics and Statistics, York University, 4700 Keele Street, Toronto, ON, Canada M3J 1P3. Tel.: +1 519 265 2841. E-mail address: [email protected] (S. Ghosh). 0307-904X/$ - see front matter Ó 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2009.08.026

S. Ghosh, D.K. Bhattacharya / Applied Mathematical Modelling 34 (2010) 1382–1395

1383

active viruses, and this requires selection of larval age for infection and standardization of inoculum in the laboratory [12– 14]. Besides, it should also confirm to the requisite quality control standards with minimum to negligible contamination of in vivo system with bacteria and other microbes [15]. Other shortcomings include sensitivity to ultraviolet light and speed of action when used as foliar applications, and lack of contact activity against the pest species in soil use. Problems with ultraviolet light can be overcome through optimization formulation, but the other problems remain to be tackled [16]. However, cost effectiveness of the control process does not only depend on the virus production in vivo, but it also depends simultaneously, on specific ecological parameters like, determining the optimum fertilization levels, optimum pesticide application rates (optimal control problems), or optimum parameter sets of given model (parameter optimization problems) for maximum profit out of the control procedure. The necessary parts for such optimization problems are a model of the given system and an objective function of the underlined profit to maximize. Though there are works on optimization [17,18] with other control measures, like sterile insect technique (SIT) [19], additional predators [20], but number of studies involving optimization of viral infection in pest control is less. In the present study, an Integrated Pest Management (IPM) model is designed to optimize the simultaneous use of chemical as well as biological measures to control the insect pest dynamics regulating crop damage. In this paper we aim to derive optimal rate of insecticide and viral pesticide application under maximum profit in the procedure. Though it apparently seems that restricted release of viral pesticide and minimum use of chemical insecticide would enhance crop production in the system, but still question remains on economic viability of the control procedure. The objective of this paper is to describe an economic optimization model incorporating a detailed analysis of the system with a pesticide control scheme, keeping all the complexities intact under dynamic conditions. Using Pontryagin Maximum Principle (PMP), we try to find out the optimal biomass of the species under optimal use of viral pesticide and chemical insecticide. The paper is organized as follows: Section 2 describes the model formulation and preliminary analyses, Section 3 describes the optimization of the model system and lastly in Section 4, we end up with discussion of the result with numerical simulation of the model.

2. Model and preliminaries Four populations are considered in the model, namely the susceptible pest population SðtÞ, infected pest population IðtÞ, predators or natural enemy PðtÞ and virus population VðtÞ as free-living stage, where NðtÞ ¼ SðtÞ þ IðtÞ denotes total pest population. Detailed model formulation and necessary assumptions were well stated in the paper [21]. To avoid some unnecessary complications in the present model, and to make the model more tractable analytically, we do not consider the predation function as Holling type II. We further assume that the system is under integrated control, i.e., a compatible mixture of chemical insecticide and viral pesticide are used for control of insect pest. A review by Jacques and Morris [22] showed that out of 10 insecticide-virus combinations, 9 resulted in an additive effect on insect mortality. We assume in this case that the chemical insecticide affects susceptible and infected pests and also it affects the natural enemies in the system, as an environmental damage due to chemical toxicity. But, it has no effect on the viruses as such. Hence, the modified form of the model in paper [21] is as follows:

  dS SþI ¼ rS 1   kSV  qS m1 uS; dt K dI ¼ kSV  nI  lIP  qI m2 uI; dt dP ¼ PðdP  P P þ C P lIÞ  qP m3 uP; dt dV ¼ lV V þ jnI; dt

ð1Þ ð2Þ ð3Þ ð4Þ

where ‘u’ denotes the rate of spraying of pesticide and j; ðj 2 Rþ Þ is called ‘virus replication parameter’. qi ’s are the damaging P coefficients for susceptible, infected pest and predator and mi ði ¼ 1; 2; 3Þ, 3i¼1 mi ¼ 1 are the fractions of the total amount of pesticide u ¼ uðtÞ used on species S, I and P, respectively. The dimensionless form of the above model is given by

ds ¼ asf1  ðs þ iÞg  v s  q1 m1 us; dt di ¼ v s  gi  ip  q2 m2 ui; dt dp ¼ pðd  p þ c0 iÞ  q3 m3 up; dt dv ¼ lv þ jgi; dt where q1 ¼ qS =kK, q2 ¼ qI =kK and q3 ¼ qP =kK, q1 ; q2 ; q3 > 0.

ð5Þ ð6Þ ð7Þ ð8Þ

1384

S. Ghosh, D.K. Bhattacharya / Applied Mathematical Modelling 34 (2010) 1382–1395

In the current paper we have stated only the modified results of the stability analysis of the underlined model as it is discussed earlier in detail [21,23]. These can be easily obtained as the chemical control considered in the Eqs. (1)–(3) is linear in fashion. 2.1. Equilibria The model (5)–(8) has the following equilibria as earlier: (i) vanishing equilibrium, E0 ¼ ð0; 0; 0; 0Þ, (ii) disease-free equilibrium, E1 ¼ ð1  q1 um1 =a; 0; 0; 0Þ,  Þ, where (iii) predator-free equilibrium, E ¼ ðs; i; 0; v

l ðg þ q2 um2 Þ; jg     l lðg þ q2 um2 Þ i ¼ a 1  q1 um1 ; al þ jg jg j g  v ¼ i; l s ¼

ð9Þ ð10Þ ð11Þ



(iv) the interior equilibrium E ¼ ðs ; i ; p ; v  Þ, where

   q um1 jg  s ¼ 1  1  1þ i; a al

 1 þ jlg ððd þ q3 um3 Þ  ðg þ q2 um2 ÞÞ  1  q1 um a   0 ; i ¼  1 þ jalg þ cjgl p ¼

1



ðc0 i  d  q3 um3 Þ;

ð12Þ ð13Þ ð14Þ

 jg  v ¼ i : l

ð15Þ

For existence of boundary and predator-free steady state, we assume the following conditions:

u<

a ; q1 m1

j>

lð g þ q2 um2 Þ : 1 g 1  q1 um a

ð16Þ

Moreover, positivity of s and p implies that

j < lQ and i >

d þ q3 um3 ; c0

ð17Þ

where

Q¼



  að1  i Þ  q1 um1 : gi 

Also from positivity of i another threshold based on

j>

j may be given as

l½ðg þ q2 um Þ  ðd þ q3 um3 Þ

2  : 1 g 1  q1 um a

ð18Þ

Observation 1. From Eq. (9), it may be noted that whenever j is increasing, the endemic value of the susceptible pest s is decreasing. Ecologically this means that through lysis in the infected pest, when more viruses are being released into the environment, it is causing a high rate of infection in pest and thus resulting in a decrease in the susceptible population. ðgþq2 um2 Þ    However, j ! glð1q um =aÞ :¼ j0 þ, implies s ! ð1  q1 um1 =aÞ, i ! 0þ, and v ! 0þ, in other words, equilibrium E tends to E1 , 1

1

and E collapses to E1 , when

j ¼ j0 . Thus j0 is the threshold value for force of infection.

Observation 2. From Eq. (16)–(18) we obtain the expressions for j0 , jmax and jmin where the interior equilibrium would remain feasible. Note that for feasibility of j0 , jmin or jmax , we define the upper limit umax of rate of spraying of pesticide ‘u’, which is given by (16).

S. Ghosh, D.K. Bhattacharya / Applied Mathematical Modelling 34 (2010) 1382–1395

1385

Thus we have the following proposition, Proposition 1. If u < umax , then for j 2 ð1; j0 , the system has only two equilibria, namely the vanishing equilibrium E0 and disease-free equilibrium E1 . Whenever j > j0 , one more equilibrium E exists. As j ! j0 þ, E ! E1 . Moreover, the positive interior equilibrium E becomes feasible, if j remains in ðjmin ; jmax Þ and also u < umax . 2.2. Stability analysis It is clear that the vanishing equilibrium E0 is locally asymptotically stable if u > a=q1 m1 . Whenever, u < a=q1 m1 , E0 becomes unstable and the axial or disease-free equilibrium E1 becomes feasible. The Variational matrix at E1 is given by

0 B B JE1 ¼ B @

1 ð1  q1 um Þ a

ða  q1 um1 Þ

ða  q1 um1 Þ

0

0

ðg þ q2 um2 Þ

0

0

0

ðd þ q3 um3 Þ

0

jg

0

1

1 ð1  q1 um Þ C C a C A 0

l

which gives a characteristic equation in k with two negative roots ða  q1 um1 Þ, ðd þ q3 um3 Þ and other two are given by

k2 þ kfl þ g þ q2 um2 g þ flðg þ q2 um2 Þ  jgð1  q1 um1 =aÞg ¼ 0: This implies (i) whenever (ii) whenever (iii) whenever

j < j0 , it has two roots with negative real parts, j ¼ j0 , it has one negative real root and one zero root, j > j0 , it has one root with negative real part and other root with positive real part.

Hence we have the following proposition: Proposition 2. The vanishing equilibrium E0 becomes stable, whenever u > a=q1 m1 . u < a=q1 m1 makes the vanishing equilibrium unstable and disease-free equilibrium E1 feasible. Moreover, the disease-free equilibrium E1 is locally asymptotically stable only when j 2 ð1; j0 Þ, i.e., when E is not feasible. At j ¼ j0 , E1 is critically stable and for j > j0 , E becomes unstable and E is feasible in this case. The Variational matrix at predator-free equilibrium E is given by

0

as as 0 B v ðg þ q2 um2 Þ i B JE ¼ B @ 0 0 ðd þ q3 m3 uÞ þ c0i 0

jg

1 s s C C C; 0 A l

0

which gives the characteristic equation in W as



 h



W þ d þ q3 m3 u  c0i

i

W3 þ W2 d1 þ Wd2 þ d3 ¼ 0;

 Þ þ lðg þ q2 um2 Þ  sjg, d3 ¼ ½aðg þ q2 um2 þ v  Þl þ jgðv   asÞs. where d1 ¼ as þ ðg þ q2 um2 Þ þ l, d2 ¼ asðg þ q2 um2 þ l þ v Under the assumption

i < d þ q3 m3 u ; c0

ð19Þ 0

ci this equation has clearly one negative real root, namely, d þ cþk  and other three roots are given by the second part of the 0i above equation. Hence, in a similar way as in paper [21], we can prove in terms of m ¼ s, that

Proposition 3. A single Hopf-bifurcation occurs at m ¼ m0 2 ð0; 1Þ for decreasing m, i.e., the predator-free equilibrium E is locally asymptotically stable in ðm0 ; 1Þ and unstable in ð0; m0 Þ. Moreover, such value m ¼ m0 2 ð0; 1Þ is unique. Next, the Variational matrix at endemic equilibrium E is given by

0

as

B v B JE ¼ B @ 0 0

as ðp þ g þ q2 um2 Þ

s

0 

i

cp

p

jg

0

0 



1

s C C C; 0 A l

Using Lyapunov–LaSalle theorem, we can prove that

1386

S. Ghosh, D.K. Bhattacharya / Applied Mathematical Modelling 34 (2010) 1382–1395

Proposition 4. The interior equilibrium E is locally asymptotically stable if,

4laðp þ g þ q2 um2 Þ : j <   v  g 4as þ as ðp þ g þ q2 um2 Þ

ð20Þ 

In this case, the stable manifold W S for E is given by W S ¼ fðs; i; p; v Þ 2 R4þ : ðs < s ; i < i ; p > p ; v > v  Þ p < p ; v < v  Þg. Proposition 5. The positive equilibrium E of the system (5)–(8) enters into Hopf-bifurcation when

S  ðs > s ; i > i ;

j passes through j .

Observation 3. In the above discussion on local stability of the equilibrium E , we have shown that there exists a j 2 ðl; 1Þ such that for j < j , E is stable. Further increase of j, system bifurcates towards a periodic solution of small amplitude; i.e., an oscillatory nature in the system can be observed. Thus we can see that there exists a lower threshold viz., j0 (Proposition 1) of virus replication parameter j, whenever j exceeds j0 the system becomes endemic and there exists an upper threshold j (Proposition 4) as well. If j exceeds j , the system exhibits an oscillatory nature in the population level. This definitely indicates that the overdosage of chemical as well as viral pesticide destabilizes the system and affect the overall performance of integrated pest management programme. 3. Optimization of net profit So far we have seen that in a large range of different control parameters defined by suitable thresholds in the earlier propositions, four dimensional model of the insect–virus–predator system can be controlled in the desired manner. But during the process of applying different control measures like, spraying of pesticide as chemical means and release of egt-virus as biological means, there is an obvious question of incurring some cost and allied profit in the whole process. Thus the objective is to quantify the units to express the net profit during the given time of experiment. In other words, this amounts to construct an economic model out of the given dynamic model of pest control. In this case the problem reduces to an optimal control problem under two control parameters (i) u, the rate of spraying of chemical pesticide and (ii) j, the rate of release of viral pesticide. Our task is then to formulate an optimal policy to find out the restrictions on the economic parameters of the model, when it is endemic and stable. We rewrite the model (5)–(8) as follows:

ds ¼ asf1  ðs þ iÞg  v s  u1 ; dt

ð21Þ

di ¼ v s  gi  ip  u2 ; dt

ð22Þ

dp ¼ pðd  p þ c0 iÞ  u3 ; dt

ð23Þ

dv ¼ lv þ u4 ; dt

ð24Þ

where ui ’s are the total effects (gain/loss) to the growth equations of species sðtÞ; iðtÞ; pðtÞ and v ðtÞ, respectively, due to simultaneous release of virus particles and spraying off chemical pesticide. We still assume that j lies between j0 , j and u < umax , i.e., the system is endemic and stable under integrated control. 3.1. Formulation of the profit function The formulation of the profit function p ¼ pðs; i; p; v Þ (in dimensionless form), which is to be maximized over the set of control parameters, is explained as follows: Let c: cost per unit chemical pesticide and cS : cost due to laboratory preparation of egt virus. p1 : projected price out of crop saving due to killing of susceptible pest. p2 : projected price out of crop saving due to killing of infected pest. p3 : indirect price of unit predator measured in terms of price of crop saved via killing of pest. p4 : indirect price per unit amount of egt virus measured in terms of crop saved via killing of pest. D1: The net profit for killing of susceptible pest by the use of pesticide = projected profit due to killing of susceptible pest  over the cost of pesticide of amount m1 u on crop ¼ p1 u1  cm1 u ¼ p1  qc s u1 . 1

D2: The net profit for killing of infected pest by the use of pesticide = projected profit due to killing of infected pest over the  cost of pesticide of amount m2 u ¼ p2 u2  cm2 u ¼ p2  qc i u2 . 2

S. Ghosh, D.K. Bhattacharya / Applied Mathematical Modelling 34 (2010) 1382–1395

1387

D3: Net expenditure for using pesticide of amount m3 u on the predator over the loss incurred due to killing of predators  ¼ cm3 u  p3 u3 ¼ qc p  p3 u3 . 3

D4: Net profit for killing of susceptible pest by the use of viral pesticide = projected profit due to killing of susceptible pest over the cost of laboratory preparation of egt virus ¼ p4 u4  cS . Hence, from the above discussion the profit function

p¼



p is given by

     c c c u1 þ p2  u2  p1   p3 u3 þ p4 u4  cS : q1 s q2 i q3 p

ð25Þ

After the formulation of the profit function, our next task is to find out the functions ui ðtÞ which drive the dynamical sys tem (21)–(24) from its initial state to a steady state optimal solution ðs ; i ; p ; v  Þ so as to maximize the integral

Z

J¼

T

pðs; i; p; v ; ui ; tÞdt;

ð26Þ

0

where T is the total time of applying both measures on model (21)–(24). We may point out that this steady state is a singular extremal in this case, because the control variables ui appear linearly in the system of Eqs. (21)–(24) and objective function (25). Applying the Pontryagin’s Maximum Principle on the constructed Hamiltonian H we obtain the optimal steady state  solution ðs ; i ; p ; v  Þ and corresponding control vector u0i . Moreover, in this optimal control problem the main objective is to maximize J. Thus the Legendre condition requires that along the singular solution, the matrix ðaij Þ whose  generalized   @H typical element aij ¼ @u@ i D2 @u ði; j ¼ 1; 2; 3; 4Þ, must be negative semi-definite at ðs ; i ; p ; v  Þ. This condition of semij definiteness of the above matrix imposes some restrictions on the economic parameters of the model. This means that the cost and the allied profit in the process must have some necessary limitations to get the maximum profit out of the whole process. In this connection we may point out that in the application of optimal control theory, necessary conditions are more useful than the sufficient conditions. This is because, it is extremely difficult to apply sufficient conditions in a real world problem. There are, in fact several sets of necessary conditions in optimal control theory in ecological problems, which in turn reflects the complexity in ecological systems. For example, a set of necessary conditions in this regard may be found in the book of Goh [24]. The most useful set of necessary conditions for singular control consists of the generalized Legendre conditions. A most general form of this condition was obtained by Goh [25]. By applying this condition we have the following theorem: Theorem 1. Let the problem be to find out functions ui ðtÞ which drive the dynamical system (21)–(24) from initial state ðs0 ; i0 ; p0 ; v 0 Þ to a steady state optimal solution ðs ; i ; p ; v  Þ so as to maximize the integral

Z

J¼

T

pðs; i; p; v ; ui ; tÞdt;

ð27Þ

0

where T is the total time of applying both measures on model (21)–(24). Then, a necessary condition that the profit function locally maximum at this steady state is given by

M2 þ 4csv 3

q2 i

ap1 2 4acp1 sv N < ; 3 p3 q2 i O2 þ

ac2 p1 s2 2

ðq2 i Þ2

<

4cs 2

q2 i

ð28Þ OðM þ

1

p3

O2 N 2 Þ

ð29Þ



at ðs ; i ; p ; v  Þ where B1 ¼ p2  p1 , M ¼ ap1  qcvi2 , N ¼ c0 p3  p2 , O ¼ B1  qc i. 2

Proof. We transform the integral occurring in J by choosing

g 1 ¼ g 1 ðsÞ;

g 2 ¼ g 2 ðiÞ;

g 3 ¼ g 3 ðpÞ;

g 4 ¼ g 4 ðv Þ

such that g 1 ðsð0ÞÞ ¼ 0; g 2 ðið0ÞÞ ¼ 0; g 3 ðpð0ÞÞ ¼ 0; g 4 ðv ð0ÞÞ ¼ 0, with

 c s_ ; q1 s   d c _ ½g ðiÞ ¼ p2  i; dt 2 q2 i   d c _ ½g ðpÞ ¼ p3  p; dt 3 q3 p d ½g ðsÞ ¼ dt 1



p1 

d ½g ðv Þ ¼ p4 v_ : dt 4

p is

2

1388

S. Ghosh, D.K. Bhattacharya / Applied Mathematical Modelling 34 (2010) 1382–1395

Then

Z

T

½pðs; i; p; v ; ui ; tÞdtÞ    Z0 T  o c c n ¼ p1  v s  gi  i_ fasf1  ðs þ iÞg  v s  s_ g þ p2  q1 s q2 i  0   c 0 þ p3  fpðd  ep þ c iÞ  p_ g þ p4 flv þ v_ g  cS dt q3 p  Z T Z T d d d d Rðs; i; p; v ; tÞdt  ðg 1 ðsÞÞ þ ðg 2 ðiÞÞ þ ðg 3 ðpÞÞ  ðg 4 ðv ÞÞ dt ¼ dt dt dt dt 0 Z0 T ¼ Rðs; i; p; v ; tÞdt  ½g 1 ðsðTÞÞ þ g 2 ðiðTÞÞ þ g 3 ðpðTÞÞ  g 4 ðv ðTÞÞ;

J¼

0

where

Rðs; i; p; v ; tÞ ¼ A/ þ Bw þ Cd þ D  cS

ð30Þ

c c c ; B ¼ p2  ; C ¼ p3 p  ; D ¼ p4 lv ; q1 q2 i q3 / ¼ af1  ðs þ iÞg  v ; w ¼ v s  ip  gi; d ¼ d  ep þ c0 i: A ¼ p1 s 

Now we consider the Hamiltonian

H ¼ Rðs; i; p; v ; tÞ þ

3 X

ki ðGi  ui Þ þ k4 ðG4 þ u4 Þ;

ð31Þ

i¼1

where G1 ¼ asf1  ðs þ iÞg  v s, G2 ¼ v s  ip  gi, G3 ¼ pðd  ep þ c0 iÞ, G4 ¼ lv and ki are co-state variables to be determined suitably.For steady state solution, we have

Gi  ui ¼ 0;

i ¼ 1; 2; 3 and G4 þ u4 ¼ 0:

ð32Þ 

If we suppose that there exists for which ðs ; i ; p ; v Þ gives a steady state optimal solution of (21)–(24); then from Pontryagin’s maximum principle, it follows that at this steady state, we have ðu01 ; u02 ; u03 ; u04 Þ

4 X

@H @R k_ 1 ¼  ¼  @s @s

4 X

4 X

ð33Þ

ki

@Gi ; @i

ð34Þ

ki

@Gi ; @p

ð35Þ

ki

@Gi @v

ð36Þ

i¼1

@H @R ¼  k_ 4 ¼  @v @v

4 X i¼1



@Gi ; @s

i¼1

@H @R k_ 3 ¼  ¼  @p @p



ki

i¼1

@H @R ¼  k_ 2 ¼  @i @i



and

@H ¼0 @ui

ð37Þ

at ðu01 ; u02 ; u03 ; u04 Þ. This implies that ki ¼ 0. The control vector u appears linearly in the Hamiltonian. Therefore ðu01 ; u02 ; u03 ; u04 Þ is a singular control variable. Now, the necessary conditions of optimality further ensure that at this singular control

  @H D ¼ 0; @u  i   @ @H D ¼ 0; @uj @ui

ð38Þ ð39Þ

which, in turn shows that k_ i ¼ 0. This further implies through (33)–(37) that

@R ¼ 0; @s

@R ¼ 0; @i

@R ¼ 0; @p

@R ¼ 0; @v

which by (30) reduces to,

@R vc ¼ 2ap1 s  ap1 i  þ B1 v þ C 1 ¼ 0; @s q2 i

ð40Þ

S. Ghosh, D.K. Bhattacharya / Applied Mathematical Modelling 34 (2010) 1382–1395

@R cv s ¼ ap1 s þ 2 þ C 2 ¼ 0; @i q2 i

1389

ð41Þ

@R ¼ ðc0 p3  p2 Þi  2ep3 p þ C 3 ¼ 0; @p @R cs ¼ B1 s  þ C 4 ¼ 0; @v q2 i

ð42Þ ð43Þ 0

where B1 ¼ p2  p1 , C 1 ¼ ap1 þ qac , C 2 ¼ qac  gp2  cc , C 3 ¼ qc þ qec  dp3 , C 4 ¼ qc þ p4 l, all are positive. q3 1 1 2 3 1 Solving (40)–(43), we have the following equation in i:

a4 i4 þ a3 i3 þ a2 i2 þ a1 i þ a0 ¼ 0; a4 ¼

where

ðc0 p3  p2 Þ2 B31 q32 , 3 ¼ cp3 q22 B21 C 2 þ 3cq2 B1 ðc0 p3

a

ð44Þ

2 p3 q32 B21 ðap1 C 4 þ B1 C 2 Þ þ q22 B21 ðc0 p3  p2 Þq2 B1 C 3  3cðc0 p3  p2 Þ, 2  p2 Þcðc0 p3  p2 Þ  q2 B1 C 3 , 1 ¼ 2 c p3 q2 ðap1 C 4 þ 3B1 C 2 Þþ 2ap1 q2 C 4 Þ, 0 ¼ 2c2 q2 C 1 C 4  c3 C 3 ðc0 p3  p2 Þ  2 c3 p3 C 2 .



a2 ¼ 2acp1 q22 B1 C 4

ð1 þ 2p3 Þ  6 a  c2 ðc0 p3  p2 Þ3q2 B1 C 3  a  cðc0 p3  p2 Þ þ 2cq2 C 4 ðq2 B1 C 1 þ acp1 þ clC P kK 2 > 0, as r, the Now we have, a4 > 0, and the constant coefficient of (44), a0 < 0, since C 1 ; C 3 ; C 4 and C 2 ¼ qcr  np kK  qP S growth rate of the pest population is high and C P the conversion factor is low. However, by DesCartes’ rule of sign, the above equation has at least one positive root and consequently, under the above  conditions, we get the optimal steady state equilibrium or nontrivial bionomic equilibrium ðs ; i ; p ; v  Þ of the model (21)–(24) which determines the upper and lower threshold values of costs. Moreover, using these solutions and model (21)–(24), we value of the singular control variable ðu01 ; u02 ; u03 ; u04 Þ.  can  find out the 2 @H @ @2 R @2 R ¼  Now @u D ði; j ¼ 1; . . . ; 4Þ. So it remains to show that the matrix @s@i is negative semi-definite at @uj @s@i i ðs ; i ; p ; v  Þ. From (40)–(43) we have

0

2ap1

B B ap1 þ v c2 B q2 i J¼B B 0 @ B1  qc i

ap1 þ qv ci2

0

B1  qc i

vs  2c q i3

c0 p3  p2

cs q2 i2

c0 p3  p2

2p3

0

cs q2 i2

0

0

2

2

2

2

1 C C C C: C A

For local maximum of profit function p, this matrix J is to be negative semi-definite i.e., all principle minors of the above  matrix are negative and positive in alternative order at ðs ; i ; p ; v  Þ. On simplification this condition finally reduces to

M2 þ 4csv q2 i

3

ap1 2 4acp1 sv N < ; 3 p3 q2 i O2 þ

ac2 p1 s2 2

ðq2 i Þ2

<

  1 2 2 O M þ O N 2 p3 q2 i 4cs



at ðs ; i ; p ; v  Þ. These are precisely the conditions (28) and (29). Hence the theorem is proved. h 4. Numerical results and discussion In this section, we use a numerical experiment on the dimensionless form of the system (5)–(8) to verify the results and interpretations we have made throughout this paper. This numerical experiment and simulation under various choice of control parameters, i.e., spraying of pesticide u and release of genetically engineered virus j, evidences that our model predicts the actual feature for integrated management in agricultural pest control. We carry out the following numerical simulations of our model with parameters estimated as in the paper [21]. Through time series graph it is shown how the population of all species in this system evolves under different choice of control parameters, ‘j’ and ‘u’. Fig. 1 shows that under the release of viral pesticide only, though the infection process starts at the beginning, it tends to the disease-free equilibrium (1, 0, 0, 0) soon. This indicates that the production of virus particles at this lower value of ‘j’ is not high enough for persistence of the infection into the system. The analytical result in Section 2 (Proposition 1), however, also supports this phenomenon. At a bit higher value of j, the system becomes endemic. Due to the increase in the value of j, the virus development increases, and as a result the level of susceptible pest decreases and that of infected one increases. This we notice in the Fig. 2. But further increase in the value of j, is not economically viable; so we have to make a balance in the whole control process, i.e. an integrated approach is encouraged. Together with the release of viral pesticide, spraying of chemical pesticide is also allowed into the system in a balanced dose. Fig. 3 shows that the spraying of chemical pesticide simultaneously with the release of viral pesticide in a lower value, rapidly decreases the level of susceptible pest, compared to the Fig. 2, and the system approaches the equilibrium (0.96775, 0, 0, 0) which is ‘disease-free’ (as infected population is zero!). So, let the value of j be increased so that the system enters into endemic state which is needed to suppress the pest population. Keeping the value of u, same as in Fig. 3, increase in the value of j, raises the pest mortality (Fig. 4). Now further increase in the value of j enhances the pest mortality as the virus development increases (Fig. 5). So the system remains

1390

S. Ghosh, D.K. Bhattacharya / Applied Mathematical Modelling 34 (2010) 1382–1395

Fig. 1. Time series graph of all species populations, at (j ¼ 107 ) and u ¼ 0, i.e., no spraying of chemical pesticide. The black, red, green and blue curves depict the susceptible pest SðtÞ, infected pest IðtÞ, predator PðtÞ and virus VðtÞ, respectively. Clearly the solution approaches the asymptotically stable equilibrium ð1; 0; 0; 0Þ. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2. Time series graph of (A) susceptible and infected pest and that of (B) predator and virus, whenever j ¼ 1010 and u ¼ 0. Curve descriptions are as in Fig. 1. Due to the increase in j, the virus development increases (B) and as a result the level of susceptible pest decreases and that of infected pest increases (A).

S. Ghosh, D.K. Bhattacharya / Applied Mathematical Modelling 34 (2010) 1382–1395

1391

Fig. 3. Time series graph of all species populations whenever the system is under both the control measures viz. viral pesticide (j ¼ 107 ) as well as chemical insecticide (u ¼ 2:5). Curve descriptions are as in Fig. 1. Due to the spraying of chemical pesticide the level of susceptible pest as well as infected pest decreases compared to Fig. 2, though at the same time the predator population and virus population go to extinction. This happens, because the required amount of infected inoculums in the site for initiation of next infection cycle is not sufficient.

Fig. 4. Time series graph of (A) susceptible and infected pest and that of (B) predator and virus, whenever the system is under both the control measures (j ¼ 1010 ; u ¼ 2:5). Due to the increase in the value of j, virus development increases (lower part) and effectively the level of susceptible pest decreases and that of infected increases (A), compared to Figs. 5 and 3.

1392

S. Ghosh, D.K. Bhattacharya / Applied Mathematical Modelling 34 (2010) 1382–1395

endemic and the solution approaches asymptotically stable equilibrium. But the system enters into bifurcation and oscillates, when the virus replication parameter j goes beyond some threshold value. Fig. 6A shows that at a higher value of j, the system consisting of susceptible and infected pest, converges to an endemic equilibrium with susceptible and infected population density 0.115234, 0.16933, respectively. But a little increase in the value of j, changes the picture dramatically. The solution bifurcates towards a limit cycle Fig. 6B. Moreover, up to certain threshold of j, the time period of the cycle increases. This however, we can observe through the results in Section 2 (Proposition 5). Thus we see that a minimum value of j is required to start the infection process into the system and we have to increase the value of j so that the system remains stable and infection persists into the system. But due to more increase in the value of j, we see that the system moves from stable to unstable state. So the natural question arises regarding proper, i.e., feasible as well as economic, control of pest model by minimizing the use of viral pesticide and applying some additional amount of chemical pesticide. For this purpose we have the Fig. 7. This figure shows the behavior of the solution when the system is under the application of chemical insecticide together with the viral pesticide. Here we see that when there is no spraying of chemical pesticide but a high dose of viral pesticide is applied, the system oscillates with a larger time period, but in the same value of j, if we apply chemical pesticide u, we see that the system oscillates in a smaller period, and finally the system comes back into a stable situation if we increase the value of u. Thus application of chemical pesticide indeed minimizes the effect of oscillation in the system and makes the system stable. So we may conclude that the integrated pest management policy with minimum use of pesticide as chemical control and simultaneous but restricted release of viral pesticide as biological control is the most favorable one. In the present work, we have seen the effectiveness of viral infection as a biological control in agro-ecosystem. It is definitely possible to eradicate pest population rapidly through heavy amount of release of viral pesticide. Our aim is not to eradicate whole pest population, but to keep it under injury level and also to increase the biomass of crop to a maximum

Fig. 5. Evolution of populations is shown at higher value of viral pesticide (j ¼ 1040 ). Curves are described as in Fig. 4. Solution approaches towards asymptotically stable equilibrium. Note that increase in the value of ‘j’, enhances the pest mortality as the development of virus increases.

S. Ghosh, D.K. Bhattacharya / Applied Mathematical Modelling 34 (2010) 1382–1395

1393

Fig. 6. Evolution of susceptible and infected pests at j ¼ 1075 (A) where the solution converges to the equilibrium (SðtÞ ¼ 0:12917, IðtÞ ¼ 0:143525), and that at j ¼ 1080 (B) where the solution bifurcates towards a limit cycle. Thus oscillatory nature comes into the system at higher values of j.

Fig. 7. Behavior of the solution when the system is under the application of chemical insecticide in different amount, together with the viral pesticide. The large cycle (in black) indicates that the system oscillates whenever j ¼ 1080 and there is no spraying of pesticide i.e.,u ¼ 0. The cycle of smaller period (in red) depicts the system oscillates for j ¼ 1080 and u ¼ 2, the cycle of more smaller period (in green) depicts the system oscillates for j ¼ 1080 and u ¼ 4 and finally the system comes back into stable situation (in blue) whenever j ¼ 1080 and u ¼ 6. Thus application of chemical pesticide indeed minimizes the effect of oscillation in the system and makes the system stable. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

1394

S. Ghosh, D.K. Bhattacharya / Applied Mathematical Modelling 34 (2010) 1382–1395

level so that the economical viability of the whole control procedure is maintained. However, model analysis cannot reveal what should be the best proportion of virus replication parameter ‘j’ and chemical pesticide ‘u’ from economic point of view and at this stage the application of optimization theory is needed. In fact, we have shown by forming a suitable objective function and performing optimal analysis that one can find optimal biomass of the species whenever optimal control is applied to the system. The essence of management is to make optimal decision subject to the realistic constraints. This is none other than an optimization problem. However, in practice there are enormous difficulties in quantifying the variables, the objectives and the constraints in a given decision problem. Here, we have focused our attention on the task of formulating an optimal policy when a decision problem has already been defined in a mathematical form. In this regard, we may mention that the numerical solution of the optimal control problem is much harder than that for the standard optimal control problem. This is because of the presence of mixed end conditions in the control problem. In this case, the steady state optimization plays an important role. It has added advantage of being easily implemented to the real world problem. To keep the state of the system at such an optimal steady state,if possible, one can, employ a globally stable control policy whenever the state is displaced from the optimal steady state. Another important theoretical development in the 1960s was the derivation of necessary conditions for singular control [26,27]. Singular control occurs in an optimal control problem in which one or more of the control variables appear linearly in system dynamics and the objective function. In ecological problems, the concept of singular control is important because control variables often occur linearly in the system dynamics. References in this connection may be given by the work of [28–30]. However, we have performed the optimal analysis of the economic model under two control parameters and have found out the necessary conditions on the parameters to make the profit maximum. In this connection, we may point out that though we have performed optimal analysis of the multi species community under two control parameters, but the techniques that we have developed in the theorem (Theorem 1, Section 3), could be applied for any finite number of singular control in the multi species community. Acknowledgements Authors are thankful to the reviewers for their valuable comments and suggestions for improvement of the manuscript. References [1] K.A. Jones, B. Zelazny, U. Ketunuti, A. Cherry, D. Grazywacz, A survey of virus control of insect-pests in South-East Asia and Western Pacific, in: F.R.H. Hunter-Fujita, P.F. Entwistle, M.F. Evans, N.E. Crook (Eds.), Insect Virus and Pest Management, John Wiley and Sons, NY, USA, 1998, pp. 244–248. [2] K. Narayanan, Insect defence: its impact on microbial control oh insect pests, Curr. Sci. 86 (2004) 800–814. [3] L.G. Coping, Baculoviruses in Crop Protection, Agro-Business Report DS 85, Richmond, UK, 1993, pp. 1–3. [4] D.R. O’Reilly, L.K. Miller, Improvement of baculovirus pesticide by deletion of EGT gene, Biotechnology 9 (1991) 1086–1089. [5] D. Grazywacz, D.J. Mckinley, G. Moawad, A. Cherry, The in vivo production of (Spodoptera littoralis) nuclear polyhedrosis virus, J. Virol. Method 71 (1998) 115–122. [6] M.R. Bell, In vivo production of nuclear polyhedrosis virus utilizing tobacco budwarm and a multicellular larval rearing container, J. Enomol. Sci. 26 (1991) 69–75. [7] A.J. Cherry, R.J. Rabindra, M.A. Parnell, N. Geetha, J.S. Kennedy, D. Grazywacz, Field evaluation of (Helicoverpa armigera) nucleopolyhedrovirus formulation for chickpea pod borer, (H. armigera) (Hubn.) on chickpea (Cicer arietinum var) Shoba in Southern India, Crop Prot. 19 (2000) 51–60. [8] V. Kalia, S. Chaudhari, G.T. Gujar, Optimization of production of nuclepolyhedrovirus of (Helicoverpa armigera) throughout larval stages, Phytoparasitica 29 (2001) 1–6. [9] C.M. Kumar, N. Sathiah, R.J. Rabindra, Optimizing the time of harvest of nucleopolyhedrovirus infected (Spodoptera litura) (Fabricius) larvae under (in vivo) production system, Curr. Sci. 88 (2005) 1682–1684. [10] M. Shapiro, C.D. Owens, R.A. Bell, H.A. Wood, Simplified efficient sytem for (in vivo) mass production of gypsy moth nucleopolyhedrosis virus, J. Econ. Entomol. 74 (1981) 110–111. [11] S. Subramanian, G. Santharam, N. Sathiah, J.S. Kennedy, R.J. Rabindra, Influence of incubation temperature on productivity and quality of (Spodoptera litura) nucleopolyhedro virus, Biol. Control. 37 (2006) 367–374. [12] A.J. Cherry, M.A. Parnell, D. Grzywacz, K.A. Jones, The optimization of in vivo nuclear polyhedrosis virus production in Spodoptera exempta (Walker) and Spodoptera exigua (Hubner), J. Invertbr. Pathol. 70 (1997) 50–58. [13] T.R. Shieh, Industrial production of viral pesticides, Adv. Virus Res. 36 (1989) 315–343. [14] X. Sun, X. Sun, B. Bai, W.V.D. Werf, J.M. Vlak, Z. Hu, Production of polyhedral inclusion bodies from (Helicoverpa armigera) larvae infected with wild type and recombinant HaSNPV, Biocontrol. Sci. Technol. 15 (2005) 353–366. [15] M. Shapiro, J.L. Robertson, R.A. Bell, Quantitative and qualitative differences in gypsy moth (Lepidoptera: Lymantriidae) nucleopolyhedrosis virus produced in different-aged larvae, J. Econ. Entomol. 79 (1986) 1174–1177. [16] N. Sathiah, S. Palaniswamy, Project Report on Establishment of Biopesticides Production Unit for (Helicoverpa armigera) Nuclear Polyhedrosis Virus, Tamil Nadu Agricultural University Press, 2004. pp. 1–17. [17] D.K. Bhattacharya, S. Karan, On bionomic model of integrated pest management of a single pest population, J. Differ. Equat. Dyn. Syst. 12 (3/4) (2004) 301–330. [18] T.L. Vincent, Pest management programs via optimal control theory, Biometrics 31 (1) (1975) 1–10. [19] S. Bhattacharyya, D.K. Bhattacharya, An improved integrated pest management model under 2-control parameters (sterile male and pesticide), Math. Biosci. 209 (2007) 256–281. [20] S. Bhattacharyya, D.K. Bhattacharya, A more realistic approach to pest-management problem, Bull. Math. Biol. 69 (2007) 1277–1310. [21] S. Bhattacharyya, D.K. Bhattacharya, Pest control through viral diseases: mathematical modeling and analysis, J. Theor. Biol. 238 (2006) 177–197. [22] R.P. Jacques, O.N. Morris, Compatibility of pathogens with other methods of pest control and crop protection, in: H.D. Burges (Ed.), Microbial Control of Pests and Plant Diseases, Academic Press, London, UK, 1981. [23] S. Ghosh, S. Bhattacharyya, D.K. Bhattacharya, The role of viral infection in pest control: a mathematical study, Bull. Math. Biol. 69 (2007) 2649–2691. [24] B.S. Goh, Management and Analysis of Biological populations, Elsevier Scientific Publishing Company, Amsterdam, Oxford, NY, 1980. [25] B.S. Goh, J. Leitman, T.L. Vincent, Optimal control of preypredator systems, Math. Biosci. 19 (1974) 263–286.

S. Ghosh, D.K. Bhattacharya / Applied Mathematical Modelling 34 (2010) 1382–1395 [26] [27] [28] [29] [30]

D.J. Bell, D.E. Jacobson, Singular Control Problems, Academic Press, NY, 1975. H.J. Kelley, R.E. Kopp, H.G. Moyer, Singular extremals, in: G. Leitmann (Ed.), Topics in Optimization, Academic Press, NY, 1967, pp. 63–101. N.G. Becker, Control of a pest population, Biometrics 26 (1970) 365–375. C.W. Clark, Economically optimal policies for the utilization of biologically renewable resources, Math. Biosci. 12 (1971) 245–260. C.W. Clark, Mathematical Bioeconomics, second ed., Wiley, New York, 1990.

1395