Optimum Control for Spread of Pollutants through Forest Resources

Nita H. Shah^*, Moksha H. Satia, Bijal M. Yeolekar

Department of Mathematics, Gujarat University, Ahmedabad, India.

DOI: 10.4236/am.2017.85047                

Abstract

Pollution has become the most critical factor spread by forest resources through wood-based and non-wood based industries. In other words, pollution is omnipresent. In this paper, the major pollutants caused due to wood and non-wood based industries are discussed which are the primary resources of the forest in spreading the pollution. In order to study the impact of industrialization and associated pollution on forest resources, the system of non-linear ordinary differential equations is formulated. The controls are advised on both types of industries to reduce the pollution.

Keywords

Mathematical Model, Pollutants, Pollution, Forest Resources, Wood and Non-Wood Based Industries, System of Non-Linear Ordinary Differential Equation, Control

1. Introduction

The forest resources mean a large area covered by trees. It means the various types of vegetation automatically growing on forest land where forest is considered as if it grew trees in the past, or will grow trees in the future. Wood-based industries are the branch of production and employment based on the fabrication, processing and preparation of products from raw materials and merchandises of wood and wood-based pulp products. Non-wood based industries are the branch of production and employment based on fabrication, processing and preparation of products from water, energy, chemicals etc. The term “pollution” is a substance into the environment which has harmful or poisonous effects on living beings. “Pollutants” are the components of pollution. It is observed that most of the pollution is associated with man-made industries. Wood and non- wood based industries affect the environment through pollutants emitted from them merged in air as well as absorbed by forest region can be harmful for nature and can kill human organisms, essential microbes etc.

The forest resources are significant for human and some organisms. But after industrial revolution in 18^th century, industries are also growing very speedy [1] . This growth of wood and non-wood based industries has reduced the density of forest region. Damodar Valley, Nowamundi, Saranda are example of reduced forest resources [2] . Once Damodar Valley was covering 65% of forest area, now-a-days it is surrounded by only 0.05% [3] . In past, few years, the temperature of the environment is increasing due to the emission of pollutants which is examined by scientists and ecologists. This gives opposite impact on humans and environment [4] [5] . Absorption of pollutants by the plants is harmful and which affected the growth of forest resources [6] [7] [8] [9] .

This motivated to formulate the system in which the effect of industrialization on the forest resources is analyzed. Some researchers have studied the mathematical model for the effects of industrialization and pollution on forest resources. [10] studied the models for the effect of toxicant in single-species and predator-prey system. [11] analysed the modeling effect of an intermediate toxic product formed by uptake of a toxicant on plant biomass. [12] introduced the effects of industrialization and pollution on resource biomass with the help of a mathematical model. [13] have performed the modeling effects of industrialization, population and pollution on renewable resources. [14] deliberated modeling effects of primary and secondary toxicants on renewable resources.

In this paper, a mathematical model is formulated with hereditary transmission of SIRS model in Section 2. The stability analysis of the transmission model is derived in Section 3. Sensitivity analysis is carried out in Section 4. Optimal control for the forest resources is discussed in Section 5. In Section 6, the model validated with numerical simulation and analysis.

2. Mathematical Modeling

In the society, there are different types of industries and pollution. Everybody in the society plays a role to decrease the pollution. Therefore, in the proposed model, five discrete compartments viz. the density of forest resources (F), the density of wood based industries (W), the density of non-wood based industries (I), the pollutants through wood based industries (P_W) and the pollutants through non-wood based industries (P_I) are considered. u₁ is the rate which decreases wood based industries to control the usage of forest resources. u₂ and u₃ are the control rates which decreases pollutants due to wood and non-wood based industries, respectively.

The notations and parametric values for the dynamical model are exhibited in Table 1.

Using these notations and assumptions which are required for formulating the mathematical model, the transmission diagram of forest resources is shown in Figure 1.

The dynamics of forest resources transmission in wood and non-wood based

Table 1. Notation and parametric values.

