Applied Science and Innovative Research ISSN 2474-4972 (Print) ISSN 2474-4980 (Online) Vol. 1, No. 1, 2017 www.scholink.org/ojs/index.php/asir 1 Harvesting on Facultative Mutualist Prey Species in Presence of a Predator Saroj Kumar Chattopadhyay1* 1 Principal, Chandraketugarh Sahidullah Smriti Mahavidyalaya Berachampa, North 24-Parganas, West Bengal, India * Saroj Kumar Chattopadhyay, E-mail: saroj.scc@gmail.com Received: January 15, 2017 Accepted: January 24, 2017 Online Published: February 5, 2017 doi:10.22158/asir.v1n1p1 URL: http://dx.doi.org/10.22158/asir.v1n1p1 Abstract This paper proposes a model with two preys of facultative mutualist type and one predator. Linear predation functions are considered and preys are only considered to be harvested. The stability of the model is analyzed theoretically and numerically in this paper. The optimal harvest policy is studied and the solution is derived in the interior equilibrium case using Pontryagin’s maximum principle. Finally, some numerical simulations are discussed. Keywords mutualism, facultative, obligate, stability, optimal equilibrium 1. Introduction Mutualism is an interaction in which species help one another. Janzen (1985) has argued that most of the mutualisms can be classified into one of the four classes: seed-dispersal mutualism, pollination mutualisms, digestive mutualisms, and protective mutualisms. In this paper we are interested to discuss the protective mutualisms of fish species. This type of mutualisms may be facultative mutualism or obligate mutualism. In facultative mutualism the interaction between the species is helpful but not essential but in obligate mutualism neither mutualist can survive without the other. Many fish species form protective mutualisms. A particularly well-known example involves tropical anemone fishes, or clown fishes, and their anemones. Clown fishes are immune to stinging nematocysts of giant sea anemones and will and nest amongst their tentacles. Horse mackerels appear to have a similar relationship with Portuguese man-of war jellyfish. In recent past many works on mutualism have done (Lengeler et al., 1999; Wallin, 1923, 1927; Margulis, 1970, 1981). Those are not harvesting models. Simultaneously, some extraordinary harvesting models are studied by Clark (1985, 1990) and some other ecologists and scientists (Mesteron-Gibbons, 1998; Kot, 2001; Kar & Chaudhuri, 2004; Strobele & Wacker, 1995; Dai & Tang, www.scholink.org/ojs/index.php/asir Applied Science and Innovative Research Vol. 1, No. 1, 2017 2 Published by SCHOLINK INC. 1998) also studied the management and behaviour dynamic harvesting models. Mark Kot (2001) discussed a two species protective mutualism model. In which there are two species with population sizes 1N and 2N , and each species grows logistically in the absence of other. The model is )1( 1 2121 11 1 k NN Nr dt dN   (1) )1( 2 1212 22 2 k NN Nr dt dN   (2) where 12 and 21 are the measures of the strength of positive effect of species 2 on the species 1 and of species 1 on species 2. This model is for facultative mutualism so far 2121 ,,, kkrr are all positive. That is each species can, in other words, survive without its mutualist. The species may surpass their carrying capacity or may undergo unlimited growth what has been called “an orgy of mutual benefaction” (May, 1981) depending upon the values of the strength of one species on another. The equations (1) and (2) may be used for obligate mutualism if we take 2121 ,,, kkrr all negative. In this case neither species can survive on its own; each species is banking on the other to save it. The facultative models are generally more stable than obligate models. We state the definition of facultative or co-operative models. Definition 1.1: The system ),( 21 1 NNf dt dN  ),( 21 2 NNg dt dN  defined on D 2R is co-operative if 0,0 12       N g N f for all DNN ),( 21 . In this paper, we study the problem of harvesting two facultative species in the presence of a predator species which feeds on both the facultative prey species. The predator species is not harvested. The problem is clearly stated in the section 2. We have examined the equilibrium of the system and the conditions of their existence in section 3. The local stability of the steady state solutions is examined in section 4. We derive an optimal harvesting policy in section 5. Numerical examples are discussed in section6. The paper ends with a brief conclusion in section 7. 2. Formulation of the Model The governing equations of our model are, 111311 1 2121 11 1 )1( NEqNNa K NN Nr dt dN     (3) 222322 2 1212 22 2 )1( NEqNNa K NN Nr dt dN     (4) www.scholink.org/ojs/index.php/asir Applied Science and Innovative Research Vol. 1, No. 1, 2017 3 Published by SCHOLINK INC. 332223111 3 dNNNaNNa dt dN   (5) where 21 , NN are population sizes of the prey species and 3N the population size of the predator at any time t. Here, )0(),0( 21  rr are intrinsic growth rate of the first two species. Since we are not making a case study in respect of a specific prey-predator community, we have opted logistic growth rate. The parameters )0(),0( 21  KK are carrying capacities of first two species; 21 , aa , both positive, are predation rates on which the third species feeds on the first two species respectively; 21 , EE are the harvesting efforts given on the first two species, the third species is not harvested; 21, qq (both positive) are catchability co-efficients of 21,NN respectively. The catch rate functions 111 NEq and 222 NEq are based on CPUE (catch per unit effort) hypothesis. The parameters )0(),0( 21   are known as conversion factors and d is the mortality rate of predator species. Here 12 measures the strength of the positive effect of the species 2N on the species 1N and 21 measures the strength of the positive effect of the species 1N on the species 2N . Considering the definition (1.1) we see that in absence of predator 02 1 121 2    N K r N f  , and 01 2 212 1    N K r N g  . Thus our model is co-operative. We hypothetically consider that the predator does not disturb the facultative mutualism of prey species. 3. The Equilibria of the Model and the Existence Conditions The biological equilibria or the steady state solutions are obtained by solving 0 . 1 N , 0 . 2 N , 0 . 3 N where . iN (i=1,2,3) are time derivatives of iN (i=1,2,3) respectively. Solving these equations we get the points 6543210 ,,,,,, PPPPPPP of equilibrium. The co-ordinates of these points and corresponding conditions of existence are given bellow. The point )0,0,0(0P is trivial which always exists, the point ))1(,,0( 2 22 2222 2 22 1 a Eq aK d a r a d P   exists if 222 2 222 aK dr Eqr   , (6) )0),(,0( 222 2 2 2 Eqr r K P  exists if 222 Eqr  >0, (7) )0,0),(( 111 1 1 3 Eqr r K P  exists if 111 Eqr  >0, (8) ))1(,0,( 1 11 1111 1 11 4 a Eq aK d a r a d P   exists if 111 Eqr  > 111 1 aK dr  , (9) )0,,( 215 NNP where )]()([ 1 1 22212 2 2 111 1 1 2112 1 Eqr r K Eqr r K N      , )]()([ 1 1 222 2 2 11121 1 1 2112 2 Eqr r K Eqr r K N      , and it exists if www.scholink.org/ojs/index.php/asir Applied Science and Innovative Research Vol. 1, No. 1, 2017 4 Published by SCHOLINK INC. 111 Eqr  > 0 , 222 Eqr  > 0 and 12112  . (10) The other equilibrium point is interior equilibrium and is ),,( * 3 * 2 * 16 NNNP where 2 11 2 22 22 1 1 212 2 2 121 1 1 21 1 2 2 2222221212 1 1 21121 2 2 * 1 )( )()( a K r a K r K r K r aa da K r Eqraada K r Eqra N      (11) 2 11 2 22 22 1 1 212 2 2 121 1 1 21 2 1 1 1111121121 2 2 12212 2 1 * 2 )( )()( a K r a K r K r K r aa da K r Eqraada K r Eqra N      (12) 2 11 2 22 22 1 1 212 2 2 121 1 1 21 2112 21 21 1112 1 1 22 1 1 22211 2 2 2221 2 2 111 * 3 )( )1())(())(( a K r a K r K r K r aa KK rdr a K r a K r Eqra K r a K r Eqr N      (13) The interior equilibrium point exists if                )1())(())(( )()( )()( ,, 2112 21 21 1 1 111222222 2 2 112221111 1111212 1 1 121 2 2 1 2 1222 2222211 2 2 212 1 1 2 2 2111 222111    KK rdr K r aaEqr K r aaEqr Eqraada K r da K r aEqr Eqraada K r da K r aEqr EqrEqr (14) 4. Local Stability Analysis To analyze the local stability at the different equilibria we consider the following community matrix and use variational principle.                       dNaNaNaNa Nav K Nr Na K Nr v J 222111322311 2222 2 2212 11 1 1121 11    (15) where 1 11 1131 1 2121 111 1 K Nr EqNa K NN rv          (16) 2 22 2232 2 1212 222 1 K Nr EqNa K NN rv          (17) At )0,0,0(0P the community matrix reduces to               d Eqr Eqr J 00 00 00 222 111 0 www.scholink.org/ojs/index.php/asir Applied Science and Innovative Research Vol. 1, No. 1, 2017 5 Published by SCHOLINK INC. whose eigenvalues are 111 1 0 Eqr  , 222 2 0 Eqr  , d3 0 . Thus origin is a stable node if , 1 1 1 q r E  2 2 2 q r E  that is if the harvesting efforts are more than the corresponding biotechnical productivity of the two facultative prey species. We also see that if the system be stabilized at origin then no other equilibrium of the system will be found. At ))1(,,0( 2 22 2222 2 22 1 a Eq aK d a r a d P   the Jacobian matrix reduces to                      0 00)1( 122111 2222 2 222 212 1111 221 12 1 1 zaza d aK dr aK dr zaEq aK d r J      Where 2 22 2222 2 1 )1( a Eq aK d a r z   . One of the eigenvalues of 1J is 1111 221 12 1 1 1 )1( zaEq aK d r     which may be positive or negative depending upon the values of the parameters. In the inspection of other two eigenvalues we see that they are the roots of the quadratic 012 222 22  zda aK dr    . The sum of whose roots= 0 222 2  aK dr  , the product of the roots= 012  zda with the assumption that 1P exists. Thus other two eigenvalues are one positive and one negative. Hence 1P is an unstable equilibrium, whatever the sign of 1111 221 12 1 1 1 )1( zaEq aK d r     may be. At )0),(,0( 222 2 2 2 Eqr r K P  the community matrix becomes                       dEqr r Ka Eqr r Ka EqrEqr EqEqr rK K r J )(00 )()()( 00))(1( 222 2 222 222 2 22 22222221 11222 21 122 1 2    . Whose eigenvalues are, )()( 22212 12 21 111 1 2 Eqr Kr Kr Eqr   , )( 222 2 2 Eqr  , dEqra r K  )( 2222 2 221 2   . www.scholink.org/ojs/index.php/asir Applied Science and Innovative Research Vol. 1, No. 1, 2017 6 Published by SCHOLINK INC. Now, when 2P exists, 0)( 222  Eqr . Thus the above eigenvalues are negative i.e., the equilibrium at 2P is asymptotically stable if 0        )(,)( 111 1221 12 222 2 222 rEq Kr Kr aK dr MinEqr  and 0)( 111  Eqr (18) At )0,0),(( 111 1 1 3 Eqr r K P  the Jacobian matrix (15) reduces to                       dEqr r Ka EqEqr rK K r Eqr r Ka EqrEqr J )(00 0))(1(0 )()()( 111 1 111 22111 12 121 2 111 1 11 11112111 3    whose eigenvalues are, 0)( 111 1 3  Eqr , )()( 11121 21 12 222 2 3 Eqr Kr Kr Eqr   , dEqra r K  )( 1111 1 113 3   . Thus the eigenvalues are negative real numbers that is the equilibrium at 3P is asymptotically stable if        )(,)(0 222 2112 21 111 1 111 rEq Kr Kr aK dr MinEqr  and 0)( 222  Eqr . (19) At ))1(,0,( 1 11 1111 1 11 4 a Eq aK d a r a d P   ),0,( 444 zxP , say, the Jacobian matrix (15) becomes                        0 0)()1(0 422411 2 111 1 1 111 222 112 21 2 1111 121 111 1 4 zaza aK dr a Eqr aEq aK d r d aK dr aK dr J       . One of whose eigenvalues is )()1( 2 111 1 1 111 222 112 21 2 1 4 aK dr a Eqr aEq aK d r       and other two eigenvalues are the roots of 041 111 12  dza aK dr    . The sum of whose roots = 0 111 1  aK dr  , and product of the roots = 041 dza under the condition of existence of 4P . Thus by Routh-Hurwitz rule the equilibrium at 4P is asymptotically stable if 01 4  when the equilibrium at 4P exists. Thus 4P is a stable equilibrium if, )()( 111 1 111 1 2 112 212 222 aK dr Eqr a a aK dr Eqr    where 111 1 111 )( aK dr Eqr   (20) At the equilibrium point )0,,( 215 NNP where www.scholink.org/ojs/index.php/asir Applied Science and Innovative Research Vol. 1, No. 1, 2017 7 Published by SCHOLINK INC. )]()([ 1 1 22212 2 2 111 1 1 2112 1 Eqr r K Eqr r K N      , )]()([ 1 1 222 2 2 11121 1 1 2112 2 Eqr r K Eqr r K N      , with 111 Eqr  > 0, 222 Eqr  > 0 and 12112  , the Jacobian matrix (15) reduces to                       dNaNa Na K Nr K Nr Na K Nr N K r J 222111 22 2 22 2 2212 11 1 1121 1 1 1 5 00    One of the eigenvalues of 5J is dNaNa  222111 1 5  , and other two are the roots of the quadratic 0)1()( 21 21 21 2112 2 22 1 112  NN KK rr K Nr K Nr  . The sum of whose roots = 0)( 2 22 1 11  K Nr K Nr , and the product of the roots = 0)1( 21 21 21 2112  NN KK rr  under the conditions of existence of 5P . Thus by Routh-Hurwitz rule the equilibrium point 5P is asymptotically stable if 0222111  dNaNa  , that is if )1())(())(( 2112 2 222 12 2 211 22221 1 122 1 111 111              d r Ka r Ka Eqr r Ka r Ka Eqr (21) together with the condition of existence of 5P . The community matrix at ),,( * 3 * 2 * 16 NNNP is                      0* 322 * 311 * 22 2 * 22 2 * 2212 * 11 1 * 1121* 1 1 1 6 NaNa Na K Nr K Nr Na K Nr N K r J    . (22) The eigenvalues )3,2,1(6 ii of 6J are the roots of the following characteristic cubic 0. 66 2 6 3  JJadjTraceJTrace  By Routh-Hurwitz rule, this cubic has roots with negative real parts if 0,0 66  JJTrace , and 666 .. JJadjTraceJTrace  . Now,  6J )( * 3 * 2 * 1 NNN 0)()( 1 2 122112 2 2 12112 2 22 1 1         aaa K r aaa K r ,  6JTrace 0)( 2 * 22 1 * 11  K Nr K Nr , and  6.JadjTrace * 2 * 12112 21 212 11 * 3 * 1 2 22 * 3 * 2 )1( NN KK rr aNNaNN   . Therefore the www.scholink.org/ojs/index.php/asir Applied Science and Innovative Research Vol. 1, No. 1, 2017 8 Published by SCHOLINK INC. interior equilibrium at ),,( * 3 * 2 * 16 NNNP is asymptotically stable if    , that is if )( )1()( 2 2122 1 1211 21 * 3 * 2 * 1 * 2 * 12112 21 21* 2 2 2* 1 1 1* 3 * 2 2 22 2 2* 3 * 1 2 11 1 1 22 K r K r aaNNN NN KK rr N K r N K r NNa K r NNa K r     (23) 5. Optimal Harvesting Policy Once the process of harvesting the resources is started, the problem of management of the fisheries can be viewed in terms of rent maximization. Now we shall use Pontryagin’s maximum principle to solve our optimization problem and obtain optimal harvest policies )(,)( 21 tEtE such that our objective functional is maximized. Importance of discount rate cannot be underestimated in addressing environmental and resources issues. In fact the optimal stock size for a given fishery will vary depending on discount rate. Now we assume that 1p Constant price per unit of biomass of 1N species, 2p Constant price per unit of biomass of 2N species, 21 ,cc are constant costs of fishing for the two facultative prey species per unit effort. Our objective is to study the optimal harvest policy, i.e., ],0[)( max ii EtE  , i=1, 2 and   ,0t to maximize the profits of harvesting agencies and to keep the populations at an optimum level. The objective functional representing the present value J of a continuous time stream of revenue is given by     0 22112222111121 )[),( dtEcEcNEqpNEqpeEEJ t (24) where  is instantaneous annual rate of discount. In this section we will find an optimal harvest policy solving the following optimization problem, Maximize ),( 21 EEJ Subject to (3)-(5) and 2,1;)(0 max  iEtE ii The Hamiltonian for this control problem is taken as H= ][ 221122221111 EcEcNEqpNEqpe t  + 1 [ 111311 1 2121 11 )1( NEqNNa K NN Nr     ] + 2 222322 2 1212 22 )1([ NEqNNa K NN Nr     ] + 3 ][ 332223111 dNNNaNNa   (25) where 3,2,1,)(  itii  are adjoint variables. We have, 1E H   = )(][ 11111111 tNqcNqpe t   (26) www.scholink.org/ojs/index.php/asir Applied Science and Innovative Research Vol. 1, No. 1, 2017 9 Published by SCHOLINK INC. and 2E H   = )(][ 22222222 tNqcNqpe t   (27) The optimal control 2,1),( itEi must clearly satisfy the condition  0)(0,0)()( max  twhenandtwhenEtE iiii  (28) Since )(ti causes )(tEi to switch between the levels 0 and maxE , )(ti (i=1,2) are called switching functions. Depending on the switching functions )(ti , the optimal control )(tEi is a bang-bang control switching from one extreme level to other one. Once )(ti (i=1,2) vanishes ,the Hamiltonian functions H becomes independent of the control variable )(tEi (i=1,2) and its optimal value cannot be determined by the above procedure. It is then called singular control max** )(0,)( iii EtEtE  . Hence our optimal harvest policy becomes          0)( 0)(0 0)( )( * max ii i ii i fortE tfor tforE tE    for i=1, 2. For singular control )(ti =0; i=1, 2. Then from (26) and (27) we find that the values of 1 and 2 are )( 11 1 11 Nq c pe t   (29) )( 22 2 22 Nq c pe t   (30) Now the adjoint equations are  dt d 1 1N H   , 2 2 N H dt d     , 3 3 N H dt d     (31) Then using (25) and the third equation of (31) we find )]([ 2221113222111 3 dNaNaNaNa dt d    Using (29), (30) we write,  dt d 3  11 11 1 1 )( Na Nq c pe t 22 22 2 2 )( Na Nq c pe t  3   11 11 1 1 )[( Na Nq c p e t   ])( 22 22 2 2 Na Nq c p  (32) We here consider that the constant of integration vanishes so that the shadow prices )3,2,1( ie t i  of three species are bounded. Using (29) and (30) in other two equations of (31) we find that, ][)( 31132122 2 2 11 1 1 111 11 1 1 NaNr K rN K Eqpe Nq c pe tt       , and, www.scholink.org/ojs/index.php/asir Applied Science and Innovative Research Vol. 1, No. 1, 2017 10 Published by SCHOLINK INC. ][)( 32231211 1 1 22 2 2 222 22 2 2 NaNr K rN K Eqpe Nq c pe tt       . Now using the values of the adjoint variables 321 ,,  in the above equations we get, ])()()()([)( 2 22 2 221 11 1 11 311 221 22 2 2 2 2 1 11 1 1 1 1 111 11 1 1        N Nq c paN Nq c pa Na N Nq c p K r N Nq c p K r Eqp Nq c p    and, ])()()()([)( 2 22 2 221 11 1 11 322 2 22 2 2 2 2 121 11 1 1 1 1 222 22 2 2        N Nq c paN Nq c pa Na N Nq c p K r N Nq c p K r Eqp Nq c p    Therefore, ))(())(( 2 221232211 22 2 2 31 2 11 1 11 11 1 1111 K NrNNaa Nq c p NNa K Nr Nq c pEqp        , and ))(())(( 31 2 22 2 22 22 2 2 1 112131212 11 1 1222       NNa K Nr Nq c p K NrNNaa Nq c pEqp  From these we get the following expressions of harvesting efforts, )])(())([( 1 2 221232211 22 2 2 31 2 11 1 11 11 1 1 11 1 K NrNNaa Nq c p NNa K Nr Nq c p qp E        (32) )])(())([( 1 31 2 22 2 22 22 2 2 1 112131212 11 1 1 22 2       NNa K Nr Nq c p K NrNNaa Nq c p qp E  (33) Hence solving the steady state equations together with (32) and (33) we get an optimal equilibrium solution   321 ,, NNN and optimal harvesting efforts 1E and 2E . 6. Numerical Simulation Example 1. For numerical analysis we first consider the following set of values of parameters, 2.0,4.0,3.0,12,05.0,02.0,2.0 ,100,5.2,10,04.0,01.0,1.0,100,09.2 21222`21 221111211   dEqa KrEqaKr   Then the equilibrium points are 0P (0,0,0), 1P (0,25,63.75), 2P (0,76,0), 3P (80.86124,0,0), 4P (66.6666,0,29.6666), 5P (90.2666,94.0533,0), 6P (38.2647,10.6507,91.2527). www.scholink.org/ojs/index.php/asir Applied Science and Innovative Research Vol. 1, No. 1, 2017 11 Published by SCHOLINK INC. 0 10 20 30 40 50 60 70 80 90 100 0 20 40 60 80 100 120 time p o p u la ti o n N 1 Species N 2 Species N 3 Species Figure 1. Solution Curves of the Species 0 50 100 150 200 250 300 350 400 450 500 10 20 30 40 50 60 70 80 time p o p u la ti o n N 1 Species N 3 Species Figure 2. Solution Curves When the Prey Species N2 is Absent www.scholink.org/ojs/index.php/asir Applied Science and Innovative Research Vol. 1, No. 1, 2017 12 Published by SCHOLINK INC. 0 50 100 150 200 250 300 350 400 450 500 10 20 30 40 50 60 70 time p o p u la ti o n N 2 Species N 3 Species Figure 3. Solution Curves When the Prey Species N1 is Absent 0 20 40 60 0 20 40 60 80 0 20 40 60 80 100 120 N 1 -SpeciesN 2 -Species N 3 -S p e c ie s Figure 4. Phase Diagram in Terms of the Values of the Parameters Taken in Example 1 In Figure 1, we see the solution curves which exhibits the stability of the system at 6P and in Figure 2, Figure 3 we see that in absence of one prey the other survives in presence of the predator also. Thus, this numerical example supports our hypothesis that the facultative mutualism is not disturbed in the www.scholink.org/ojs/index.php/asir Applied Science and Innovative Research Vol. 1, No. 1, 2017 13 Published by SCHOLINK INC. presence of the predator. In Figure 4, we see that for the considered values of the parameters the system is also globally stable. Example 2. We consider the following set of values of parameters for numerical analysis of optimal equilibrium. 04.0,6,5,15,10,2.0,4.0,3.0,05.0 ,02.0,2.0,100,5.2,04.0,01.0,1.0,100,09.2 2121212 2`2122111211     ccppdq aKrqaKr Then optimal harvesting efforts are, 9632.131 E , 4487.322 E and corresponding optimal equilibrium point is )823.48,082.6,529.50(),,( 321  NNN . From the Figure 7 and Figure 8, it can be realized that the mutualism between the prey species remains facultative even when optimal harvesting efforts are used. 0 20 40 60 80 100 120 140 160 180 200 0 10 20 30 40 50 60 70 time p o p u la ti o n N 1 Species N 2 Species N 3 Species Figure 5. Optimal Solution Curves www.scholink.org/ojs/index.php/asir Applied Science and Innovative Research Vol. 1, No. 1, 2017 14 Published by SCHOLINK INC. 0 20 40 60 80 100 0 20 40 60 80 0 20 40 60 80 100 N 1 -SpeciesN 2 -Species N 3 -S p e c ie s Figure 6. Phase Diagram Corresponding to the Optimal Harvesting Efforts 0 50 100 150 200 250 300 350 400 450 500 10 12 14 16 18 20 22 24 26 time p o p u la ti o n N 2 Species N 3 Species Figure 7. Solution Curves Showing Survival of N2 Species in Absence of N1 Species www.scholink.org/ojs/index.php/asir Applied Science and Innovative Research Vol. 1, No. 1, 2017 15 Published by SCHOLINK INC. 0 10 20 30 40 50 60 70 80 90 100 10 20 30 40 50 60 70 time p o p u la ti o n N 1 Species N 3 Species Figure 8. Solution Curves Showing Survival of N1 Species in Absence of N2 Species 7. Conclusion In this paper, we have presented a mutualism model with independent harvesting efforts on mutualist prey species in presence of a predator. During equilibrium analysis and numerical simulation we have seen that the predator survives even when one species is extinct. The important fact is that the mutualism remains facultative in presence of the predator. This is due to the survival of one species in absence of other and in presence of the predator. We derive the optimal harvesting policy and numerically the optimal equilibrium, optimal harvesting efforts are obtained. Numerical illustrations show that, the mutualism remains facultative even when optimal harvesting efforts are used on prey species. This paper includes simple linear predation functions and does not include obligate mutualism between the prey species. We used Mat Lab for numerical calculations and graphs. Acknowledgement The research of Dr. Saroj Kumar Chattopadhyay is financed by the University Grants Commission of India (F.NO.PSW-184/14-15 (ERO)). References Clark, C. W. (1985). Bioeconomic Modelling and Fisheries Management. Wiley, New York. Clark, C. W. (1990). Mathematical Bioeconomics: The Optimal Management of Renewable Resources. John Wiley & sons, New York. Dai, G., & Tang, M. (1998). Coexistence region and global dynamics of harvested predator-prey system. www.scholink.org/ojs/index.php/asir Applied Science and Innovative Research Vol. 1, No. 1, 2017 16 Published by SCHOLINK INC. SIAM J. Appl. Math., 13, 193-210. https://doi.org/10.1137/S0036139994275799 Janzen, D. H. (1985). The natural history of mutualisms. In D. H. Boucher (Ed.), The Biology of Mutualism (pp. 40-99). Oxford University Press, Oxford. Kar, T. K., & Chaudhuri, K. S. (2004). Harvesting in a two prey one predator system: A bioeconomic model. ANZIAM J., 45, 443-456. https://doi.org/10.1017/S144618110001347X Kot, M. (2001). Elements of Mathematical Ecology. Cambridge University Press. https://doi.org/10.1017/cbo9780511608520 Lengeler, J. W., Drews, G., & Schlegel, H. G. (1999). Biology of the Prokaryotes. Blackwel Science, Stuttgart. Margulis, L. (1970). Origin of Eukaryotic Cells. Yale University Press, New Haven. Margulis, L. (1981). Symbiosis in Cell Evolution. W. H. Freeman, San Francisco. May, R. M. (1981). Theoretical Ecology: Principles and Applications. Sinauer Associates, Sunderland, MA. Mesterton, G. M. (1988). On the optimal policy for combined harvesting of predator and prey. Natural Resources Modelling, 3, 63-90. Strobele, W. J., & Wacker, H. (1995). The economics of harvesting predator-prey system. J. Econ., 61, 65-81. https://doi.org/10.1007/BF01231484 Wallin, I. E. (1923).The mitochondria problem. American Naturalist, 57, 255-261. https://doi.org/10.1086/279919 Wallin, I. E. (1927). Symbionticism and the origin of Species. Williams & Wilkins, Baltimore, MD. https://doi.org/10.5962/bhl.title.11429