Depletion of Forest Resources and Wildlife Population with Habitat Complexity: A Mathematical Model

Rachana Pathak

Department of Mathematics, Faculty of Engineering, University of Lucknow (New Campus), Lucknow, India.

DOI: 10.4236/oje.2018.811034                

Abstract

Nowadays, management and regulation of natural resources like agriculture, fisheries, forestry and wildlife is one of the popular topics in research. The evolution of humankind is largely dependent on the quality of the environment and the resources it provides; but numerous human-induced factors, and climate change may drastically alter the conditions of human sustainability. This paper deals with effect of numerous human-induced activities on the depletion of forestry resources and wildlife population with habitat complexity. A nonlinear mathematical model is proposed and analyzed. In modeling process, we assume that the growth rate of wildlife population wholly depends on forestry biomass. It is depleted by human-induced activities. Local stability analysis of the mathematical model along with the persistence of the system is checked by using theory of nonlinear ordinary differential equations and Butler-McGhee lemma. Analytical results obtained are justified numerically through numerical simulation. Important parameters are investigated and variation of variables with change in these parameters is determined.

Keywords

Forestry Biomass, Depletion, Wildlife Population, Human-Induced Activities, Stability, Persistence

1. Introduction

India is known for its rich heritage of biological diversity, having already documented over 91,000 species of animals and 45,500 species of plants in its ten bio-geographic regions. Nearly 6500 native plants are still used prominently in indigenous healthcare systems. Thousands of locally-adapted crop varieties, grown traditionally since ancient times, and nearly 140 native breeds of farm livestock, continue to thrive in its diversified farming systems. The country is recognized as one of the eight Vavilovian Centres of Origin and Diversity of Crop Plants, having more than 300 wild ancestors and close relatives of cultivated plants still growing and evolving under natural conditions. But Forests are being destroyed and degraded at alarming rates. Deforestation comes in many forms, including fires, clear-cutting for agriculture, ranching and development, unsustainable logging for timber, and degradation due to climate change. This impacts people’s livelihoods and threatens a wide range of plant and animal species. Forests are disappearing at an alarming rate―18.7 million acres of forests are lost annually, equivalent to 27 soccer fields every minute.

Mathematical modeling is an essential tool in studying a diverse range of such ecological situations. Proposed work deals with the study of mathematical model on the sustainable management of renewable resources in India. The depletion of forest biomass by human population and industrialization has been investigated both theoretically as well as experimentally by many researchers ( [1] - [7] ). In particular, [8] [9] have proposed a mathematical model for the depletion of forest resources where they have shown that the forest resources may be doomed to extinction under growing population living in the forest or by industrialization which wholly depends on forest resources. [8] [9] results have also been found by other researchers showing their concern about biodiversity caused by deforestation. [10] has proposed and analyzed a mathematical model for conservation of forestry biomass and wildlife population. In above paper, we assumed that the growth rate of wildlife conservation is proportional to the depletion of forestry biomass due to wildlife population. We have also studied the affect of illegal trade in forestry biomass and wildlife population. [11] has proposed and analysed the mathematical model and shown the effect of population and industrialization on the forestry biomass and wildlife population. This model is motivated by our food chain model given by [12] . We proposed a nonlinear mathematical model and analyzed to see the effect of harvesting on food chain model with habitat complexity. We again study this model in the view of effect of human-induced activities on the forestry biomass and wildlife population with habitat complexity.

The rest of this paper is organized as follows: in Section 2, we introduce our mathematical model. In Section 3 and 4, our model is analyzed with regard to equilibria and their stabilities, respectively. In Section 5, we find the conditions for persistence of the system. In Section 6, we present a numerical example to illustrate the applicability of results obtained and also investigate the occurrence of Hopf bifurcation for certain set of parameters. We conclude with a short discussion in Section 7.

2. Mathematical Model

The mathematical model for effecting human-induced activities on the forestry biomass and wildlife population with habitat complexity can be described by the following system of nonlinear differential equations.