Figure 1. Forest resources transmission diagram.

industries with associated pollutants is described as follows:

$\frac{\text{d}F}{\text{d}t}=B-\left(\beta W+g\right)F-{\beta }_{1}FW+{\epsilon }_1{P}_W-{\gamma }_{1}F-{\gamma }_{2}F+{\epsilon }_2{P}_I+u_{1}W-\mu F$ (1)

$\frac{\text{d}W}{\text{d}t}=\left(\beta W+g\right)F+{\beta }_{1}FW-{\delta }_{1}W+{\delta }_{2}I-{\eta }_{1}W-u_{1}W-u_{2}W-\mu W$ (2)

$\frac{\text{d}I}{\text{d}t}=QI+{\delta }_{1}W-{\delta }_{2}I-{\eta }_{2}I-u_{3}I-\mu I$ (3)

$\frac{\text{d}{P}_W}{\text{d}t}={\eta }_{1}W-{\epsilon }_1{P}_W+{\gamma }_{1}F+u_{2}W-{\mu }_W{P}_W$ (4)

$\frac{\text{d}{P}_I}{\text{d}t}={\eta }_{2}I-{\epsilon }_2{P}_I+{\gamma }_{2}F+u_{3}I-{\mu }_I{P}_I$ (5)

Equations (1) to (5) is described as system (1) in the model.

With $F+W+I+P_W+P_I=N$ and $F>0;W,I\ge 0;P_w\ge 0;P_I\ge 0$

Adding all the above system of differential equations gives,

#Math_26# (6)

This gives,

$\underset{t\to \infty }{\lim }\sup \left(F+W+I+P_W+P_I\right)\le \frac{B}{\mu }$ (7)

Therefore, the feasible region for system (1) is

$\Lambda =\left\{\left(F+W+I+P_W+P_I\right)/F+W+I+P_W+P_I\le \frac{B}{\mu },F>0;W,I,P_W,P_I\ge 0\right\}$ . (8)

Thus, the equilibrium state of the system (1) is $X_0=\left(\frac{B}{\mu },0,0,0,0\right)$

Next, the basic reproduction number $R_0$ can be calculated using the next generation matrix.

Let $X^{\prime }={\left(W,F,I,P_W,P_I\right)}^{\prime }$ , where dash denotes derivative. So,

$X^{\prime }=\frac{\text{d}X}{\text{d}t}=\mathbb{F}\left(X\right)-V\left(X\right)$ (9)

where $\mathbb{F}\left(X\right)$ denotes the rate of appearance of new individual in compartment and $V\left(X\right)$ represents the rate of transfer of culture, which is given by

$\begin{array}{l}\mathbb{F}\left(X\right)=\left[\begin{array}{c}\left(\beta +{\beta }_1\right)FW\\ 0\\ 0\\ 0\\ 0\end{array}\right],\\ V\left(X\right)=\left[\begin{array}{c}-gF+{\delta }_{1}W-{\delta }_{2}I+{\eta }_{1}W+u_{1}W+u_{2}W+\mu W\\ -B+\beta WF+gF+{\beta }_{1}WF-{\epsilon }_1{P}_W+{\gamma }_{1}F+{\gamma }_{2}F-{\epsilon }_2{P}_I-u_{1}W+\mu F\\ -QI-{\delta }_{1}W+{\delta }_{2}I+{\eta }_{2}I+u_{3}I+\mu I\\ -{\eta }_{1}W+{\epsilon }_1{P}_W-{\gamma }_{1}F-u_{2}W+{\mu }_W{P}_W\\ -{\eta }_{2}I+{\epsilon }_2{P}_I-{\gamma }_{2}F+u_{3}I+{\mu }_I{P}_I\end{array}\right]\end{array}$

Now,

$D\mathbb{F}\left(X_0\right)=\left[\begin{array}{cc}f& 0\\ 0& 0\end{array}\right],DV\left(X_0\right)=\left[\begin{array}{cc}v& 0\\ J_1& J_2\end{array}\right]$

