- Email: [email protected]
Species extinction problem: genetic vs ecological factors
Applied Mathematical Modelling 25 (2001) 937±951
www.elsevier.com/locate/apm
Species extinction problem: genetic vs ecological factors Ranjit Kumar Upadhyay a, Vikas Rai b, S.R.K. Iyengar
c,*
a
Department of Applied Mathematics, Indian School of Mines, Dhanbad 826 004, India National Institute of Immunology, Aruna Asaf Ali Marg, New Delhi 110 067, India Department of Mathematics, Indian Institute of Technology, Haus Khas, New Delhi 110 016, India b
c
Received 19 May 1998; received in revised form 2 January 2001; accepted 12 February 2001
Abstract Conservation biologists have been facing an intriguing question: whether it is genetic or ecological factors which govern the ecological systems. In this paper, we have constructed a few model systems describing real ecological situations and analysed them using a methodology designed for the purpose. Simulation experiments suggest that both these factors should be given equal weightage in working out strategies for any conservation eort. We conclude that the complex ecosystems are safe places for species belonging to the higher life forms, especially, generalist predators. On the contrary, simple (small) ecosystems cannot harbour these species for long. Another useful observation is that the vertebrate predators should be preferred to their invertebrate counterparts while aiming at conserving endangered prey species. Ó 2001 Elsevier Science Inc. All rights reserved. Keywords: Species extinction problem; Mathematical modelling; Genetic factors; Ecological factors; Ecosystem function; Predator±prey communities
1. Introduction The inadequacies of the current management systems for endangered species become apparent when one considers the following fact. Between 1973 (the year in which enactment of Endangered species took place) and 1990 as few as 16 species were delisted [25]. Over the same period, as many as 26 species (either listed in IUCN list or under review to be listed) became extinct [25]. Such an unaordable level of extinction calls for an in-depth study of the conservation strategies and plans. As a ®rst step, it is necessary to study the basic principles on which these eorts were designed. Many factors determine the abundance or rarity of species. The biological information that is most crucial for the design of the conservation schemes has been the subject of many investigations [6±8,16,17]. However, no decisive conclusions could be drawn from these investigations. Most of the investigators advocated an approach that is either ecological or genetic in emphasis. The proponents of the population genetic approach feel that understanding the organization of genetic diversity should be the key to the persistence of species as the genetic variation is a
*
Corresponding author. E-mail address: [email protected] (S.R.K. Iyengar).
0307-904X/01/$ - see front matter Ó 2001 Elsevier Science Inc. All rights reserved. PII: S 0 3 0 7 - 9 0 4 X ( 0 1 ) 0 0 0 3 4 - 8
938
R.K. Upadhyay et al. / Appl. Math. Modelling 25 (2001) 937±951
requisite for evolutionary adaption [9,11]. The proponents of the other approach argue that biotic interactions and habitat requirements of a species should be in the centre of any sound conservation programme [3,4]. Much of the theory that is being used in conservation biology was not developed for management applications, instead, has been inherited from parent disciplines like community ecology and population biology [18]. With [19] has opined that the debate over the utility of the theory in conservation research has not been able to appreciate the fact that the problem does not lie with the theory as such, but with the application of theory and the failure to understand and address the underlying assumptions that may restrict its use for practical applications [13,20]. Keeping the above facts in mind, we have designed four model ecosystems (see Figs. 2±4) to investigate whether ecological factors are more important than genetic ones. In the present paper, we assume that the growth rate parameter of a species is a genetic parameter. It is based on the possibility that population genetic processes may in¯uence vital rates through changes in the number and organization of alleles within and among individuals. Theoretical investigations [5] have suggested that reduced heterozygosity can result in decreased population growth due to inbreeding depression (reduced viability and fecundity). Moreover, allele richness could contribute to the population growth through its eect on the ability of a species to respond to changes in its selective environment [10]. The rest of the parameters in the model systems are ecological in nature. We have suggested a methodology to study and compare the in¯uences of genetic and ecological factors on a species persistence. The model systems are decoupled into two subsystems and conditions of sound health (a stable limit cycle in the phase space) are obtained. We have performed two-dimensional (2D) parameter scans to identify the parameter regimes in which dierent species extinct. The parameters chosen for the 2D scans are the parameters which control the dynamics of the system. The basis for performing the 2D scans is the belief that changes in physical conditions may bring corresponding changes in at least two parameters at a time. In the 2D scans, phase space studies are carried out in two dimensions by ®xing a genetic parameter on the X-axis and exploring the whole range of ecological parameters on the Y-axis. The analysis is repeated by giving an increment to the value of the genetic parameter along the X-axis. Many combinations of the genetic and ecological parameters are studied. Model systems are described in Section 2. The methodology used is presented in Section 3. Results are summarized in Section 4. In Section 5, we present a discussion on the relative importance of genetic and ecological factors in determining the fate of a species. 2. The model systems The design principles of the model ecosystems considered here are described in Fig. 1. These ecosystems are designed by assembling the simplest ecological units, i.e., a predator±prey interaction. Two such units are taken and are networked with each other under dierent schemes. The rationale for designing such model networking systems is that the natural ecosystems are complex ecological networks which are assembled in a very intricate way. Four such models are designed for analysis so that the study may provide some conclusions independent of the model analysed. We now describe the four model ecosystems. Assume that there exist two ecological units as shown in Fig. 1. Prey 1 and prey 2 are similar as far as their genetic make-up is concerned. Predator 2 is dierent from predator 1 in the sense that it is a sexually reproducing species and is a generalist one, i.e., it can live on food other than its favourite food (prey 2). Thus, the two basic units involve biologically distinct predator species.
R.K. Upadhyay et al. / Appl. Math. Modelling 25 (2001) 937±951
939
Fig. 1. Relationship between food-chain species and pseudo-prey.
We now combine these two units so that the resulting system is more rich in complexity. An ecologically sound way of networking the two systems is given in Fig. 2. This networking assumes that predator 2 is a predator for predator 1. This is a realistic assumption and supports the earlier assertion that predator 2 is a generalist one. The system depicted in Fig. 2 can be mathematically modelled as follows: dH1 WP1 H1 a1 H1 b1 H12 ; dt H1 D dP1 W1 H1 P1 W2 P 1 P 2 a2 P1 ; dt H1 D1 P1 D2 dH2 H2 W3 H 2 P 2 AH2 1 ; dt K H2 D3 dP2 W4 P22 cP2 ; dt H2 P1
1a 1b 1c 1d
where H1 ; H2 ; P1 ; P2 are prey 1, prey 2, specialist predator and generalist predator, respectively. Parameter a1 is the rate of self-reproduction of species H1 ; b1 measures the intensity of competition among individuals of species H1 for space, food, etc. W = H1 D is the per capita rate of removal of species H1 by P1 , W is the maximum value that the function WH1 = H1 D can attain, D represents that value of H1 at which the per capita predation rate attains half of its maximum value W . Biologically, it measures the agility of the prey to evade attacks by the predator. It also signi®es the extent to which the environment protects prey species from being captured by the predators. P1 is a specialist predator, i.e., H1 is the only food for it. a2 is the rate at which P1 dies out in the absence of its prey H1 . W1 H1 = H1 D1 denotes the gain in the specialist predator population due to proportionate loss in its prey. W1 represents the conversion eciency of the predator P1 in relation to the prey population H1 . Therefore, P1 dies out exponentially in the absence of H1 . Similar explanations hold for the other ratios. A and K are, respectively, the rate of
Fig. 2. Schematic representation of predator±prey interactions described in model 1.
940
R.K. Upadhyay et al. / Appl. Math. Modelling 25 (2001) 937±951
Fig. 3. A real world ecosystem exemplifying models 1 and 2.
self-reproduction and carrying capacity for prey 2. The ®rst term in the evolution equation for P2 describes the growth of P2 due to self-reproduction. W3 has similar meaning as W. c is the growth rate of the generalist predator P2 due to sexual reproduction. The last term in the evolution equation (1d) describes how loss in species P2 depends on per capita availability of its preys (H2 and P1 ). W4 measures the loss of predator P2 due to scarcity of its favourite food (H2 ) at any time. An example of a terrestrial ecosystem exemplifying this interaction network is shown in Fig. 3. In this model, a1 ; A; c are the genetic parameters and b1 ; K are the ecological parameters. There is a possibility of interaction between prey 1 and prey 2 in the designed model. This closes the loop which starts at prey 1. Transforming this possibility into ecological terms gives us a new ecosystem which is dierent from the earlier one in the sense that this also includes interspeci®c competition between prey species. To include this competition, Eqs. (1a) and (1c) need to be modi®ed. They now become dH1 WP1 H1 a1 H1 b1 H12 b2 H1 H2 ; dt H1 D dH2 H2 W3 H2 P2 AH2 1 B2 H1 H2 : dt K H2 D3
2a 2b
Eqs. (2a), (1b), (2b) and (1d) describe model 2. The ecosystem depicted in Fig. 3 when food supply is not in abundance serves as a real world example for this model system. Let us now form a simple food chain from the model system 1. If we assume that the second prey species in Fig. 2 does not exist, then we obtain the system as dH1 WP1 H1 a1 H1 b1 H12 ; dt H1 D dP1 W1 H1 P1 W2 P 1 P 2 a2 P1 ; dt H1 D1 P1 D2 dP2 W3 P22 cP2 : dt P1
3a 3b 3c
These equations describe model 3. An example of the real world ecological situation described by these equations is given in Fig. 4. It may be noted that this ecosystem is simple in form but rich in biology.
Fig. 4. A real world food chain represented by model 3.
R.K. Upadhyay et al. / Appl. Math. Modelling 25 (2001) 937±951
941
The biological richness of model 3 can be further enhanced by assuming that predator 2 is a vertebrate predator. This assumption is valid as we have earlier taken this predator to be a sexually reproducing species. The representative three species system can be written by modifying Eq. (3b) as dP1 dt
a2 P1
W1 H1 P1 H1 D1
W2 P12 P2 ; P12 D2
4
where W2 P12 = P12 D2 is the per capita functional response of the vertebrate predator P2 and was ®rst introduced by Takahashi [15]. The ecological role of per capita functional response was well described in the classic text by May [14]. Model 4 consists of Eqs. (3a), (3c) and (4). The parameters a1 ; a2 and c are genetic and b1 is the ecological parameter. A real world example representing this model is the food chain consisting of, say, rodent as prey, snake as specialist predator and peacock as the generalist (vertebrate) predator. The model ecosystems 3 and 4 are of lower dimensionality but are rich in biological complexity. The main objective behind designing these two 3D systems is to compare the role of dimensionality with that of biological complexity in the functioning of an ecosystem. For most sexual species over most population densities, reproduction is determined primarily by female numbers. Because of the ``female-biased demography'' of the sexual species the population growth rate is linear in population size which is well below the carrying capacity. However, the growth of a sexually reproducing population is proportional to square of the number of individuals present, in situations when the population is at very low densities. In many species, population decline to low numbers may be due to other causes, e.g., inbreeding depression, etc. [1,2,12]. This is caused due to the modi®cation of the environment of organisms physically or chemically by social interaction or by density dependent mating success. One such example is the conditioning of the medium by some aquatic microorganisms, by releasing substances that stimulate growth of conspeci®cs. In very sparse populations, either it is very dicult to ®nd a mate or social interaction necessary for reproduction is lacking. We have analysed the variants of all the model systems designed by replacing the linear term by a quadratic one in the rate equation for the bottom predator in order to know to what extent the Allee eect (an ecological phenomenon) in¯uences the persistence of species. This eect is observed in sexually reproducing species at very low population densities and is responsible for quadratic growth of the population due to reproduction. Recently, model systems constructed on the same principles have yielded interesting results [23,24]. The rationale for considering the reproduction of sexual species in two parts (at very low densities and the normal ones) is that no single analytic function is able to model it completely. Since the second subsystem containing the bottom predator should also form a Kolmogorov system (see Section 3), the last equation for the 4D model systems (models 1 and 2) is to be modelled as dP2 cP22 dt
W4 P22 ; H2 P1 D4
4a
and for the 3D systems it takes the form dP2 cP22 dt
W3 P22 : P1 D3
4b
Real ecosystems are organized in food webs [21,22]) which consist of chains with specialist predators on the top. These top predators are, in turn, predated by generalist predators. These facts
942
R.K. Upadhyay et al. / Appl. Math. Modelling 25 (2001) 937±951
have been the guiding principles in designing the model ecosystems which represent some typical real world ecological situations studied in the present paper. 3. The methodology We now derive the conditions for sound health of the ecosystems and study the eect of changes which may disturb their sound health. The present methodology assumes that two healthy constituent subsystems make a healthy complete system. Consider now, the four species ecosystems depicted in Fig. 2. We decouple the system into two constituent units (see Fig. 1). The evolution equations which describe these two units are dH1 WP1 H1 a1 H1 b1 H12 ; dt H1 D dP1 W1 H1 P1 a2 P1 ; dt H1 D1 and
dH2 H2 AH2 1 dt K 2 dP2 W4 P 2 cP2 : dt H2
5a 5b
W3 H 2 P 2 ; H2 D3
6a 6b
Eqs. (5a), (5b) constitute subsystem I and Eqs. (6a), (6b) constitute subsystem II. Each of the two systems can be written in the form dX XF X ; Y ; 7a dt dY YG X ; Y : 7b dt The relations among dierent system parameters ful®lling the demands of nine Kolmogorov conditions [14] (®ve constraints and four requirements) can now be obtained. It may be noted that some of the inequalities can be relaxed and replaced by equalities in certain situations [14]. According to the Kolmogorov theorem, if a system satis®es all the nine conditions simultaneously then it can admit either stable limit cycle or stable equilibrium point solutions. From the linear stability theory, we know that conditions on the parametric values can be obtained in a manner such that stable equilibrium point solutions are observed in the phase space of the system. If parametric values are chosen in such a way that those conditions obtained from the application of Kolmogorov theorem are satis®ed and the conditions of linear stability are violated, then the system yields stable limit cycle solution in the phase space. This corresponds to a situation of sound health for an ecosystem [14]. For the subsystem I, the local stability conditions are obtained as a2 D1 b1 D a1 > 0: D1 > 0; and 2b1 W1 a 2 The Kolmogorov constraints give rise to the following inequalities: W1 > a 2 ;
and a1 =b1 > D1 a2 = w1
a2 :
To perform numerical experiments, we use the following set of parameter values: a1 2; b1 0:05; W 1; D 10; a2 1; W1 2 and D1 10:
R.K. Upadhyay et al. / Appl. Math. Modelling 25 (2001) 937±951
943
In this case, the second local stability condition is violated. There exist many other sets of parameters satisfying the above requirements. We have chosen the above set for further study. For the subsystem II, Kolmogorov constraints are satis®ed automatically. The linear stability conditions for the subsystem II can be obtained as follows. The equilibrium point H2 ; P2 , which is the intersection of the isoclines dH2 =dt 0, and dP2 =dt 0 is given by q 1 cH H2 q q2 4A2 KD3 W42 ; q A K D3 W4 KcW3 ; and P2 2 : 2AW4 W4 If F1 ; F2 represent the right-hand sides of Eqs. (6a) and (6b), respectively, then the community matrix is obtained as a11 a12 ; M a21 a22 where
a11 a21
oF1 oH2 oF2 oH2
A 1
2H2 K
W4 P2 2 H2 2
;
W3 D3 P2
; 2
H2 D3 oF2 a22 c oP2
a12
oF1 oP2
2W4 P2 H2
W3 H2 ; H2 D3
W4 P2 H2
c:
The eigenvalues of the matrix M are the roots of k2 ak b 0;
with a
a11 a22 ; and b a11 a22
a12 a21 :
8
The subsystem is locally stable, if the eigenvalues are negative or have negative real parts. A necessary and sucient condition is a > 0 and b > 0. Simplifying a and b, we obtain these conditions as W3 D3 P2 2H2 c>A 1 ; 9 K H2 D3 2 " # W3 P2 W3 D3 P2 2H2 A 1 > 0: 10 H2 D3 K H2 D3 2 Subsystem II would admit stable limit cycle solutions when one or both of the above constraints, Eqs. (9) and (10), are violated. We ®nd that for the following values of the parameters, the subsystem II has a limit cycle solution: A 1:5; K 100; W3 0:74; D3 10; c 0:5; W4 0:2. For this set of values, we obtain the equilibrium point as H2 19:08; P2 47:7. In this case, Eq. (9) is not satis®ed. If the Allee eect is included in the subsystem II and the Kolmogorov theorem is applied, then either a stable equilibrium point or a stable limit cycle in its phase space is obtained, when the following constraints are satis®ed: c W4 = H2 D4 or c < W4 =D4 ;
11a
K > W4
11b
cD4 =c:
A Kolmogorov system can either admit a stable equilibrium or a stable limit cycle solution. Following the linear stability analysis, we ®nd that the subsystem II would be locally stable if:
944
R.K. Upadhyay et al. / Appl. Math. Modelling 25 (2001) 937±951
2H2 W3 D3 P2 > 0; 1 K H2 D3 2 W3 W4 H2 P2 > 0; 2 H2 D3 H2 D4
A
12a 12b
where H2 and P2 are the equilibrium points H2
1 W4 c
cD4 ;
and
P2
H2
D3 1
H2 K
A : W3
When the values of the parameters of this subsystem are chosen in such a way that the constraints (11a) and (11b) are satis®ed and one or both of the inequalities (12a) and (12b) are violated, then the subsystem admits a limit cycle solution. We ®nd that for the following values of the parameters the subsystem has a limit cycle solution: A 1;
K 50;
W3 1;
D3 20;
c 0:0062;
W4 0:2;
D4 20:
For this set of values, we obtain the equilibrium point as H2 12:26, and P2 24:35: In this case, Eq. (12a) is not satis®ed. There exist many other sets of parameter values in all the above cases. Let us now couple the two units. Essentially, this amounts to ®nding out the values of W2 and D2 . Since the parameters D1 and D2 play similar roles, their values may be taken to be same. W2 can be varied between 0.1 and 1. This choice of W2 is guided by the values of W ; W1 and W3 which were obtained from the application of the Kolmogorov theorem and the linear stability theory. Simulation experiments of the ®rst 4D system (model 1) are done using the above parameter values. The same methodology is used for studying the second 4D system. The subsystems in this case are identically same as for the model system 1. Therefore, the same parametric values are used for computations. In the case of the 3D systems, we assume that there exists a second prey (other than predator 1, called as a pseudo-prey) for predator 2. It is a realistic assumption as it is in conformity with the biology of all the other species in the ecosystem. The two subsystems obtained in this case are identical with those described by Eqs. (5a), (5b) and (6a), (6b), respectively . Application of the Kolmogorov theorem and the local stability criteria yields sets of parametric values for the simulation experiments. It may be mentioned that in the ®nal selection of parametric values, pseudo-prey terms are omitted. Parametric values of W2 and D2 are chosen in the same way as before. Since the second 3D system diers from the ®rst one in only one essential way (predator 2 is now a vertebrate), the same methodology is used for its analysis. 4. Results We have performed 2D parameter scans to identify the parameter regimes in which dierent species extinct. The parameters chosen for the 2D scans are the parameters which control the dynamics of the system. The basis for performing the 2D scans is the belief that changes in physical conditions may bring corresponding changes in at least two parameters at a time. In the 2D scans, phase space studies are carried out in two dimensions by ®xing a genetic parameter on the X -axis and exploring the whole range of ecological parameters on the Y -axis. The analysis is repeated by giving an increment to the value of the genetic parameter along the X -axis. Many combinations of the genetic and ecological parameters are studied. An ecosystem is considered to
R.K. Upadhyay et al. / Appl. Math. Modelling 25 (2001) 937±951
945
be in good health when all its constituent species coexist on a stable limit cycle (slc) attractor in the phase space. All the model systems are solved using the ODE Workbench package (AIP, New York). Twodimensional simulations have been done in the following cases. All the simulations were performed in the screen area 100 6 X 6 100 100 6 Y 6 100. Model 1: The following are the parameters in this model. Genetic parameters: · a1 : rate of self-reproduction for the prey species H1 , · A: rate of self-reproduction for the prey species H2 , · c: linear growth rate of the generalist predator P2 . Ecological parameters: · b1 : strength of intraspeci®c competition which measures the intensity of competition among individuals of prey species H1 for space, food, etc. · K: carrying capacity of its environment (for prey 2). Application of the Kolmogorov theorem and the local stability criteria yields sets of parametric values for the simulation experiments. Six groups of genetic vs ecological parameters can be formed for the simulation studies. The parameter values which are common for all the simulation experiments are W 1, D 10, A2 1, W1 2, D1 10, W2 0:5, D2 10, W3 0:74, D3 10 and W4 0:2. The following combinations of the genetic and ecological parameters are studied: (i) a1 ; b1 : 0:25 6 a1 6 2:5; 0:01 6 b1 6 2:0: The values of the other parameters are A 1:5; c 0:5; K 100. The ranges of a1 and b1 are divided into 9 and 49 subintervals, respectively. From the simulation studies, it is found that for all values of the parameters in the given ranges, P1 (the specialist predator) becomes extinct and the other species H2 ; P2 exist on the slc attractor in the phase plane. (ii) a1 ; K : 0:25 6 a1 6 2:5; 50 6 K 6 200: The values of the other parameters are A 1:5; c 0:5; b1 0:05: The ranges for a1 and K are divided into 9 and 6 subintervals, respectively. It is found that for all values of the parameters in the given ranges, P1 (the specialist predator) becomes extinct and the other species H2 ; P2 rest on a stable focus for smaller values of K 50; 75 and for other values of K, these rest on the stable limit cycle attractor in the phase plane. To study the eect of the parameter W2 (maximum value of the per capita functional response of the invertebrate generalist predator P2 ), a small value of W2 0:05 is chosen. In this case, it is found that P1 becomes extinct for a1 2 0:25; 0:75; while H2 (prey 2) becomes extinct for the remaining values of a1 . The other species H1 rests on a stable focus for all values of K and a1 . (iii) A; b1 : 0:25 6 A 6 2:5; 0:03 6 b1 6 2. The values of the other parameters are W2 0:5; a1 2; c 0:5 and K 100. In this case, P1 (the specialist predator) and H2 (prey 2) become extinct and the other species rests on a stable focus for A 2 0:25; 1:25: For A 1:5; 1:75 and for all b1 ; P1 becomes extinct and H2 ; P2 exist on slc (sound health condition). For A 2 2; 2:5 and for all b1 ; P1 becomes extinct and H2 ; P2 exist on a stable focus. (iv) A; K : 0:25 6 A 6 2:5; 50 6 K 6 200: The values of the other parameters are a1 2; b1 0:05; W2 0:5 and c 0:5. For A 2 0:25; 0:75; H2 (prey 2) becomes extinct for all values of K and the other species H1 ; P1 rest on a stable focus. For A 2 1; 2:5; P1 (specialist predator) becomes extinct and the other species H2 ; P2 rest on a stable focus and slc attractor, respectively. A test run was also made with a small value of W2 0:05. In this case, for all values of A and K in their ranges, H2 becomes extinct and H1 ; P1 rest on slc. (v) c; b1 : 0:25 6 c 6 2; 0:01 6 b1 6 2:5. The values of the other parameters are a1 2; A 1:5; K 100 and W2 0:5. For all values of c and b1 in the given ranges, P1
946
R.K. Upadhyay et al. / Appl. Math. Modelling 25 (2001) 937±951
(specialist predator) becomes extinct and the other species H2 ; P2 rest on a stable focus except for c 0:5. In this case, these species rest on the slc attractor. (vi) c; K : 0:25 6 c 6 2; 50 6 K 6 200. The values of the other parameters are a1 2; b1 0:05; W2 0:5 and A 1:5. For all values of c and K in the given ranges, P1 becomes extinct and the other species H2 ; P2 rest on a stable focus except at the value c 0:5: For c 0:5 and 100 6 K 6 200; H2 ; P2 rest on the slc attractor, indicating sound health conditions. Again, a test run was made with W2 0:05. In this case, for c 2 0:25; 1, H2 becomes extinct and H1 ; P1 rest on an slc attractor. For c 2 1:25; 2; H2 becomes extinct and H1 ; P1 rest on a stable focus. Model 2: As in model 1, a1 ; A and c are the genetic parameters with the same meanings. Ecological parameters: · b2 ; B2 : decay rates of the prey species H1 and H2 , respectively, due to interspeci®c competition between prey species. · K: carrying capacity of the environment of prey 2. Nine groups of genetic vs ecological parameters can be formed. Simulation experiments are done in all the cases. The values of the parameters that are common for all experiments are b1 0:05;
W 1;
D 10;
a2 1;
W1 2;
D2 10;
W3 0:74;
D3 10 and W4 0:2: The following ranges of parameters were used to form groups of genetic vs ecological parameters: 0:25 6 a1 6 2:5; 0:25 6 A 6 2:5;
0:0005 6 b2 6 0:5;
0:0005 6 B2 6 0:5;
25 6 K 6 250;
0:25 6 c 6 2:5:
The results are presented in Tables 1±3. From the tables, it is found that the specialist predator P1 is the most likely candidate to become extinct. Possible extinction of the species H1 (prey 1) is observed in the groups a1 ; b2 ; c; b2 while extinction of the species H2 (prey 2) is observed in the group A; B2 . In the remaining groups, the species behave in a similar way. Model 3: Genetic parameters: a1 ; c are as de®ned in model 1 and a2 is the rate at which the specialist predator P1 dies out in the absence of its prey H1 . Ecological parameter: b1 as in model 1. Three groups of genetic vs ecological parameters can be found. The parametric values which are common for all the simulation experiments are W 1;
D 10;
W1 2;
D1 10;
W2 0:5;
D2 10
and W3 0:2:
The results are presented in Figs. 5±7. In the ®gures, x1 is the prey 1 (H1 ), x2 is the specialist predator (P1 ) and x3 is the generalist predator (P2 ). Points shown in these ®gures represent extinction events. Extinction events corresponding to species H1 are located in the parametric region which is spanned by higher values of b1 and lower values of a1 . The top predator (generalist predator) is the most likely candidate to become extinct. The specialist predators are least prone to extinction. Less severe intraspeci®c competition among prey species and higher rate of reproduction of the top predators favour extinction of these species. Model 4: As in model 3, a1 ; c; a2 are the genetic parameters and b1 is the ecological parameter. Three groups of genetic vs ecological parameters are used for the simulation studies. The results are presented in Figs. 8±10. In the ®gures, x2 and x3 mean the specialist predator P1 and generalist predator P2 , respectively. Points shown in these ®gures represent extinction events. It is interesting to note that the prey species H1 and H2 always survive and never become extinct.
R.K. Upadhyay et al. / Appl. Math. Modelling 25 (2001) 937±951
947
Table 1
Simulation experiments for the group a1 ; b2 a Genetic parameter a1
Ecological parameter b2
Outcome
0.25±0.5
0.0005±0.01 0.03±0.07 0.09±0.5 0.0005±0.03 0.25±0.5 0.05±0.2 0.0005±0.05 0.25±0.45 0.07±0.2 0.5 0.0005±0.07 0.3±0.5 0.09±0.25 0.0005±0.1 0.3±0.5 0.15±0.25 0.0005±0.1 0.15±0.5 0.0005±0.1 0.3±0.5 0.15±0.25
P1 becomes extinct H1 ; P1 becomes extinct H1 becomes extinct P1 becomes extinct P1 becomes extinct H1 ; P1 become extinct P1 becomes extinct P1 becomes extinct H1 ; P1 become extinct H1 becomes extinct P1 becomes extinct P1 becomes extinct H1 ; P1 become extinct P1 becomes extinct P1 becomes extinct H1 ; P1 become extinct P1 becomes extinct H1 becomes extinct P1 becomes extinct P1 becomes extinct H1 ; P1 become extinct
0.75±1 1.25
1.5±1.75 2 2.25 2.5
a
The parametric values which were common in all the experiments are: W 1; D 10; a2 1; W1 2; W2 0:5; b1 0:05; D1 10; D2 10; W3 0:74; D3 10; and W4 0:2: The ®xed values for the other parameters are: A 1:5; K 100; B2 0:005 and c 0:5. Table 2
Simulation experiments for the group A; B2 a Genetic parameter A
Ecological parameter B2
Outcome
0.25±0.75 1
0.0005±0.5 0.0005±0.005 0.006±0.5 0.0005±0.01 0.02±0.5 0.0005±0.02 0.03±0.5 0.0005±0.03 0.04±0.5 0.0005±0.04 0.05±0.5
H2 becomes extinct P1 becomes extinct H2 becomes extinct P1 becomes extinct H2 becomes extinct P1 becomes extinct H2 becomes extinct P1 becomes extinct H2 becomes extinct P1 becomes extinct H2 becomes extinct
1.25±1.5 1.75 2±2.25 2.5 a
The parametric values which were common in all the experiments are: W 1, D 10, a2 1, W1 2, W2 0:05, b1 0:05, D1 10, D2 10, W3 0:74, D3 10; and W4 0:2: The ®xed values for the other parameters are: a1 2; b2 0:005; K 100; and c 0:5.
Lower values of the parameters a1 and b1 lead to extinction of the middle predator. But higher intraspeci®c competition among individuals of prey species b1 in combination with lower intrinsic death rates a2 favours extinction of the middle predator. The top predator is the species which faces the highest risk of extinction. A comparison of Figs. 5±7 and Figs. 8±10 suggests that the vertebrate predators feeding on alternative preys should be preferred to their invertebrate counterparts if preserving the prey species is the goal to be achieved. 5. Conclusions In this paper, we have attempted to study the in¯uence of genetic and ecological factors on the species extinction. We have designed four model systems to study these aspects. Simulation
948
R.K. Upadhyay et al. / Appl. Math. Modelling 25 (2001) 937±951
Table 3
Simulation experiments for the group c; b2 a Genetic parameter c
Ecological parameter b2
Outcome
0.25
0.0005±0.03 0.05±0.1 0.15±0.5 0.0005±0.15 0.2±0.5 0.0005±0.09 0.1±0.5 0.0005±0.15 0.35±0.5 0.0005±0.3 0.45±0.5 0.0005±0.45 0.0005±0.5
P1 becomes extinct H1 becomes extinct P1 becomes extinct P1 becomes extinct H1 becomes extinct P1 becomes extinct H1 becomes extinct P1 becomes extinct H1 becomes extinct P1 becomes extinct H1 becomes extinct P1 becomes extinct P1 becomes extinct
0.5 0.75 1 1.25 1.5 1.75±2.5 a
The parametric values which were common in all the experiments are: W 1, D 10, a2 1, W1 2, W2 0:5, b1 0:05, D1 10, D2 10, W3 0:74, D3 10; and W4 0:2: The ®xed values for the other parameters are: a1 2; A 1:5; B2 0:005 and K 100.
Fig. 5. Model system 3. 2D scan diagram between a1 ; b1 parameter space with W 1, D 10, a2 1, c 0:5, W1 2, D1 10, W2 0:5, D2 10 and W3 0:2:
Fig. 6. Model system 3. 2D scan diagram in a2 ; b1 parameter space with a1 2; c 0:5 (remaining parameter values are same as in Fig. 5).
R.K. Upadhyay et al. / Appl. Math. Modelling 25 (2001) 937±951
949
Fig. 7. Model system 3. 2D scan diagram in c; b1 parameter space with a1 2; a2 1 (remaining parameter values are same as in Fig. 5).
Fig. 8. Model system 4. 2D scan diagram in a1 ; b1 parameter space with a2 1, c 0:5, W 1, D 10, W1 2, D1 10, D2 100, W3 0:2 and W2 0:5.
Fig. 9. Model system 4. 2D scan diagram in a2 ; b1 parameter space with a1 2; c 0:5 (remaining parameter values are same as in Fig. 8).
experiments using 2D parameter scans are carried out. An important conclusion is that the specialist predators face high risk of extinction depending on the complexity of the system. The generalist predators ®nd safe habitats in complex ecosystems. An interesting observation which emanates from the present study is that the interspeci®c competition between species has an asymmetric eect on the two participating species, i.e., it
950
R.K. Upadhyay et al. / Appl. Math. Modelling 25 (2001) 937±951
Fig. 10. Model system 4. 2D scan diagram in c; b1 parameter space with a1 2; a2 1 (remaining parameter values are same as in Fig. 8).
favours one over the other. Use of vertebrate predators is advisable when preservation of the prey species is the goal. In simple (small) ecosystems the top predator species face the highest risk of extinction. The general observation from the present studies is that genetic and ecological parameters have equal in¯uence on the persistence or extinction of species. Therefore, conservationists should consider that the genetic and ecological factors are equally important. Acknowledgements The authors are grateful to Physics Academic Software (AIP, New York) for providing the ODE Workbench package which was used for numerical computations. The authors express their sense of gratitude to the reviewers whose comments and suggestions have improved the paper. References [1] W.C. Allee, A.E. Emerson, O. Park, T. Park, K.P. Schmidt, Principles of Animal Ecology, Saunders, London, 1949. [2] H.G. Andrewartha, L.C. Birch, The Distribution and Abundance of Animals, University of Chicago Press, Chicago, 1954. [3] P.F. Brussard, The role of ecology in biological conservation, Ecological Applications 1 (1991) 6±12. [4] M.A. Burgman, H.R. Akcakaya, S.S. Loew, The use of extinction models for species conservation, Biological Conservation 43 (1988) 9±25. [5] D. Charlesworth, B. Charlesworth, Inbreeding depression and its evolutionary consequences, Annual Review of Ecology and Systematics 18 (1987) 237±268. [6] E. Duey, A.S. Watts, The Scienti®c Management of Animal and Plant Communities for Conservation, Blackwell Scienti®c Publications, Oxford, 1971. [7] I.R. Franklin, Evolutionary change in small populations, in: Conservation Biology: An Evolutionary Ecological Perspective, Sinauer, Sunderland, MA, 1980, pp. 135±150. [8] D.A. Falk, K. Holsinger, Genetics and Conservation of Rare Plants, Oxford University Press, New York, 1991. [9] J.L. Hamrick, M.J.W. Godt, D.A. Murawski, M.D. Loveless, Correlations between species traits and allozyme diversity: implications for conservation biology, in: Genetics and Conservation of Rare Plants, Oxford University Press, New York, 1991, pp. 75±86. [10] R.K. Koehn, T.J. Hilbish, The adaptive importance of genetic variation, American Scientist 75 (1987) 134±141. [11] R. Lande, G.F. Barrowclough, Eective population size, genetic variation and their use in population management, in: Viable Populations for Conservation, Cambridge University Press, Cambridge, 1987, pp. 87±124. [12] R. Lande, Genetics and demography in biological conservation, Science 241 (1988) 1455±1460. [13] J.H. Lawton, Ecological experiments with model systems, Science 269 (1995) 328±331.
R.K. Upadhyay et al. / Appl. Math. Modelling 25 (2001) 937±951
951
[14] R.M. May, Stability and Complexity in Model Ecosystems, Princeton University Press, Princeton, NJ, 1973. [15] F. Takahashi, Reproduction curve with two equilibrium points. A considaration on the ¯uctuation of insect population, Res. Pop. Ecol. 6 (1964) 28±36. [16] C.M. Schonewald-Cox, S.M. Chambers, B. MacBryde, L. Thomas, Genetics and Conservation: A Reference for Managing Wild Animal and Plant Populations, Benjamin±Cummings, London, 1983. [17] M.E. Soule', Where do we go from here? in: Viable Populations for Conservation, Sinauer, Sunderland, MA, 1987, pp. 175±183. [18] D. Simberlo, The contribution of population and community biology to conservation science, Annual Review of Ecology and Systematics 19 (1988) 473±511. [19] K.A. With, The theory of conservation biology, Conservation Biology 11 (6) (1997) 1436±1440. [20] D.F. Doak, L.S. Mills, A useful role for theory in conservation, Ecology 75 (1994) 615±626. [21] S.L. Pimm, Food Webs, Chapman & Hall, London, 1982. [22] S.L. Pimm, The Balance of Nature, University of Chicago Press, Chicago, IL, 1991. [23] R.K. Upadhyay, V. Rai, Why chaos is rarely observed in natural populations? Chaos, Solitons & Fractals 8 (12) (1997) 1933±1939. [24] R.K. Upadhyay, S.R.K. Iyengar, V. Rai, Chaos: An ecological reality? International Journal of Bifurcation and Chaos 8 (6) (1998) 1325±1333. [25] D.W. Schemske, B.C. Husband, M.H. Ruckelshaus, C. Goodwillie, I.M. Parker, J.G. Bishop, Evaluating approaches to the conservation of rare and endangered plants, Ecology 75 (3) (1994) 584±606.
www.elsevier.com/locate/apm
Species extinction problem: genetic vs ecological factors Ranjit Kumar Upadhyay a, Vikas Rai b, S.R.K. Iyengar
c,*
a
Department of Applied Mathematics, Indian School of Mines, Dhanbad 826 004, India National Institute of Immunology, Aruna Asaf Ali Marg, New Delhi 110 067, India Department of Mathematics, Indian Institute of Technology, Haus Khas, New Delhi 110 016, India b
c
Received 19 May 1998; received in revised form 2 January 2001; accepted 12 February 2001
Abstract Conservation biologists have been facing an intriguing question: whether it is genetic or ecological factors which govern the ecological systems. In this paper, we have constructed a few model systems describing real ecological situations and analysed them using a methodology designed for the purpose. Simulation experiments suggest that both these factors should be given equal weightage in working out strategies for any conservation eort. We conclude that the complex ecosystems are safe places for species belonging to the higher life forms, especially, generalist predators. On the contrary, simple (small) ecosystems cannot harbour these species for long. Another useful observation is that the vertebrate predators should be preferred to their invertebrate counterparts while aiming at conserving endangered prey species. Ó 2001 Elsevier Science Inc. All rights reserved. Keywords: Species extinction problem; Mathematical modelling; Genetic factors; Ecological factors; Ecosystem function; Predator±prey communities
1. Introduction The inadequacies of the current management systems for endangered species become apparent when one considers the following fact. Between 1973 (the year in which enactment of Endangered species took place) and 1990 as few as 16 species were delisted [25]. Over the same period, as many as 26 species (either listed in IUCN list or under review to be listed) became extinct [25]. Such an unaordable level of extinction calls for an in-depth study of the conservation strategies and plans. As a ®rst step, it is necessary to study the basic principles on which these eorts were designed. Many factors determine the abundance or rarity of species. The biological information that is most crucial for the design of the conservation schemes has been the subject of many investigations [6±8,16,17]. However, no decisive conclusions could be drawn from these investigations. Most of the investigators advocated an approach that is either ecological or genetic in emphasis. The proponents of the population genetic approach feel that understanding the organization of genetic diversity should be the key to the persistence of species as the genetic variation is a
*
Corresponding author. E-mail address: [email protected] (S.R.K. Iyengar).
0307-904X/01/$ - see front matter Ó 2001 Elsevier Science Inc. All rights reserved. PII: S 0 3 0 7 - 9 0 4 X ( 0 1 ) 0 0 0 3 4 - 8
938
R.K. Upadhyay et al. / Appl. Math. Modelling 25 (2001) 937±951
requisite for evolutionary adaption [9,11]. The proponents of the other approach argue that biotic interactions and habitat requirements of a species should be in the centre of any sound conservation programme [3,4]. Much of the theory that is being used in conservation biology was not developed for management applications, instead, has been inherited from parent disciplines like community ecology and population biology [18]. With [19] has opined that the debate over the utility of the theory in conservation research has not been able to appreciate the fact that the problem does not lie with the theory as such, but with the application of theory and the failure to understand and address the underlying assumptions that may restrict its use for practical applications [13,20]. Keeping the above facts in mind, we have designed four model ecosystems (see Figs. 2±4) to investigate whether ecological factors are more important than genetic ones. In the present paper, we assume that the growth rate parameter of a species is a genetic parameter. It is based on the possibility that population genetic processes may in¯uence vital rates through changes in the number and organization of alleles within and among individuals. Theoretical investigations [5] have suggested that reduced heterozygosity can result in decreased population growth due to inbreeding depression (reduced viability and fecundity). Moreover, allele richness could contribute to the population growth through its eect on the ability of a species to respond to changes in its selective environment [10]. The rest of the parameters in the model systems are ecological in nature. We have suggested a methodology to study and compare the in¯uences of genetic and ecological factors on a species persistence. The model systems are decoupled into two subsystems and conditions of sound health (a stable limit cycle in the phase space) are obtained. We have performed two-dimensional (2D) parameter scans to identify the parameter regimes in which dierent species extinct. The parameters chosen for the 2D scans are the parameters which control the dynamics of the system. The basis for performing the 2D scans is the belief that changes in physical conditions may bring corresponding changes in at least two parameters at a time. In the 2D scans, phase space studies are carried out in two dimensions by ®xing a genetic parameter on the X-axis and exploring the whole range of ecological parameters on the Y-axis. The analysis is repeated by giving an increment to the value of the genetic parameter along the X-axis. Many combinations of the genetic and ecological parameters are studied. Model systems are described in Section 2. The methodology used is presented in Section 3. Results are summarized in Section 4. In Section 5, we present a discussion on the relative importance of genetic and ecological factors in determining the fate of a species. 2. The model systems The design principles of the model ecosystems considered here are described in Fig. 1. These ecosystems are designed by assembling the simplest ecological units, i.e., a predator±prey interaction. Two such units are taken and are networked with each other under dierent schemes. The rationale for designing such model networking systems is that the natural ecosystems are complex ecological networks which are assembled in a very intricate way. Four such models are designed for analysis so that the study may provide some conclusions independent of the model analysed. We now describe the four model ecosystems. Assume that there exist two ecological units as shown in Fig. 1. Prey 1 and prey 2 are similar as far as their genetic make-up is concerned. Predator 2 is dierent from predator 1 in the sense that it is a sexually reproducing species and is a generalist one, i.e., it can live on food other than its favourite food (prey 2). Thus, the two basic units involve biologically distinct predator species.
R.K. Upadhyay et al. / Appl. Math. Modelling 25 (2001) 937±951
939
Fig. 1. Relationship between food-chain species and pseudo-prey.
We now combine these two units so that the resulting system is more rich in complexity. An ecologically sound way of networking the two systems is given in Fig. 2. This networking assumes that predator 2 is a predator for predator 1. This is a realistic assumption and supports the earlier assertion that predator 2 is a generalist one. The system depicted in Fig. 2 can be mathematically modelled as follows: dH1 WP1 H1 a1 H1 b1 H12 ; dt H1 D dP1 W1 H1 P1 W2 P 1 P 2 a2 P1 ; dt H1 D1 P1 D2 dH2 H2 W3 H 2 P 2 AH2 1 ; dt K H2 D3 dP2 W4 P22 cP2 ; dt H2 P1
1a 1b 1c 1d
where H1 ; H2 ; P1 ; P2 are prey 1, prey 2, specialist predator and generalist predator, respectively. Parameter a1 is the rate of self-reproduction of species H1 ; b1 measures the intensity of competition among individuals of species H1 for space, food, etc. W = H1 D is the per capita rate of removal of species H1 by P1 , W is the maximum value that the function WH1 = H1 D can attain, D represents that value of H1 at which the per capita predation rate attains half of its maximum value W . Biologically, it measures the agility of the prey to evade attacks by the predator. It also signi®es the extent to which the environment protects prey species from being captured by the predators. P1 is a specialist predator, i.e., H1 is the only food for it. a2 is the rate at which P1 dies out in the absence of its prey H1 . W1 H1 = H1 D1 denotes the gain in the specialist predator population due to proportionate loss in its prey. W1 represents the conversion eciency of the predator P1 in relation to the prey population H1 . Therefore, P1 dies out exponentially in the absence of H1 . Similar explanations hold for the other ratios. A and K are, respectively, the rate of
Fig. 2. Schematic representation of predator±prey interactions described in model 1.
940
R.K. Upadhyay et al. / Appl. Math. Modelling 25 (2001) 937±951
Fig. 3. A real world ecosystem exemplifying models 1 and 2.
self-reproduction and carrying capacity for prey 2. The ®rst term in the evolution equation for P2 describes the growth of P2 due to self-reproduction. W3 has similar meaning as W. c is the growth rate of the generalist predator P2 due to sexual reproduction. The last term in the evolution equation (1d) describes how loss in species P2 depends on per capita availability of its preys (H2 and P1 ). W4 measures the loss of predator P2 due to scarcity of its favourite food (H2 ) at any time. An example of a terrestrial ecosystem exemplifying this interaction network is shown in Fig. 3. In this model, a1 ; A; c are the genetic parameters and b1 ; K are the ecological parameters. There is a possibility of interaction between prey 1 and prey 2 in the designed model. This closes the loop which starts at prey 1. Transforming this possibility into ecological terms gives us a new ecosystem which is dierent from the earlier one in the sense that this also includes interspeci®c competition between prey species. To include this competition, Eqs. (1a) and (1c) need to be modi®ed. They now become dH1 WP1 H1 a1 H1 b1 H12 b2 H1 H2 ; dt H1 D dH2 H2 W3 H2 P2 AH2 1 B2 H1 H2 : dt K H2 D3
2a 2b
Eqs. (2a), (1b), (2b) and (1d) describe model 2. The ecosystem depicted in Fig. 3 when food supply is not in abundance serves as a real world example for this model system. Let us now form a simple food chain from the model system 1. If we assume that the second prey species in Fig. 2 does not exist, then we obtain the system as dH1 WP1 H1 a1 H1 b1 H12 ; dt H1 D dP1 W1 H1 P1 W2 P 1 P 2 a2 P1 ; dt H1 D1 P1 D2 dP2 W3 P22 cP2 : dt P1
3a 3b 3c
These equations describe model 3. An example of the real world ecological situation described by these equations is given in Fig. 4. It may be noted that this ecosystem is simple in form but rich in biology.
Fig. 4. A real world food chain represented by model 3.
R.K. Upadhyay et al. / Appl. Math. Modelling 25 (2001) 937±951
941
The biological richness of model 3 can be further enhanced by assuming that predator 2 is a vertebrate predator. This assumption is valid as we have earlier taken this predator to be a sexually reproducing species. The representative three species system can be written by modifying Eq. (3b) as dP1 dt
a2 P1
W1 H1 P1 H1 D1
W2 P12 P2 ; P12 D2
4
where W2 P12 = P12 D2 is the per capita functional response of the vertebrate predator P2 and was ®rst introduced by Takahashi [15]. The ecological role of per capita functional response was well described in the classic text by May [14]. Model 4 consists of Eqs. (3a), (3c) and (4). The parameters a1 ; a2 and c are genetic and b1 is the ecological parameter. A real world example representing this model is the food chain consisting of, say, rodent as prey, snake as specialist predator and peacock as the generalist (vertebrate) predator. The model ecosystems 3 and 4 are of lower dimensionality but are rich in biological complexity. The main objective behind designing these two 3D systems is to compare the role of dimensionality with that of biological complexity in the functioning of an ecosystem. For most sexual species over most population densities, reproduction is determined primarily by female numbers. Because of the ``female-biased demography'' of the sexual species the population growth rate is linear in population size which is well below the carrying capacity. However, the growth of a sexually reproducing population is proportional to square of the number of individuals present, in situations when the population is at very low densities. In many species, population decline to low numbers may be due to other causes, e.g., inbreeding depression, etc. [1,2,12]. This is caused due to the modi®cation of the environment of organisms physically or chemically by social interaction or by density dependent mating success. One such example is the conditioning of the medium by some aquatic microorganisms, by releasing substances that stimulate growth of conspeci®cs. In very sparse populations, either it is very dicult to ®nd a mate or social interaction necessary for reproduction is lacking. We have analysed the variants of all the model systems designed by replacing the linear term by a quadratic one in the rate equation for the bottom predator in order to know to what extent the Allee eect (an ecological phenomenon) in¯uences the persistence of species. This eect is observed in sexually reproducing species at very low population densities and is responsible for quadratic growth of the population due to reproduction. Recently, model systems constructed on the same principles have yielded interesting results [23,24]. The rationale for considering the reproduction of sexual species in two parts (at very low densities and the normal ones) is that no single analytic function is able to model it completely. Since the second subsystem containing the bottom predator should also form a Kolmogorov system (see Section 3), the last equation for the 4D model systems (models 1 and 2) is to be modelled as dP2 cP22 dt
W4 P22 ; H2 P1 D4
4a
and for the 3D systems it takes the form dP2 cP22 dt
W3 P22 : P1 D3
4b
Real ecosystems are organized in food webs [21,22]) which consist of chains with specialist predators on the top. These top predators are, in turn, predated by generalist predators. These facts
942
R.K. Upadhyay et al. / Appl. Math. Modelling 25 (2001) 937±951
have been the guiding principles in designing the model ecosystems which represent some typical real world ecological situations studied in the present paper. 3. The methodology We now derive the conditions for sound health of the ecosystems and study the eect of changes which may disturb their sound health. The present methodology assumes that two healthy constituent subsystems make a healthy complete system. Consider now, the four species ecosystems depicted in Fig. 2. We decouple the system into two constituent units (see Fig. 1). The evolution equations which describe these two units are dH1 WP1 H1 a1 H1 b1 H12 ; dt H1 D dP1 W1 H1 P1 a2 P1 ; dt H1 D1 and
dH2 H2 AH2 1 dt K 2 dP2 W4 P 2 cP2 : dt H2
5a 5b
W3 H 2 P 2 ; H2 D3
6a 6b
Eqs. (5a), (5b) constitute subsystem I and Eqs. (6a), (6b) constitute subsystem II. Each of the two systems can be written in the form dX XF X ; Y ; 7a dt dY YG X ; Y : 7b dt The relations among dierent system parameters ful®lling the demands of nine Kolmogorov conditions [14] (®ve constraints and four requirements) can now be obtained. It may be noted that some of the inequalities can be relaxed and replaced by equalities in certain situations [14]. According to the Kolmogorov theorem, if a system satis®es all the nine conditions simultaneously then it can admit either stable limit cycle or stable equilibrium point solutions. From the linear stability theory, we know that conditions on the parametric values can be obtained in a manner such that stable equilibrium point solutions are observed in the phase space of the system. If parametric values are chosen in such a way that those conditions obtained from the application of Kolmogorov theorem are satis®ed and the conditions of linear stability are violated, then the system yields stable limit cycle solution in the phase space. This corresponds to a situation of sound health for an ecosystem [14]. For the subsystem I, the local stability conditions are obtained as a2 D1 b1 D a1 > 0: D1 > 0; and 2b1 W1 a 2 The Kolmogorov constraints give rise to the following inequalities: W1 > a 2 ;
and a1 =b1 > D1 a2 = w1
a2 :
To perform numerical experiments, we use the following set of parameter values: a1 2; b1 0:05; W 1; D 10; a2 1; W1 2 and D1 10:
R.K. Upadhyay et al. / Appl. Math. Modelling 25 (2001) 937±951
943
In this case, the second local stability condition is violated. There exist many other sets of parameters satisfying the above requirements. We have chosen the above set for further study. For the subsystem II, Kolmogorov constraints are satis®ed automatically. The linear stability conditions for the subsystem II can be obtained as follows. The equilibrium point H2 ; P2 , which is the intersection of the isoclines dH2 =dt 0, and dP2 =dt 0 is given by q 1 cH H2 q q2 4A2 KD3 W42 ; q A K D3 W4 KcW3 ; and P2 2 : 2AW4 W4 If F1 ; F2 represent the right-hand sides of Eqs. (6a) and (6b), respectively, then the community matrix is obtained as a11 a12 ; M a21 a22 where
a11 a21
oF1 oH2 oF2 oH2
A 1
2H2 K
W4 P2 2 H2 2
;
W3 D3 P2
; 2
H2 D3 oF2 a22 c oP2
a12
oF1 oP2
2W4 P2 H2
W3 H2 ; H2 D3
W4 P2 H2
c:
The eigenvalues of the matrix M are the roots of k2 ak b 0;
with a
a11 a22 ; and b a11 a22
a12 a21 :
8
The subsystem is locally stable, if the eigenvalues are negative or have negative real parts. A necessary and sucient condition is a > 0 and b > 0. Simplifying a and b, we obtain these conditions as W3 D3 P2 2H2 c>A 1 ; 9 K H2 D3 2 " # W3 P2 W3 D3 P2 2H2 A 1 > 0: 10 H2 D3 K H2 D3 2 Subsystem II would admit stable limit cycle solutions when one or both of the above constraints, Eqs. (9) and (10), are violated. We ®nd that for the following values of the parameters, the subsystem II has a limit cycle solution: A 1:5; K 100; W3 0:74; D3 10; c 0:5; W4 0:2. For this set of values, we obtain the equilibrium point as H2 19:08; P2 47:7. In this case, Eq. (9) is not satis®ed. If the Allee eect is included in the subsystem II and the Kolmogorov theorem is applied, then either a stable equilibrium point or a stable limit cycle in its phase space is obtained, when the following constraints are satis®ed: c W4 = H2 D4 or c < W4 =D4 ;
11a
K > W4
11b
cD4 =c:
A Kolmogorov system can either admit a stable equilibrium or a stable limit cycle solution. Following the linear stability analysis, we ®nd that the subsystem II would be locally stable if:
944
R.K. Upadhyay et al. / Appl. Math. Modelling 25 (2001) 937±951
2H2 W3 D3 P2 > 0; 1 K H2 D3 2 W3 W4 H2 P2 > 0; 2 H2 D3 H2 D4
A
12a 12b
where H2 and P2 are the equilibrium points H2
1 W4 c
cD4 ;
and
P2
H2
D3 1
H2 K
A : W3
When the values of the parameters of this subsystem are chosen in such a way that the constraints (11a) and (11b) are satis®ed and one or both of the inequalities (12a) and (12b) are violated, then the subsystem admits a limit cycle solution. We ®nd that for the following values of the parameters the subsystem has a limit cycle solution: A 1;
K 50;
W3 1;
D3 20;
c 0:0062;
W4 0:2;
D4 20:
For this set of values, we obtain the equilibrium point as H2 12:26, and P2 24:35: In this case, Eq. (12a) is not satis®ed. There exist many other sets of parameter values in all the above cases. Let us now couple the two units. Essentially, this amounts to ®nding out the values of W2 and D2 . Since the parameters D1 and D2 play similar roles, their values may be taken to be same. W2 can be varied between 0.1 and 1. This choice of W2 is guided by the values of W ; W1 and W3 which were obtained from the application of the Kolmogorov theorem and the linear stability theory. Simulation experiments of the ®rst 4D system (model 1) are done using the above parameter values. The same methodology is used for studying the second 4D system. The subsystems in this case are identically same as for the model system 1. Therefore, the same parametric values are used for computations. In the case of the 3D systems, we assume that there exists a second prey (other than predator 1, called as a pseudo-prey) for predator 2. It is a realistic assumption as it is in conformity with the biology of all the other species in the ecosystem. The two subsystems obtained in this case are identical with those described by Eqs. (5a), (5b) and (6a), (6b), respectively . Application of the Kolmogorov theorem and the local stability criteria yields sets of parametric values for the simulation experiments. It may be mentioned that in the ®nal selection of parametric values, pseudo-prey terms are omitted. Parametric values of W2 and D2 are chosen in the same way as before. Since the second 3D system diers from the ®rst one in only one essential way (predator 2 is now a vertebrate), the same methodology is used for its analysis. 4. Results We have performed 2D parameter scans to identify the parameter regimes in which dierent species extinct. The parameters chosen for the 2D scans are the parameters which control the dynamics of the system. The basis for performing the 2D scans is the belief that changes in physical conditions may bring corresponding changes in at least two parameters at a time. In the 2D scans, phase space studies are carried out in two dimensions by ®xing a genetic parameter on the X -axis and exploring the whole range of ecological parameters on the Y -axis. The analysis is repeated by giving an increment to the value of the genetic parameter along the X -axis. Many combinations of the genetic and ecological parameters are studied. An ecosystem is considered to
R.K. Upadhyay et al. / Appl. Math. Modelling 25 (2001) 937±951
945
be in good health when all its constituent species coexist on a stable limit cycle (slc) attractor in the phase space. All the model systems are solved using the ODE Workbench package (AIP, New York). Twodimensional simulations have been done in the following cases. All the simulations were performed in the screen area 100 6 X 6 100 100 6 Y 6 100. Model 1: The following are the parameters in this model. Genetic parameters: · a1 : rate of self-reproduction for the prey species H1 , · A: rate of self-reproduction for the prey species H2 , · c: linear growth rate of the generalist predator P2 . Ecological parameters: · b1 : strength of intraspeci®c competition which measures the intensity of competition among individuals of prey species H1 for space, food, etc. · K: carrying capacity of its environment (for prey 2). Application of the Kolmogorov theorem and the local stability criteria yields sets of parametric values for the simulation experiments. Six groups of genetic vs ecological parameters can be formed for the simulation studies. The parameter values which are common for all the simulation experiments are W 1, D 10, A2 1, W1 2, D1 10, W2 0:5, D2 10, W3 0:74, D3 10 and W4 0:2. The following combinations of the genetic and ecological parameters are studied: (i) a1 ; b1 : 0:25 6 a1 6 2:5; 0:01 6 b1 6 2:0: The values of the other parameters are A 1:5; c 0:5; K 100. The ranges of a1 and b1 are divided into 9 and 49 subintervals, respectively. From the simulation studies, it is found that for all values of the parameters in the given ranges, P1 (the specialist predator) becomes extinct and the other species H2 ; P2 exist on the slc attractor in the phase plane. (ii) a1 ; K : 0:25 6 a1 6 2:5; 50 6 K 6 200: The values of the other parameters are A 1:5; c 0:5; b1 0:05: The ranges for a1 and K are divided into 9 and 6 subintervals, respectively. It is found that for all values of the parameters in the given ranges, P1 (the specialist predator) becomes extinct and the other species H2 ; P2 rest on a stable focus for smaller values of K 50; 75 and for other values of K, these rest on the stable limit cycle attractor in the phase plane. To study the eect of the parameter W2 (maximum value of the per capita functional response of the invertebrate generalist predator P2 ), a small value of W2 0:05 is chosen. In this case, it is found that P1 becomes extinct for a1 2 0:25; 0:75; while H2 (prey 2) becomes extinct for the remaining values of a1 . The other species H1 rests on a stable focus for all values of K and a1 . (iii) A; b1 : 0:25 6 A 6 2:5; 0:03 6 b1 6 2. The values of the other parameters are W2 0:5; a1 2; c 0:5 and K 100. In this case, P1 (the specialist predator) and H2 (prey 2) become extinct and the other species rests on a stable focus for A 2 0:25; 1:25: For A 1:5; 1:75 and for all b1 ; P1 becomes extinct and H2 ; P2 exist on slc (sound health condition). For A 2 2; 2:5 and for all b1 ; P1 becomes extinct and H2 ; P2 exist on a stable focus. (iv) A; K : 0:25 6 A 6 2:5; 50 6 K 6 200: The values of the other parameters are a1 2; b1 0:05; W2 0:5 and c 0:5. For A 2 0:25; 0:75; H2 (prey 2) becomes extinct for all values of K and the other species H1 ; P1 rest on a stable focus. For A 2 1; 2:5; P1 (specialist predator) becomes extinct and the other species H2 ; P2 rest on a stable focus and slc attractor, respectively. A test run was also made with a small value of W2 0:05. In this case, for all values of A and K in their ranges, H2 becomes extinct and H1 ; P1 rest on slc. (v) c; b1 : 0:25 6 c 6 2; 0:01 6 b1 6 2:5. The values of the other parameters are a1 2; A 1:5; K 100 and W2 0:5. For all values of c and b1 in the given ranges, P1
946
R.K. Upadhyay et al. / Appl. Math. Modelling 25 (2001) 937±951
(specialist predator) becomes extinct and the other species H2 ; P2 rest on a stable focus except for c 0:5. In this case, these species rest on the slc attractor. (vi) c; K : 0:25 6 c 6 2; 50 6 K 6 200. The values of the other parameters are a1 2; b1 0:05; W2 0:5 and A 1:5. For all values of c and K in the given ranges, P1 becomes extinct and the other species H2 ; P2 rest on a stable focus except at the value c 0:5: For c 0:5 and 100 6 K 6 200; H2 ; P2 rest on the slc attractor, indicating sound health conditions. Again, a test run was made with W2 0:05. In this case, for c 2 0:25; 1, H2 becomes extinct and H1 ; P1 rest on an slc attractor. For c 2 1:25; 2; H2 becomes extinct and H1 ; P1 rest on a stable focus. Model 2: As in model 1, a1 ; A and c are the genetic parameters with the same meanings. Ecological parameters: · b2 ; B2 : decay rates of the prey species H1 and H2 , respectively, due to interspeci®c competition between prey species. · K: carrying capacity of the environment of prey 2. Nine groups of genetic vs ecological parameters can be formed. Simulation experiments are done in all the cases. The values of the parameters that are common for all experiments are b1 0:05;
W 1;
D 10;
a2 1;
W1 2;
D2 10;
W3 0:74;
D3 10 and W4 0:2: The following ranges of parameters were used to form groups of genetic vs ecological parameters: 0:25 6 a1 6 2:5; 0:25 6 A 6 2:5;
0:0005 6 b2 6 0:5;
0:0005 6 B2 6 0:5;
25 6 K 6 250;
0:25 6 c 6 2:5:
The results are presented in Tables 1±3. From the tables, it is found that the specialist predator P1 is the most likely candidate to become extinct. Possible extinction of the species H1 (prey 1) is observed in the groups a1 ; b2 ; c; b2 while extinction of the species H2 (prey 2) is observed in the group A; B2 . In the remaining groups, the species behave in a similar way. Model 3: Genetic parameters: a1 ; c are as de®ned in model 1 and a2 is the rate at which the specialist predator P1 dies out in the absence of its prey H1 . Ecological parameter: b1 as in model 1. Three groups of genetic vs ecological parameters can be found. The parametric values which are common for all the simulation experiments are W 1;
D 10;
W1 2;
D1 10;
W2 0:5;
D2 10
and W3 0:2:
The results are presented in Figs. 5±7. In the ®gures, x1 is the prey 1 (H1 ), x2 is the specialist predator (P1 ) and x3 is the generalist predator (P2 ). Points shown in these ®gures represent extinction events. Extinction events corresponding to species H1 are located in the parametric region which is spanned by higher values of b1 and lower values of a1 . The top predator (generalist predator) is the most likely candidate to become extinct. The specialist predators are least prone to extinction. Less severe intraspeci®c competition among prey species and higher rate of reproduction of the top predators favour extinction of these species. Model 4: As in model 3, a1 ; c; a2 are the genetic parameters and b1 is the ecological parameter. Three groups of genetic vs ecological parameters are used for the simulation studies. The results are presented in Figs. 8±10. In the ®gures, x2 and x3 mean the specialist predator P1 and generalist predator P2 , respectively. Points shown in these ®gures represent extinction events. It is interesting to note that the prey species H1 and H2 always survive and never become extinct.
R.K. Upadhyay et al. / Appl. Math. Modelling 25 (2001) 937±951
947
Table 1
Simulation experiments for the group a1 ; b2 a Genetic parameter a1
Ecological parameter b2
Outcome
0.25±0.5
0.0005±0.01 0.03±0.07 0.09±0.5 0.0005±0.03 0.25±0.5 0.05±0.2 0.0005±0.05 0.25±0.45 0.07±0.2 0.5 0.0005±0.07 0.3±0.5 0.09±0.25 0.0005±0.1 0.3±0.5 0.15±0.25 0.0005±0.1 0.15±0.5 0.0005±0.1 0.3±0.5 0.15±0.25
P1 becomes extinct H1 ; P1 becomes extinct H1 becomes extinct P1 becomes extinct P1 becomes extinct H1 ; P1 become extinct P1 becomes extinct P1 becomes extinct H1 ; P1 become extinct H1 becomes extinct P1 becomes extinct P1 becomes extinct H1 ; P1 become extinct P1 becomes extinct P1 becomes extinct H1 ; P1 become extinct P1 becomes extinct H1 becomes extinct P1 becomes extinct P1 becomes extinct H1 ; P1 become extinct
0.75±1 1.25
1.5±1.75 2 2.25 2.5
a
The parametric values which were common in all the experiments are: W 1; D 10; a2 1; W1 2; W2 0:5; b1 0:05; D1 10; D2 10; W3 0:74; D3 10; and W4 0:2: The ®xed values for the other parameters are: A 1:5; K 100; B2 0:005 and c 0:5. Table 2
Simulation experiments for the group A; B2 a Genetic parameter A
Ecological parameter B2
Outcome
0.25±0.75 1
0.0005±0.5 0.0005±0.005 0.006±0.5 0.0005±0.01 0.02±0.5 0.0005±0.02 0.03±0.5 0.0005±0.03 0.04±0.5 0.0005±0.04 0.05±0.5
H2 becomes extinct P1 becomes extinct H2 becomes extinct P1 becomes extinct H2 becomes extinct P1 becomes extinct H2 becomes extinct P1 becomes extinct H2 becomes extinct P1 becomes extinct H2 becomes extinct
1.25±1.5 1.75 2±2.25 2.5 a
The parametric values which were common in all the experiments are: W 1, D 10, a2 1, W1 2, W2 0:05, b1 0:05, D1 10, D2 10, W3 0:74, D3 10; and W4 0:2: The ®xed values for the other parameters are: a1 2; b2 0:005; K 100; and c 0:5.
Lower values of the parameters a1 and b1 lead to extinction of the middle predator. But higher intraspeci®c competition among individuals of prey species b1 in combination with lower intrinsic death rates a2 favours extinction of the middle predator. The top predator is the species which faces the highest risk of extinction. A comparison of Figs. 5±7 and Figs. 8±10 suggests that the vertebrate predators feeding on alternative preys should be preferred to their invertebrate counterparts if preserving the prey species is the goal to be achieved. 5. Conclusions In this paper, we have attempted to study the in¯uence of genetic and ecological factors on the species extinction. We have designed four model systems to study these aspects. Simulation
948
R.K. Upadhyay et al. / Appl. Math. Modelling 25 (2001) 937±951
Table 3
Simulation experiments for the group c; b2 a Genetic parameter c
Ecological parameter b2
Outcome
0.25
0.0005±0.03 0.05±0.1 0.15±0.5 0.0005±0.15 0.2±0.5 0.0005±0.09 0.1±0.5 0.0005±0.15 0.35±0.5 0.0005±0.3 0.45±0.5 0.0005±0.45 0.0005±0.5
P1 becomes extinct H1 becomes extinct P1 becomes extinct P1 becomes extinct H1 becomes extinct P1 becomes extinct H1 becomes extinct P1 becomes extinct H1 becomes extinct P1 becomes extinct H1 becomes extinct P1 becomes extinct P1 becomes extinct
0.5 0.75 1 1.25 1.5 1.75±2.5 a
The parametric values which were common in all the experiments are: W 1, D 10, a2 1, W1 2, W2 0:5, b1 0:05, D1 10, D2 10, W3 0:74, D3 10; and W4 0:2: The ®xed values for the other parameters are: a1 2; A 1:5; B2 0:005 and K 100.
Fig. 5. Model system 3. 2D scan diagram between a1 ; b1 parameter space with W 1, D 10, a2 1, c 0:5, W1 2, D1 10, W2 0:5, D2 10 and W3 0:2:
Fig. 6. Model system 3. 2D scan diagram in a2 ; b1 parameter space with a1 2; c 0:5 (remaining parameter values are same as in Fig. 5).
R.K. Upadhyay et al. / Appl. Math. Modelling 25 (2001) 937±951
949
Fig. 7. Model system 3. 2D scan diagram in c; b1 parameter space with a1 2; a2 1 (remaining parameter values are same as in Fig. 5).
Fig. 8. Model system 4. 2D scan diagram in a1 ; b1 parameter space with a2 1, c 0:5, W 1, D 10, W1 2, D1 10, D2 100, W3 0:2 and W2 0:5.
Fig. 9. Model system 4. 2D scan diagram in a2 ; b1 parameter space with a1 2; c 0:5 (remaining parameter values are same as in Fig. 8).
experiments using 2D parameter scans are carried out. An important conclusion is that the specialist predators face high risk of extinction depending on the complexity of the system. The generalist predators ®nd safe habitats in complex ecosystems. An interesting observation which emanates from the present study is that the interspeci®c competition between species has an asymmetric eect on the two participating species, i.e., it
950
R.K. Upadhyay et al. / Appl. Math. Modelling 25 (2001) 937±951
Fig. 10. Model system 4. 2D scan diagram in c; b1 parameter space with a1 2; a2 1 (remaining parameter values are same as in Fig. 8).
favours one over the other. Use of vertebrate predators is advisable when preservation of the prey species is the goal. In simple (small) ecosystems the top predator species face the highest risk of extinction. The general observation from the present studies is that genetic and ecological parameters have equal in¯uence on the persistence or extinction of species. Therefore, conservationists should consider that the genetic and ecological factors are equally important. Acknowledgements The authors are grateful to Physics Academic Software (AIP, New York) for providing the ODE Workbench package which was used for numerical computations. The authors express their sense of gratitude to the reviewers whose comments and suggestions have improved the paper. References [1] W.C. Allee, A.E. Emerson, O. Park, T. Park, K.P. Schmidt, Principles of Animal Ecology, Saunders, London, 1949. [2] H.G. Andrewartha, L.C. Birch, The Distribution and Abundance of Animals, University of Chicago Press, Chicago, 1954. [3] P.F. Brussard, The role of ecology in biological conservation, Ecological Applications 1 (1991) 6±12. [4] M.A. Burgman, H.R. Akcakaya, S.S. Loew, The use of extinction models for species conservation, Biological Conservation 43 (1988) 9±25. [5] D. Charlesworth, B. Charlesworth, Inbreeding depression and its evolutionary consequences, Annual Review of Ecology and Systematics 18 (1987) 237±268. [6] E. Duey, A.S. Watts, The Scienti®c Management of Animal and Plant Communities for Conservation, Blackwell Scienti®c Publications, Oxford, 1971. [7] I.R. Franklin, Evolutionary change in small populations, in: Conservation Biology: An Evolutionary Ecological Perspective, Sinauer, Sunderland, MA, 1980, pp. 135±150. [8] D.A. Falk, K. Holsinger, Genetics and Conservation of Rare Plants, Oxford University Press, New York, 1991. [9] J.L. Hamrick, M.J.W. Godt, D.A. Murawski, M.D. Loveless, Correlations between species traits and allozyme diversity: implications for conservation biology, in: Genetics and Conservation of Rare Plants, Oxford University Press, New York, 1991, pp. 75±86. [10] R.K. Koehn, T.J. Hilbish, The adaptive importance of genetic variation, American Scientist 75 (1987) 134±141. [11] R. Lande, G.F. Barrowclough, Eective population size, genetic variation and their use in population management, in: Viable Populations for Conservation, Cambridge University Press, Cambridge, 1987, pp. 87±124. [12] R. Lande, Genetics and demography in biological conservation, Science 241 (1988) 1455±1460. [13] J.H. Lawton, Ecological experiments with model systems, Science 269 (1995) 328±331.
R.K. Upadhyay et al. / Appl. Math. Modelling 25 (2001) 937±951
951
[14] R.M. May, Stability and Complexity in Model Ecosystems, Princeton University Press, Princeton, NJ, 1973. [15] F. Takahashi, Reproduction curve with two equilibrium points. A considaration on the ¯uctuation of insect population, Res. Pop. Ecol. 6 (1964) 28±36. [16] C.M. Schonewald-Cox, S.M. Chambers, B. MacBryde, L. Thomas, Genetics and Conservation: A Reference for Managing Wild Animal and Plant Populations, Benjamin±Cummings, London, 1983. [17] M.E. Soule', Where do we go from here? in: Viable Populations for Conservation, Sinauer, Sunderland, MA, 1987, pp. 175±183. [18] D. Simberlo, The contribution of population and community biology to conservation science, Annual Review of Ecology and Systematics 19 (1988) 473±511. [19] K.A. With, The theory of conservation biology, Conservation Biology 11 (6) (1997) 1436±1440. [20] D.F. Doak, L.S. Mills, A useful role for theory in conservation, Ecology 75 (1994) 615±626. [21] S.L. Pimm, Food Webs, Chapman & Hall, London, 1982. [22] S.L. Pimm, The Balance of Nature, University of Chicago Press, Chicago, IL, 1991. [23] R.K. Upadhyay, V. Rai, Why chaos is rarely observed in natural populations? Chaos, Solitons & Fractals 8 (12) (1997) 1933±1939. [24] R.K. Upadhyay, S.R.K. Iyengar, V. Rai, Chaos: An ecological reality? International Journal of Bifurcation and Chaos 8 (6) (1998) 1325±1333. [25] D.W. Schemske, B.C. Husband, M.H. Ruckelshaus, C. Goodwillie, I.M. Parker, J.G. Bishop, Evaluating approaches to the conservation of rare and endangered plants, Ecology 75 (3) (1994) 584±606.