$\begin{array}{l}\frac{\text{d}B}{\text{d}t}=r_0\left(1-\frac{B}{K}\right)B-\frac{a_{1}B{W}_x}{b_1+B}-q_{1}HB,\\ \frac{\text{d}{W}_x}{\text{d}t}=\xi W_x\left(\frac{a_{1}B}{b_1+B}-\eta \right)-\frac{a_2\left(1-c\right)W_x{W}_y}{1+a_2\left(1-c\right)h{W}_x}-q_{2}H{W}_x,\\ \frac{\text{d}{W}_y}{\text{d}t}=\theta W_y\left(\frac{a_2\left(1-c\right)W_x}{1+a_2\left(1-c\right)h{W}_x}-\mu \right)-q_{3}H{W}_y.\end{array}$ (2.1)

$B\left(0\right)=B_0\ge 0,W_x\left(0\right)=W_x{}_0\ge 0,W_y\left(0\right)=W_y{}_0\ge 0,0<\xi ,\theta \le 1,0<c<1.$

The flowchart of the model (2.1) as shown in Figure 1.

Here, B be the forestry biomass, $W_x$ be the density of wild grazer population and $W_y$ be the density of wild predator population. It is assumed that the dynamics of forestry biomass follow logrithimic growth rate. It is further assumed that the rate of depletion of B due to human-induced activities is proportional to the product of forestry biomass and human-induced activities coefficient $q_1$ . As wild grazer population is wholly dependent on forestry biomass, its growth rate is proportional to the interaction term $\left(a_{1}B{W}_x/\left(b_1+B\right)\right)$ . The depletion rate of wild grazer population due to human activities is proportional to the product of concentration of $W_x$ and human activities coefficient $q_2$ . The growth rate of wild predator population is proportional to the interaction term (i.e. $a_2\left(1-c\right)W_x{W}_y/1+a_2\left(1-c\right)h{W}_x$ ). The natural depletion rate of wild predator population is considered proportional to $W_y$ . The depletion rate of wild predator population due to human activities is proportional to the product of concentration of $W_y$ and human activities coefficient $q_3$ . $r_0$ is initial growth rate of the forestry biomass, K is the carrying capacity of forestry biomass, $a_1$ is the attack rate of grazers, $b_1$ is the half-saturation constant, $\xi$ is the forestry biomass-grazer conversion rate, $\eta$ is the zero population growth grazer intake, $a_2$ is the saturation killing rate, $\theta$ is the prey (grazer)-predator conversion rate, $\mu$ is the zero population growth predator intake, h is the handling time and $0<c<1$ is a dimension less parameter that measures the degree or strength of habitat complexity which is effected by human-induced activities like industrialization, development, framing etc. It is to be noted, when $c=0$, i.e. when there is no complexity, we get back the original Holling Type II response function. It may be pointed out that for feasibility of the model (2.1), the growth rates of wild grazer population should be positive.

To analyze the model (2.1), we need the bounds of dependent variables involved. For this we find the region of attraction in the following lemma.

Lemma 2.1: The set $\Omega =\left\{\left(B,W_x,W_y\right):0\le B+W_x+W_y\le \frac{r_{0}K}{\alpha }\right\}$,

where

$\alpha =\min \left\{q_{1}H,\eta \xi +q_{2}H,\theta \mu +q_{3}H\right\}$

is the region of attraction for all solutions initiating in the interior of the positive octant.

Figure 1. Flowchart of mathematical model (2.1).

Proof: Let $\left(B\left(t\right),W_x\left(t\right),W_y\left(t\right)\right)$ be any solution with positive initial conditions $\left(B_0,W_x{}_0,W_y{}_0\right)$ .

Define a function

$V\left(t\right)=B\left(t\right)+W_x\left(t\right)+W_y\left(t\right).$

Computing the time derivative of $V\left(t\right)$ along the solutions of a system (2.1), we get

$\frac{\text{d}V}{\text{d}t}\le r_{0}K-\alpha V,$

where

$\alpha =\min \left\{q_{1}H,\eta \xi +q_{2}H,\theta \mu +q_{3}H\right\}$