where f and $v$ are $5\times 5$ matrices defined as

$f=\left[\frac{\partial {\mathbb{F}}_i\left(X_0\right)}{\partial X_j}\right],v=\left[\frac{\partial V_i\left(X_0\right)}{\partial X_j}\right]$

$\begin{array}{l}f=\left[\begin{array}{ccccc}\frac{\left(\beta +{\beta }_1\right)B}{\mu }& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right],\\ v=\left[\begin{array}{ccccc}{\delta }_1+{\eta }_1+u_1+u_2+\mu & -g& -{\delta }_2& 0& 0\\ \frac{\left(\beta +{\beta }_1\right)B}{\mu }-u_1& g+{\gamma }_1+{\gamma }_2+\mu & 0& -{\epsilon }_1& -{\epsilon }_2\\ -{\delta }_1& 0& -Q+{\delta }_2+{\eta }_2+u_3+\mu & 0& 0\\ -{\eta }_1-u_2& -{\gamma }_1& 0& {\epsilon }_1+{\mu }_W& 0\\ 0& -{\gamma }_2& -{\eta }_2-u_3& 0& {\epsilon }_2+{\mu }_I\end{array}\right]\end{array}$

Here, $v$ is non-singular matrix, so the basic reproduction number $R_0$ is

$R_0=\text{spectralradiusofmatrix}f{v}^{-1}.$ (10)

$R_0=\left(\frac{\left(\beta +{\beta }_1\right)B{C}_1\left(B_1{D}_1{D}_2-D_2{\gamma }_1{\epsilon }_1-D_1{\gamma }_2{\epsilon }_2\right)}{\mu \left[\left(A_1{C}_1-{\delta }_1{\delta }_2\right)\left(B_1{D}_1{D}_2-D_2{\gamma }_1{\epsilon }_1-D_1{\gamma }_2{\epsilon }_2\right)+g\left(A_2{C}_1{D}_1{D}_2+C_2{D}_1{\delta }_1{\epsilon }_2+A_3{C}_1{D}_2{\epsilon }_1\right)\right]}\right)$ (11)

where

$A_1={\delta }_1+{\eta }_1+u_1+u_2+\mu ,A_2=\frac{\left(\beta +{\beta }_1\right)B}{\mu }-u_1,A_3=-{\eta }_1-u_2,B_1=g+{\gamma }_1+{\gamma }_2+\mu ,$

$C_1=-Q+{\delta }_2+{\eta }_2+u_3+\mu ,C_2=-{\eta }_2-u_3,D_1={\epsilon }_1+{\mu }_w,D_2={\epsilon }_2+{\mu }_I$

In next section, equilibrium of the forest resources transmission model is discussed.

3. Equilibrium

The equilibrium for the local and global stability of the forest transmission model are discussed here.

3.1. Local Stability

The forest resources equilibrium is locally asymptotically stable if all the eigenvalues of the matrix have positive real values [15] . The Jacobian matrix for sys-

tem (1) at $X_0=\left(\frac{B}{\mu },0,0,0,0\right)$ given by

$J=\left[\begin{array}{ccccc}-g-{\gamma }_1-{\gamma }_2-\mu & -\left(\beta +{\beta }_1\right)\frac{B}{\mu }+u_1& 0& {\epsilon }_1& {\epsilon }_2\\ g& \left(\beta +{\beta }_1\right)\frac{B}{\mu }-{\delta }_1-{\eta }_1-u_1-u_2-\mu & {\delta }_2& 0& 0\\ 0& {\delta }_1& Q-{\delta }_2-{\eta }_2-u_3-\mu & 0& 0\\ {\gamma }_1& {\eta }_1+u_2& 0& -{\epsilon }_1-{\mu }_W& 0\\ {\gamma }_2& 0& {\eta }_2+u_3& 0& -{\epsilon }_2-{\mu }_I\end{array}\right]$

Using the parametric values given in the Table 1,

$\begin{array}{l}trace\left(J\right)=-g-\left(\beta +{\beta }_1\right)W-{\gamma }_1-{\gamma }_2+\left(\beta +{\beta }_1\right)F-{\delta }_1-{\eta }_1+Q-{\delta }_2\\ \text{}\text{}\text{}-{\eta }_2-{\epsilon }_1-{\epsilon }_2-u_1-u_2-u_3-3\mu -{\mu }_W-{\mu }_I=-6.97<0\end{array}$ (12)

Hence, system (1) is locally stable.

3.2. Global Stability

The forest resources transmission model is globally stable is $\det \left(I-f{v}^{-1}\right)>0$ .

$\det \left(I-f{v}^{-1}\right)=1-R_0=1-\text{0}\text{.4960}=0.5040>0$ (13)

Therefore, system (1) is also globally stable.

4. Sensitivity Analysis

In this section, the sensitivity analysis for all parameters are discussed in Table 2.

The normalised sensitivity index of the parameters is computed by using the

following formula: ${\Upsilon }_{\alpha }^{R_0}=\frac{\partial R_0}{\partial \alpha }\cdot \frac{\alpha }{R_0}$ where $\alpha$ denotes the model parameter.

The rate of compactness degree of forest resources, the constant rate of resources, migration of wood based industries to forest region, the depletion rate of forest resources due to wood based industries, the growth rate of wood based industries due to forest resources, the rate of competitive effects of I on W, the loss of pollutants generated by wood based industries due to forest resources, the loss of pollutants generated by non-wood based industries due to forest resources and the growth rate of pollutants generated by wood based industries

Table 2. Sensitivity analysis.

have positive effect on $R_0$ which means they are helping us to save forest resources. Other parameters have negative impact on model.

5. Optimal Control

The objective of the model is to minimize the number of pollutants through wood and non-wood based industries to revive forest resources. The control functions are united to achieve the objective. The objective function for the mathematical model of forest resources in system (1) along with the optimal control is given by

$J\left(u_i,\Omega \right)={\displaystyle \underset{0}{\overset{T}{\int }}\left(A_1{F}^2+A_2{W}^2+A_3{I}^2+A_4{P}_W^2+A_5{P}_I^2+w_1{u}_1^2+w_2{u}_2^2+w_3{u}_3^2\right)\text{d}t}$ (14)

where, $\Omega$ denotes set of all compartmental variables, $A_1,A_2,A_3,A_4,A_5$ denote non-negative weight constants for $F,W,I,P_W,P_I$ compartments respectively and $w_1,w_2,w_3$ are weight constants for control variables $u_1,u_2,u_3$ respectively.

As, the weight parameters $w_1,w_2$ and $w_3$ are constants of forest resources control $\left(u_1\right)$ , wood based industries control $\left(u_2\right)$ and non-wood based industries $\left(u_3\right)$ , from which the optimal control condition is normalized. $u_1$ is the control variable for minimizing the use of forest resources. $u_2$ and $u_3$ are the control rates which minimize the density of wood and non-wood based industries respectively which automatically reduce the pollutants also. To compute the values of control variables $u_1,u_2$ and $u_3$ from $t=0$ to $t=T$ such that

$J\left(u_1\left(t\right),u_2\left(t\right),u_3\left(t\right)\right)=\min \left\{J\left(u_i{}^{*},\Omega \right)/\left(u_1,u_2,u_3\right)\in \varphi \right\}$ (15)

where $\varphi$ is a smooth function on the interval $\left[0,1\right]$ . The optimal controls denoted by $u_i{}^{*},i=1,2,3$ are found by accumulating all the integrands of Equation (14) using the lower bounds and upper bounds respectively with the results of [16] .

Now, using the pontrygin’s principle from [17] , to minimize the cost function in (14) by constructing Lagrangian function consisting of state equations and adjoint variables $A_i=\left({\lambda }_1,{\lambda }_2,{\lambda }_3,{\lambda }_4,{\lambda }_5\right)$ as