Applying a theorem of differential inequalities, we obtain $0\le V\le \frac{r_{0}K}{\alpha }$ . Therefore all solutions of system (2.1) enter into the region

$\Omega =\left\{\left(B,W_x,W_y\right):0\le B+W_x+W_y\le \frac{r_{0}K}{\alpha }\right\}.$

This completes the proof of lemma.

3. Equilibrium Analysis

There exist following four equilibria of the system (2.1)

1) $E_0\left(0,0,0\right)$,

2) $E_1\left(B_1,0,0\right)$, where $B_1=\frac{K\left(r_0-q_{1}H\right)}{r_0}$, with condition $r_0-q_{1}H>0$,

3) $E_2\left(B_2,W_x{}_2,0\right)$ and

4) $E_3\left(B^{\ast },W_x^{\ast },W_y^{\ast }\right)$ .

Existence of $E_2\left(B_2,W_x{}_2,0\right)$

Here $B_2,W_{x2}$ are the positive solutions of the system of algebraic equations given below

$\begin{array}{l}{r}_0\left(1-\frac{B_2}{K}\right)-\frac{a_1{W}_{x2}}{b_1+B_2}-q_{1}H=0,\\ \xi \left(\frac{a_1{B}_2}{b_1+B_2}-\eta \right)-q_{2}H=0.\end{array}$ (3.1)

After solving (3.1), we have conditions of existence for equilibrium $E_2$ as follows

$\xi a_1-\left(\xi \eta +q_{2}H\right)>0,$

where

$\begin{array}{l}{B}_2=\frac{b_1\left(\xi \eta +q_{2}H\right)}{\xi a_1-\left(\xi \eta +q_{2}H\right)},\\ W_x{}_2=\frac{\left(r_0-q_{1}H\right)b_1}{a_1}+\frac{\left[K\left(r_0-q_{1}H\right)-r_0{b}_1\right]}{a_1}{B}_2-\frac{r_0{B}_2^2}{K{a}_1}.\end{array}$

Existence of $E_3\left(B^{\ast },W_x^{\ast },W_y^{\ast }\right)$ .

Here $B^{*},W_x^{*},W_y^{*}$ are the positive solutions of the system of algebraic equations given below

$\begin{array}{l}{r}_0\left(1-\frac{B^{*}}{K}\right)-\frac{a_1{W}_x^{*}}{b_1+B^{*}}-q_{1}H=0,\\ \xi \left(\frac{a_1{B}^{*}}{b_1+B^{*}}-\eta \right)-\frac{a_2\left(1-c\right)W_y^{*}}{1+a_2\left(1-c\right)h{W}_x^{*}}-q_{2}H=0,\\ \theta \left(\frac{a_2\left(1-c\right)W_x^{*}}{1+a_2\left(1-c\right)h{W}_x^{*}}-\mu \right)-q_{3}H=0.\end{array}$ (3.2)

After solving (3.2), we have conditions of existence for equilibrium $E_3$ as follows

$c<1,h<\frac{\theta }{q_{3}H+\mu \theta },\frac{\xi a_1{B}^{\ast }}{b_1+B^{\ast }}>\xi \eta +q_{2}H,$

where

$W_x^{\ast }=\frac{q_{3}H+\mu \theta }{a_2\left(1-c\right)\left\{\theta -\left(q_{3}H+\mu \theta \right)h\right\}}>0,$

$W_y^{\ast }=\frac{1}{q_2\left(1-c\right)}\left[\left\{\xi \left(\frac{a_1{B}^{\ast }}{b_1+B^{\ast }}-\eta \right)-q_{2}H\right\}\left\{1+a_2\left(1-c\right)h{W}_x^{*}\right\}\right],$

and $B^{*}$ is the unique positive root of the following equation

$r_0{B}^{*2}+\left(r_0{b}_1-r_{0}K+q_{1}HK\right)B^{*}+a_1{W}_x^{*}K-r_0{b}_{1}K+q_{1}H{b}_{1}K=0.$

4. Local Stability Analysis