$\begin{array}{l}L\left(\Omega ,A_i\right)=A_1{F}^2+A_2{W}^2+A_3{I}^2+A_4{P}_W^2+A_5{P}_I^2+w_1{u}_1^2+w_2{u}_2^2+w_3{u}_3^2\\ +{\lambda }_1\left(B-\left(\beta W+g\right)F-{\beta }_{1}FW+{\epsilon }_1{P}_W-{\gamma }_{1}F-{\gamma }_{2}F+{\epsilon }_2{P}_I+u_{1}W-\mu F\right)\\ +{\lambda }_2\left(\left(\beta W+g\right)F+{\beta }_{1}FW-{\delta }_{1}W+{\delta }_{2}I-{\eta }_{1}W-u_{1}W-u_{2}W-\mu W\right)\\ +{\lambda }_3\left(QI+{\delta }_{1}W-{\delta }_{2}I-{\eta }_{2}I-u_{3}I-\mu I\right)\\ +{\lambda }_4\left({\eta }_{1}W-{\epsilon }_1{P}_W+{\gamma }_{1}F+u_{2}W-{\mu }_W{P}_W\right)\\ +{\lambda }_5\left({\eta }_{2}I-{\epsilon }_2{P}_I+{\gamma }_{2}F-u_{3}I-{\mu }_I{P}_I\right)\end{array}$ (16)

The partial derivative of the Lagrangian function with respect to each variable of the compartment gives the adjoint equation variables $A_i=\left({\lambda }_1,{\lambda }_2,{\lambda }_3,{\lambda }_4,{\lambda }_5\right)$ corresponding to the system (1) which is as follows:

$\begin{array}{l}{\stackrel{\dot{}}{\lambda }}_1=-\frac{\partial L}{\partial F}=-2W_{1}F+\left({\lambda }_1-{\lambda }_2\right)\left(\beta W+g\right)+\left({\lambda }_1-{\lambda }_2\right){\beta }_{1}W\\ +\left({\lambda }_1-{\lambda }_4\right){\gamma }_1+\left({\lambda }_1-{\lambda }_5\right){\gamma }_2+\mu {\lambda }_1\end{array}$ (17)

$\begin{array}{l}{\stackrel{\dot{}}{\lambda }}_2=-\frac{\partial L}{\partial W}=-2W_{2}W+\left({\lambda }_1-{\lambda }_2\right)\beta F+\left({\lambda }_1-{\lambda }_2\right){\beta }_{1}F+\left({\lambda }_2-{\lambda }_1\right)u_1\\ +\left({\lambda }_2-{\lambda }_3\right){\delta }_1+\left({\lambda }_2-{\lambda }_4\right){\eta }_1+\left({\lambda }_2-{\lambda }_4\right)u_2+\mu {\lambda }_2\end{array}$ (18)

$\begin{array}{c}{\stackrel{\dot{}}{\lambda }}_3=-\frac{\partial L}{\partial I}\\ =-2W_{3}I+\left({\lambda }_3-{\lambda }_2\right){\delta }_2+\left({\lambda }_3-{\lambda }_5\right)u_3+\left({\lambda }_3-{\lambda }_5\right){\eta }_2+Q{\lambda }_3+\mu {\lambda }_3\end{array}$ (19)

${\stackrel{\dot{}}{\lambda }}_4=-\frac{\partial L}{\partial P_W}=-2W_4{P}_W+\left({\lambda }_4-{\lambda }_1\right){\epsilon }_1+{\mu }_W{\lambda }_4$ (20)

${\stackrel{\dot{}}{\lambda }}_5=-\frac{\partial L}{\partial P_W}=-2W_5{P}_I+\left({\lambda }_5-{\lambda }_1\right){\epsilon }_2+{\mu }_I{\lambda }_5$ (21)

The necessary condition for Lagrangian function L to be optimal for controls are

${\stackrel{\dot{}}{u}}_1=-\frac{\partial L}{\partial u_1}=-2W_6{u}_1+\left({\lambda }_2-{\lambda }_1\right)W=0$ (22)

${\stackrel{\dot{}}{u}}_2=-\frac{\partial L}{\partial u_2}=-2W_7{u}_2+\left({\lambda }_2-{\lambda }_4\right)W=0$ (23)

${\stackrel{\dot{}}{u}}_3=-\frac{\partial L}{\partial u_3}=-2W_6{u}_3+\left({\lambda }_3-{\lambda }_5\right)I=0$ (24)

To find the values of $u_1,u_2$ and $u_3$ solve Equations (22), (23) and (24) then

$u_1=\frac{\left({\lambda }_2-{\lambda }_1\right)W}{2W_6}$ , $u_2=\frac{\left({\lambda }_2-{\lambda }_4\right)W}{2W_7}$ and $u_3=\frac{\left({\lambda }_3-{\lambda }_5\right)I}{2W_6}$ (25)

Thus, the required optimal control condition is computed as

$u_1^{\ast }=\max \left(a_1,\min \left(b_1,\frac{\left({\lambda }_2-{\lambda }_1\right)W}{2W_6}\right)\right)$ (26)

$u_2^{\ast }=\max \left(a_2,\min \left(b_2,\frac{\left({\lambda }_2-{\lambda }_4\right)W}{2W_7}\right)\right)$ (27)

$u_3^{*}=\max \left(a_3,\min \left(b_3,\frac{\left({\lambda }_3-{\lambda }_5\right)I}{2W_6}\right)\right)$ (28)

In next section the optimal control is calculated numerically to support the analytical results.

6. Numerical Simulation

Using the data given in Table 1 and Table 2, the sensitivity on model parameters is carried out.

The Figure 2 specifies that as the depletion rate of forest resources due to wood based industries increases the density of forest resources decreases.

From Figure 3, it is observed that increase in growth rate of wood based industries due to forest resources the density of forest resources decreases.

The Figure 4 indicates that with the loss of pollutants generated by non-wood based industries the forest resources increases.

Figure 2. Effects of the depletion rate of forest resources due to wood based industries on forest resources.

Figure 3. Effects of the growth rate of wood based industries due to forest resources on forest resources.

From Figure 5 it is concluded that as the constant rate of resources provided to non-wood based industries increases the forest resources decreases.

From Figure 6, one can see that forest should be controlled 12% in 50 years, wood based industries should be controlled 33% in 68 years and non-wood

Figure 4. Effects of the loss of pollutants generated by non-wood based industries on forest resources.

Figure 5. Effect of the constant rate of resources provided to non-wood based industries on forest resources.

based industries should be controlled 75% in 80 years. Even more, the figure shows that after applying this control on non-wood based industries for 163 years, the forest resources will revive for next 37 years.

With control forest resources degradation reduces at a lower rate compare to no effective majors are taken up as shown in Figure 7.

Figure 8 suggest that wood based industries can be controlled with effective majors at a lower rate compare to no control over it and when control is applied to wood based industries it decreases by 7%. Similar observation is from Figure 9 for non-wood based industries. In fact, it decreases by 4%.

Figure 6. Control variables.

Figure 7. Forest resources with control and without control.

Figure 8. Wood based industries with control and without control.

Figure 9. Non-wood based industries with control and without control.

7. Conclusions

In this paper, a mathematical model is formulated to study the spread of the forest transmission with wood and non-wood based industries. An optimal control for spread of the pollutants through the forest resources to study the effects of Wood, Non-wood based industries and the pollution emitted through them on the density of forest resources. The wood-based industries reduce the density of forest resources directly by harvesting as well as indirectly by pollutants. But non-wood based industries reduce the density of forest resources only indirectly by pollutants. Therefore, by more industries the forest resources are affected and may be wiped out.