To discuss the local stability of the system (2.1), we compute the variational matrix of the system (2.1). The entries of the general variational matrix are given by differentiating the right side of the system (2.1) with respect to $B,W_x$ and $W_y$, i.e.

$M\left(E_i\right)=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right]\text{.}$

where

$\begin{array}{l}{a}_{11}=r_0-\frac{2r_{0}B}{K}-\frac{b_1{a}_1{W}_x}{{\left(b_1+B\right)}^2}-q_{1}H,a_{12}=-\frac{a_{1}B}{b_1+B},a_{21}=\frac{b_1\xi W_x{a}_1}{{\left(b_1+B\right)}^2},\\ a_{23}=-\frac{a_2\left(1-c\right)W_x}{1+a_2\left(1-c\right)h{W}_x},a_{22}=\xi \left(\frac{a_{1}B}{b_1+B}-\eta \right)-q_{2}H-\frac{a_2\left(1-c\right)W_y}{{\left\{1+a_2\left(1-c\right)h{W}_x\right\}}^2},\\ a_{32}=\frac{\theta a_2\left(1-c\right)W_y}{{\left\{1+a_2\left(1-c\right)h{W}_x\right\}}^2},a_{33}=\theta \left[\frac{a_2\left(1-c\right)W_x}{1+a_2\left(1-c\right)h{W}_x}-\mu \right]-q_{3}H.\end{array}$

Accordingly, the linear stability analysis about the equilibrium points $E_i\mathrm{,}i\mathrm{<mo>\; =\; </mo>0,1,2,3}$ gives the following results:

1) The equilibrium point $E_0$ is unstable manifold in B-direction.

2) The equilibrium point $E_1$ is unstable manifold in W_x-direction.

3) The equilibrium point $E_2$ is unstable manifold in W_y-direction.

4) The characteristic equation corresponding to variational matrix $M\left(E_3\right)$ is given by ${\lambda }^3+A_1{\lambda }^2+A_2\lambda +A_3=0,$

where

$\begin{array}{l}{A}_1=-a_{11}-a_{22}-a_{33},\\ A_2=a_{22}{a}_{33}+a_{33}{a}_{11}+a_{11}{a}_{22}-a_{23}{a}_{32},\\ A_3=a_{23}{a}_{32}{a}_{11}+a_{12}{a}_{21}{a}_{33}-a_{11}{a}_{22}{a}_{33}.\end{array}$

Then by Routh-Hurwitz criteria equilibrium $E_3$ is locally asymptotically stable if $A_1>0,A_3>0$ and $A_1{A}_2>A_3$ and unstable if either of these conditions is not satisfied.

5. Persistence

Theorem 5.1: Let

$r_0>q_{1}H,\frac{\xi a_{1}K\left(r_0-q_{1}H\right)}{b_1{r}_0+K\left(r_0-q_{1}H\right)}>\xi \eta +q_{2}H$ and

$\mu <\min \left(\frac{\theta -q_{3}H}{\theta h},\frac{a_2\left(1-c\right)W_x}{1+a_2\left(1-c\right)h{W}_x}-\frac{q_{3}H}{\theta }\right)$ .

Then the system (2.1) is uniformly persistent.

Proof: Suppose u is a point in the positive octant and $O\left(u\right)$ is the orbit through u and $\omega$ is the omega limit set of the orbit through u Note that $\omega \left(u\right)$ is bounded.

We claim that $E_0\notin \omega \left(u\right)$ . If $E_0\in \omega \left(u\right)$, then by Butler-McGhee lemma there exists a point $u_0$ in $\omega \left(u\right)\cap M^S\left(E_0\right)$, where $M^S\left(E_0\right)$ denote the stable manifold of $E_0$ . Since $O\left(u_0\right)$ lies in $\omega \left(u\right)$, the condition $r_0>q_{1}H$ implies that $E_0$ is a saddle point and $M^S\left(E_0\right)$ is the $W_x-W_y$ plane, we conclude that $O\left(u_0\right)$ is unbounded, which is a contradiction. Again we claim that $E_1\notin \omega \left(u\right)$ . If $E_1\in \omega \left(u\right)$, then by Butler-McGhee lemma there exists a point $u_1$ in $\omega \left(u\right)\cap M^S\left(E_1\right)$, where $M^S\left(E_1\right)$ denote the stable manifold of $E_1$ .

Since $O\left(u_1\right)$ lies in $\omega \left(u\right)$, the condition $\frac{\xi a_{1}K\left(r_0-q_{1}H\right)}{b_1{r}_0+K\left(r_0-q_{1}H\right)}>\xi \eta +q_{2}H$, implies that $E_1$ is a saddle point and $M^S\left(E_1\right)$ is the $B-W_y$ plane, we conclude that $O\left(u_1\right)$ is unbounded, which is a contradiction.

Next, we show that $E_2\notin \omega \left(u\right)$ . If $E_2\in \omega \left(u\right)$, the condition $-\theta \mu +\frac{\theta a_2\left(1-c\right)W_x}{1+a_2\left(1-c\right)W_x}>q_{3}H$ is positive implies that $E_2$ is a saddle point. $M^s\left(E_2\right)$ is the $B-W_x$ plane, again we conclude that an unbounded orbit lies in $\omega \left(u\right)$, a contradiction.

Thus, $\omega \left(u\right)$ lies in the positive octant and the system (2.1) are persistent. Finally, since only the closed orbits and the equilibria form the omega limit set of the solutions on the boundary of $R_{+}^3$ and system (2.1) is deceptive, by Butler-McGhee lemma , implies that system (2.1) is uniformly persist.

6. Numerical Simulation for Bifurcation

In this section, we present numerical simulation to illustrate results obtained in previous sections. The system (2.1) is solved using fourth order Runge-Kutta Method with the help of MATLAB software package.

We shall vary c in system (2.1) so as to obtain a Hopf bifurcation. Now, we write autonomous system (2.1) in the form

$\stackrel{\dot{}}{n}=F\left(n,k\right)\text{,}$

where

$n=\left(B,W_x,W_y\right),k=\left(r_0,K,a_1,b_1,q_1,\xi ,H,\eta ,a_2,c,h,q_2,\theta ,\mu ,q_3\right)\text{.}$

We say that an ordered pair $\left(n_0,k_0\right)$ is a Hopf bifurcation point if

1) $F\left(n_0,k_0\right)=0$,

2) $J\left(n,k\right)$ has two complex conjugate eigenvalues ${\lambda }_{1,2}$ around $\left(n_0,k_0\right)$, ${\lambda }_{1,2}=a\left(n,k\right)\pm ib\left(n,k\right)$,

3) $a\left(n_0,k_0\right)=0,b\left(n_0,k_0\right)\ne 0$,

4) The third eigenvalue ${\lambda }_3\left(n_0,k_0\right)\ne 0$ .

Extensive numerical simulations are carried out for various values of parameters and for different sets of initial conditions. We take the parameters of the system (2.1) as,

$\begin{array}{l}{r}_0=3,K=50,a_1=1.0,b_1=0.5,q_1=0.5,\xi =1.0,\eta =0.8,a_2=1.0,\\ q_2=0.3,q_3=0.07,H=0.6,\theta =1.0,\mu =0.6\text{and}h=\mathrm{0.04.}\end{array}$ (6.1)

We consider the system

$\begin{array}{l}\frac{\text{d}B}{\text{d}t}=3\left(1-\frac{B}{50}\right)B-\frac{1\times B{W}_x}{0.5+B}-0.5\times 0.5\times B,\\ \frac{\text{d}{W}_x}{\text{d}t}=1\times W_x\left(\frac{1\times B}{0.5+B}-0.8\right)-\frac{1\times \left(1-c\right)W_x{W}_y}{1+1\times \left(1-c\right)\times 0.04\times W_x}-0.3\times 0.6\times W_x,\\ \frac{\text{d}{W}_y}{\text{d}t}=1\times W_y\left(\frac{1\times \left(1-c\right)W_x}{1+1\times \left(1-c\right)\times 0.04\times W_x}-0.6\right)-0.07\times 0.6\times W_y.\end{array}$ (6.2)