The stability of forest resources model discussed with numerical data. The basic reproduction number is computed as 0.4960, which shows that controls on construction of wood and non-wood based industries will be beneficial to reduce the pollution. This suggested growing more and more forest resources and putting up less number of industries per human usage.

Acknowledgements

Authors sincerely thank for the constructive comments of the reviewers. The authors thank DST-FIST file # MSI-097 for technical support to the department.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Dubey, B., Sharma, S., Sinha, P. and Shukla, J.B. (2009) Modeling the Depletion of Forestry Resources by Population and Population Pressure Augmented Industrialization. Applied Mathematical Modeling, 33, 3002-3014.
[2]	Lata, K., Dubey, B. and Misra, A.K. (2016) Modeling the Effects of Wood and Non-Wood Based Industries on Forestry Resources. Natural Resource Modeling, 29, 559-580. https://doi.org/10.1111/nrm.12111
[3]	Priyadarshi, N. (2010) Effects of Mining on Environment in the State of Jharkhand, India-Mining Has Caused Severe Damage to the Land Resources.
[4]	Intergovernmental Panel on Climate Change (2014) Climate Change 2014-Impacts, Adaptation and Vulnerability: Regional Aspects. Cambridge University Press, Cambridge, 1-32.
[5]	Patwardhan, A., Semenov, S., Schnieder, S., Burton, I., Magadza, C., Oppenheimer, M. and Sukumar, R. (2007) Assessing Key Vulnerabilities and the Risk from Climate Change. In: Climate Change, 779-810.
[6]	Chappelka, A.H. and Samuelson, L.J. (1998) Ambient Ozone Effects on Forest Trees of the Eastern United States: A Review. New Phytologist, 139, 91-108. https://doi.org/10.1046/j.1469-8137.1998.00166.x
[7]	Davison, A.W. and Barnes, J.D. (1998) Effects of Ozone on Wild Plants. New Phytologist, 139, 135-151. https://doi.org/10.1046/j.1469-8137.1998.00177.x
[8]	Treshow, M. (1968) Impact of Air Pollutants on Plant Populations. Phytopathology 58, 8.
[9]	Treshow, M. (1984) Air Pollution and Plant Life.
[10]	Freedman, H.I. and Shukla, J.B. (1991) Models for the Effect of Toxicant in Single-Species and Predator-Prey Systems. Journal of Mathematical Biology, 30, 15-30. https://doi.org/10.1007/BF00168004
[11]	Naresh, R., Sundar, S. and Shukla, J.B. (2006) Modeling the Effect of an Intermediate Toxic Product Formed by Uptake of a Toxicant on Plant Biomass. Applied Mathematics and Computation, 182, 151-160.
[12]	Dubey, B., Upadhyay, R.K. and Hussain, J. (2003) Effects of Industrialization and Pollution on Resource Biomass: A Mathematical Model. Ecological Modeling, 167, 83-95.
[13]	Dubey, B. and Narayanan, A.S. (2010) Modeling Effects of Industrialization, Population and Pollution on a Renewable Resource. Nonlinear Analysis: Real World Applications, 11, 2833-2848.
[14]	Shukla, J.B., Agrawal, A.K., Sinha, P. and Dubey, B. (2003) Modeling Effects of Primary and Secondary Toxicants on Renewable Resources. Natural Resource Modeling, 16, 99-120.
[15]	Al-Amoudi, R., Al-Sheikh, S. and Al-Tuwairqi, S. (2014) Qualitative Behavior of Solutions to a Mathematical Model of Memes Transmission. International Journal of Applied Mathematical Research, 3, 36-44.
[16]	Fleming, W.H. and Rishel, R.W. (2012) Deterministic and Stochastic Optimal Control. Vol. 1, Springer Science & Business Media, New York.
[17]	Pontriagin, L.S., Boltyanskii, V.G., Gamkrelidze, R.V. and Mishchenko, E.F. (1986) The Mathematical Theory of Optimal Process. Gordon and Breach Science Publishers, New York, 4-5.