The system (6.2) has always non-negative equilibrium point $E_1\left(45,0,0\right)$ . The system (6.2) has positive equilibra $E_2\left(B_2,W_x{}_{2}0\right)$ and $E_3\left(B^{\ast },W_x^{\ast },W_y^{\ast }\right)$ if and if $c\in \left[0,1\right)$ .

1) Take $c=0.7$, $k=\left(3,50,1,0.5,0.5,1,0.8,0.6,1,0.7,0.04,0.3,1,0.6,0.07\right)$ .

The coordinates of $E_3$ and the corresponding eigenvalues are

$n_1=\left(44.18070,2.1964,0.030138\right),$

${\lambda }_1=0+0.0742354i,{\lambda }_2=0-0.0742354i,{\lambda }_3=-\mathrm{2.60203.}$

In this way ordered pair $\left(n_0,k_0\right)$ satisfied above all conditions (i-iv). So ordered pair $\left(n_0,k_0\right)$ is a Hopf point.

2) For $c=0.4<0.7$, $k=\left(3,50,1,0.5,0.5,1,0.8,0.6,1,0.4,0.04,0.3,1,0.6,0.07\right)$ .

The coordinates of $E_3$ and the corresponding eigenvalues are

$n_1=\left(44.59410,1.098201,0.01524\right),$

${\lambda }_1=0.00006+0.0746648i,{\lambda }_2=0.00006-0.0746648i,{\lambda }_3=-\mathrm{2.65146.}$

All eigenvalues have not negative real parts, only ${\lambda }_3$ has negative real part, so $E_3$ is always saddle point at $c=0.4$ .

3) But if we take $c=0.8>0.7$, $k=\left(3,50,1,0.5,0.5,1,0.8,0.6,1,0.8,0.04,0.3,1,0.6,0.07\right)$ .

The coordinates of $E_3$ and the corresponding eigenvalues are

$n_1=\left(43.759356,3.2946,0.044661\right),$

${\lambda }_1=-0.00005+0.073786i,{\lambda }_2=-0.00005-0.073786i,{\lambda }_3=-\mathrm{2.55164.}$

All eigenvalues have negative real parts, so equilibrium point is locally asymptotically stable at $c=0.8$ .

From Figures 2-6, we analyse that habitat complexity play important role in stability of ecosystem. Human-induced activities effect the habitat complexity. We also see that from Figure 3 and Figure 4, increasing the value of coefficient of human-induced activities, forestry biomass become extinct and wild life wholly depended on the forestry biomass also become extinct. It shows numerically that the Hopf point is found when $c=0.7$ . $E_3$ is unstable when $c<0.7$ and stable when $c>0.7$ . From this result, we can say that if we conserve habitat complexity for wildlife population, we can reduce the rate of depletion of forestry biomass. Human being depend on the forestry biomass but if they use forestry biomass keeping in the range of control of habitat complexity of wildlife population, we conserve the forestry biomass as well as wildlife population. The numerical study presented here shows that, using the parameter c as control, it is possible to break the stable behavior of the system (6.2) and drive it in an unstable state. Also, it is possible to keep the population levels in a required state using the above control.

7. Conclusion

In this paper, we have proposed a nonlinear mathematical model to study the

Figure 2. Here $B\left(0\right)=44.18070,W_x\left(0\right)=2.1964,W_y\left(0\right)=0.030138$ and other parameters are same as (6.1).

Figure 3. Here $B\left(0\right)=44.18070,W_x\left(0\right)=2.1964,W_y\left(0\right)=0.030138$ and other parameters are same as (6.1).

Figure 4. Here $B\left(0\right)=44.18070,W_x\left(0\right)=2.1964,W_y\left(0\right)=0.030138$ and other parameters are same as (6.1).

Figure 5. Here $B\left(0\right)=44.59410,W_x\left(0\right)=1.098201,W_y\left(0\right)=0.01524$ and other parameters are same as (6.1).

Figure 6. Here $B\left(0\right)=43.759356,W_x\left(0\right)=3.2946,W_y\left(0\right)=0.044661$ and other parameters are same as (6.1).

effect of human activities on forestry biomass and its effect on wildlife species living in the forest. The proposed model has been analyzed by the stability theory of differential equations. The conditions of existence of equilibrium points and their stability in have been obtained. We have also obtained conditions under which system persists by using differential inequality. Our numerical study shows that habitat complexity behaves as a control; it is possible that it breaks the stable behavior of the system (2.1) and drives it to an unstable state. Also, it is possible to keep the population levels in a required state using the above control. We see that effect of human-induced activities on forestry biomass and wild grazer population and the wild predator population. The effect of depletion of forestry biomass on the survival of wildlife species has also been shown graphically.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Rawat, G.S. and Ginwal, H.S. (2009) Conservation and Management of Forest Genetic Resources in India. Forest Genetic Resources Conservation and Management: Status in Seven South and Southeast Asian Countries, 21-46.
[2]	Dubey, B. and Narayanan, A.S. (2010) Modelling Effects of Industrialization, Population and Pollution on a Renewable Resource. Nonlinear Analysis: Real World Applications, 11, 2833-2848. https://doi.org/10.1016/j.nonrwa.2009.10.007
[3]	Agrawal, M., Fatima, T. and Freedman, H.I. (2010) Depletion of Forestry Resource Biomass Due to Industrialization Pressure: A Ratio-Dependent Mathematical Model. Journal of Biological Dynamics, 4, 381-396. https://doi.org/10.1080/17513750903326639
[4]	Agrawal, M. and Devi, S. (2012) A Resource-Dependent Competition Model: Effects of Population Pressure Augmented Industrialization. International Journal of Modeling, Simulation, and Scientific Computing, 3, Article ID: 1250003(1-22).
[5]	Chaudhary, M., Dhar, J. and Sahu, G.P. (2013) Mathematical Model of Depletion of Forestry Resource: Effect of Synthetic Based Industries. International Journal of Mathematical, Computational, Physical, Electrical and Computer Engineering, 7, 639-643.
[6]	Misra, A.K., Lata, K. and Shukla, J.B. (2014) A Mathematical Model for Depletion of Forestry Resources Due to Population and Population Pressure Augmented Industrialization. International Journal of Modeling, Simulation and Scientific Computing, 5, 1-16.
[7]	Agarwal, M. and Pathak, R. (2015) Conservation of Forestry Biomass with the Use of Alternative Resource. Open Journal of Ecology, 5, 87-109. https://doi.org/10.4236/oje.2015.54009
[8]	Shukla, J.B., Pal, V.N., Misra, O.P., Agarwal, M. and Shukla, A. (1988) Effects of Population and Industrialization on the Degradation of Biomass and Its Regeneration by Afforestation: A Mathematical Model. Journal of Biomathematics, 3, 1-9.
[9]	Shukla, J.B., Freedman, H.I., Pal, V.N., Misra, O.P., Agrawal, M. and Shukla, A. (1989) Degradation and Subsequent Regeneration of a Forestry Resource: A Mathematical Model. Ecological Modelling, 44, 219-229. https://doi.org/10.1016/0304-3800(89)90031-8
[10]	Agrawal, M. and Pathak, R. (2015) Conservation of Forestry Biomass and Wildlife Population: A Mathematical Model. Asian Journal of Mathematics and Computer Research, 4, 1-15.
[11]	Agrawal, M. and Pathak, R. (2011) Modeling the Effect of Harvesting of the Vegetation Biomass and Grazer Population on Predator Population with Habitat Complexity. International Journal of Mathematical Archive, 2, 2119-2134.
[12]	Pathak, R., Agarwal, M., Takeuchi, Y. and Shukla, J.B. (2017) Modeling the Depletion of Forest Resources Due to Population and Industrialization and Its Effect on Survival of Wildlife Species. International Journal of Mathematical Modelling and Analysis of Complex Systems (IJMMACS): Nature and Society, 3, 37-65.