Modeling and Simulation Study of the Population Dynamics of Commensal-Host-Parasite System

This paper deals with the modeling and simulation study of Commensal-host species system together with the inclusion of parasite population. The model comprises of three populations viz. Host, Commensal and Parasite. The Commensal population gets benefit from Host population but the former do not do any harm to the latter. The parasite population gets benefit and also do harm to Host population. However, the Commensal population only harms the parasites. The mathematical model is comprised of a system of three first order non-linear ordinary differential equations. Mathematical analysis of the model is conducted. Positivity and boundedness of the solution have been verified and thus shown that the model is physically meaningful and biologically acceptable. Scaled model is constructed so as to reduce the number of model parameters. Equilibrium points of the model are identified and stability analysis is conducted. Simulation study is conducted in order to support the mathematical analysis. In the present model the Commensal population lies higher and the parasite population lies below respectively the host population. This fact is well supported by the mathematical analysis as well as simulation study. The results of analysis and simulation are presented and discussed lucidly in the text of the paper.


Introduction
It has been the general observation that populations can sometimes undergo drastic changes in their abundance. These changes are not only unexpected but also are equally difficult to explain in simple manner either theoretically or empirically [1][2]. Many species have been driven to extinct and many others are at the verge of extinction due to several reasons such as overexploitation, predation, and the like [3]. In nature populations may interact in different ways with one another during their lifetime and even they may undergo from one type of interaction to another.
The species populations in an ecosystem will interact in different ways and these interactions can be regarded as positive (+), negative (-) or neutral (0). For instance the interactions such as mutualism and commensalism are positive interactions while those such as intra and interspecific competitions, parasitism and predation are considered as negative interactions.
Mathematical representations have been developed to address all these kinds of interactions.
In general, in case of a two-species interaction the mathematical formulation can be expressed by a system of two equations as Here in (1) , 1, 2 but . The parameters and represent the intrinsic growth rate and carrying capacity of the species respectively. The variables and ! represent the population densities of the species and respectively. The quantity " ! is the interaction coefficient between species and .
Some classical forms of the models addressing two interacting populations are represented by these types of system of equations depending on the signs of the interaction coefficients" ! . In general, the positive sign of interaction coefficient indicates commensalism and mutualism while the negative sign of that indicates competition and predation. Also, the negative interaction includes all predator-prey and parasite-host relationships [4,5]. In general, the Host-Parasite interaction or parasitism is referred to as negative.
According to equation (1), in absence of the species , the species will grow logistically and reaches its equilibrium value or carrying capacity over the timefrom any positive initial density and continues to persist for ever. From this the following two facts can be observed: (i) the interactions between the species with are facultative. That is, species can exist with orwithout being associated to species j and (ii) density-dependent mechanisms of species alone govern the dynamics of the species .
However, it is well known that some symbiont is biologically incapable of existing without a partner: Let it be a mutualistic or a parasite or a predator. If the species cannot persist alone in absence of species , thenthe former's population size is expected to decrease exponentially down to zero. To account the absence of species , the equation (1) is modified by setting N = 0 as The literature on the mathematical modeling of the negative interactions is very abundant [6][7][8][9]. But the modeling of positive interactions, particularly commensalistic ones have generated interest among researchers in recent years.
Historically, the mutualistic interaction model that is described by Monod function form is elucidated by the following systems of ordinary differential equations: (3) Now based on these mutualistic model equations (3) -(4), in the present work it is proposed to derive an ecological model that describes commensalism.
Recall that commensalism (0, +) is a relation in which one of the species known as the Commensal gets benefit without harming or affecting the other called the host.
The reduction of mutualism to commensalism leads to setting of conditions on interaction coefficients as " '+ = 0 and " +' > 0 . Thus, the mutualistic model becomes a commensal model and the system (3) and (4) reduces to the form as This model represents increase of N + population when N ' increases. That is, the growth of N + population is somehow proportional to N ' . However, the growth of N + will neither harm nor influence in any manner that ofN ' .
The main objective of the present research is to study the dynamical behavior of a negative interaction between species assuming that the growth of commensal follows Monod functional response [10][11].
The paper is structured in a systematic manner. In section 2, the basic assumptions and mathematical formulations of the model are presented. In section 3, the results on twodimensional sub-communities are briefly presented. The modeling is extended to three-dimensional communities and incorporated mathematical analysis and other results including stability analysis of equilibrium points. Finally, numerical examples are carried out and concluding remarks of the study with a brief discussion is presented in Section 4.

Basic Assumptions and Model Formulation
Mathematical ecology and epidemiology are major fields of study that have been separately treated by different researchers [3,[12][13]. Since transmissible disease in ecological situation cannot be ignored, it is very important from both the ecological and the mathematical points of view that the dynamics of ecological populations be studied subjected to epidemiological factors. A large number of studies have been conducted in this field of research [14][15]. In the common life, the disease may spread among the host and the commensal. Mathematical model plays an ever more important role in the study of eco-epidemiology. Furthermore, it provides understanding of the underlying mechanisms that influence the spread of disease and in the process it also suggests control strategies [6,[16][17][18][19][20].
The present study is conducted in a systematic way. Initially, simple mathematical models for a two-dimensional Commensal-host and host-parasite systems are framed separately. Finally, the model is extended to threedimensional case which is of prime importance of this study [21]. The model equations for two species host-Commensal, host-parasite systems and then for three species system will be presented in terms of system of first order non-linear differential equations.
Let us start with considering the interaction between two species in which commensal would be benefited without harming the host. Both the commensal and host populations would follow logistic functional growths if the two populations live separately. However, in presence of host the logistic growth of commensal population would be influenced by Monod response function. The population dynamics of such commensal-host system can be mathematically modeled as Let us continue with considering the interaction between two species in which parasite would be benefited by harming the host. The host population grows following logistic function and the parasite population decreases proportional to its size if the two populations live separately. However, if both the host and parasite populations live together their growths will be influenced by Monod response functions. The population dynamics of such host-parasite system can be mathematically modeled as In this study, considering and carefully studying the mathematical modeling of two dimensional systems of commensal-host and host-parasite, the concept is extended to the three dimensional commensal-host-parasite systems.
Hence, the three dimensional commensal-host-parasite systems can be mathematically modeled using first order nonlinear ordinary differential equations as follows: 3 µ z k9z, y: αh9z, x: Here in (11)(12)(13), the functionsA9y, z:and k9z, y: represent the interactions between host and parasite populations where as g9x, z: and h9z, x: representthe interactions between commensal and parasite populations with varying coefficients of interactions.
Majority of the models assume that the parasite population increases the death rate of the host population and enhance the commensal population. This fact is visualized in the model equations too.

Notation
Description ' , ' Growth rate and carrying capacity of commensal population 9 ' ' ⁄ : Self inhabitation coefficient of commensal population + , + Growth rate and carrying capacity of host population 9 + + ⁄ : The interaction coefficient of commensal-host populations C + , C D The interaction coefficients of host-parasite populations E Total death rate of parasites F ' Response rate between host and commensal F + , F D Response rates between host and parasite G Interaction coefficient between commensal and parasite The model variables H9I:, J9I: and K9I: represent the sizes of commensal, host and parasite populations respectively at any time t.
The three model equations for 9dx/dt:, 9dy/dt: and 9dz/ dt: are constructed using several components, both variables and parameters, and each of which represents specific biological assumptions.
The simple schematic interactions among the model variable is shown in the figure below.

Two Dimensional Sub-systems
Now the commensal-host and parasite-host models are stated in what follows.
The ecological system under consideration consists of three species namely (i) Commensal population O9I: (ii) host populationP9I: and (iii) parasite population Q9I:.
The Commensal-Host system: The host commensal model with type I functional response and carrying capacities is expressed as [ However, this model does not include any other external forces like parasites. To fill the gap the present study modifies this model by including parasite population.
The classical Lotka -Volterra model describes exponential growth of host and exponential decay of parasite populations. Also, the model equations accommodate Holling type I functional responses. Recall that, Holling type I functional response is linear. The simple form of the Lotka -Volterra model is given in terms of ODEs as [6,22] 0 y 91 εz: The classical Lotka -Volterra model is structurally unstable since the model allows the population to grow indefinitely and this is very unnatural and unrealistic. However, it can be used as groundwork for constructing more realistic representations [23]. To overcome this instability carrying capacity is required to be incorporated. The, interaction between host-parasite system can be modified with the inclusion of carrying capacity as In the modified model the response function is still linear. But type II response is more realistic than type I functional Commensal-Host-Parasite System response. Thus, the model can be further modified by replacing type I functional response by type II response as The model (20) can be interpreted that the infected hosts are immediately eradicated as it is represented by the functional response -εf(y, z) in the host equation. This implies that the infected hosts have no ability to recover from the disease. Now, in the present work the two species systems (14), (15), (20) and (21) are combined and formulated a three species system with the incorporation of the following features: (i)carrying capacities and functional responses are assigned to both the species: Commensal and Host (ii) no recovery is possible for the Host species once they get infected. It is also possible that the commensal species will change their behavior in a long run of time and act as parasites. However, this possibility is not considered here in the present study and the authors plan to take it up as their next study.

The Three Dimensional Systems
It is well known that Epidemiology and Ecology are two major and different fields of research. There are many epidemiological models and many ecological models, but certainly few models of the eco-epidemiology [1-3, 12-13, 23-24].
The presence of logistic and functional response terms in the host equation can be interpreted s the natural growth rate and mortality rate due to parasites respectively. Here arise three cases: (i) if the logistic term dominates the functional response term the host population grows with positive growth rate (ii) if both the logistic and the functional response terms are equal then growth rate of the host population becomes zero. That is the host population size remains constant and (iii) if the functional response term is allowed to dominate logistic term then the host population grows with negative growth rate and eventually the population size extinct. This situation is not desirable. This possibility can be avoided by decreasing parasite population through the inclusion of parasite-commensal interaction [10][11].
A three dimensional host-commensal-parasite system is formulated here to describe the host-parasite interactions. The model also includes the logistic growth rates for host and commensal populations. The total mortality rate of parasites due to natural death, increment of commensal and reduction of hosts is also accommodated.
In view of the above, it is clear that the introduction of commensal population into host-parasite system plays a significant and complex role in the dynamical advancement of host-parasite interaction [11].
The objective of this study is to identify the influence of commensal on the dynamics of host-parasite interaction in the context of eco-epidemiology. The model, including all the aforementioned assumptions, for commensal-host-parasite interaction is constructed by combining systems (14), (15), (20) and (21). It can be expressed as a system of differential equations as follows: Here in this system, the initial conditions of the model variables areO T > 0, P T > 0, Q T > 0.
To reduce the number of model parameters and to determine which combinations of parameters control the behavior of the system, the system (22), (23) and (24) is dimensionalize with the following scaled parameters: After manipulation the system (22), (23) and (24) Further, all the model parameters and hence the scaled parameters are supposed to be positive constants.

Dynamical Behavior of the System
Since the state variables O, P and Q represent population sizes, positivity implies that the population sizes never become negative. The boundedness may be interpreted as a natural restriction to the growth of populations as consequences of limited resources.
Here, some basic dynamical properties of the system are discussed subjected to positive initial conditions. Let b ( denote the set of all non-negative real numbers and b D ( = cH ∈ b D : H = (O, P, Q) "fg O, P, Q ∈ b ( h . When the functions on the right-hand of the model equations (22), (23) and (24) are denoted by a vector as i = (A ' , A + , A D ) then clearly i ∈ j(b ( D ). The fundamental theorem of existence and uniqueness assures the existence and uniqueness of the solution for the model system together with the given initial conditions.

Positivity of the Solution
Here the positivity of each population size such as O(I), P(I) and Q(I) is verified. These system variables must have the positive values in order to be biologically meaningful. The positivity of these biological or system variables is tested and the results are presented in the form of a lemma as follows: Lemma 1 Every solution of system (22), (23) and (24) (22) is solved analytically and its solution is obtained as The exponential function is always non-negative and the initial commensal populationx T is assumed to be positive. Therefore, O(I) > 0 for allI ≥ 0.
Verifying the positivity of P(I): The host equation of the system, (23) and is solved analytically and its solution is obtained as y(t) = y T exp t ur + 1 − y k + − β + z α + + y v dt The exponential function is always non-negative and the initial host population y T is assumed to be positive. Therefore, P(I) > 0 for allI ≥ 0.
Verifying the positivity of Q(I): The parasite equation of the system (24) is solved analytically and its solution is obtained as The exponential function is always non-negative and the initial host population z T is assumed to be positive. Therefore, Q(I) > 0 for allI ≥ 0.
Hence, all the solutions of the system (22), (23) and (24) are positive for all t≥ 0 under the considered positive initial conditions.

Boundedness of the System
In this section, all the solutions of system (22), (23) and (24) are shown to be bounded. The boundedness of the system is presented in the form of a lemma as follows: Lemma 2 All solutions of the system (22), (23) and (24)  Here the quantityc is an arbitrarily integral constant and must be strictly positive due to logarithmic function. Now there arise three cases, namely the exponent may be negative or positive or zero. These cases as analyzed as follows: [ & % < 0 then as t → ∞ the exponential term takes a value zero and thus to get z ≤ 0 . But this, having negative parasite population, is biologically not feasible. Hence this possibility is not considered.  Here the quantity c + is an arbitrarily integral constant and must be strictly positive due to logarithmic function. Therefore, the solution of the model system (22)(23)(24) is bounded.

Existence of Equilibrium Points
The scaled system of model equations now can be expressed as: The equilibrium points of the system can be obtained by setting  It can be observed that the first two combinations lead to the same equilibrium point" T while the last combination leads to biologically infeasible solution. Thus, six equilibrium points are possible. Now the coordinates of these six possible equilibrium points as given below: (1) E T 90, 0, 0) Trivial steady state A T = k T n D − n ' k ' + n + ;n ' = θσh + m ' + , n + = θ + kβm ' (1 + σ), n D = θ + σkβ, A ' = n • − n ' n -;n • = θ + βσm ' + ,n -= m ' (θσ − σ + h + ), A + = k T k ' + nn -+ n D ;k T =2θσm ' − θσh + − θσh + m ' + m + kσ + kθm ' − kmε.

Variation matrix
To analyze the stability of the equilibrium points the community matrix called Jacobian matrix is determined as follows: Let -= h( x, y, z ), 0 = f( x, y, z )and 3 = g( x, y, z ). Also, the components of the Jacobian matrix are given
Thus E D is stable when these conditions are satisfied otherwise unstable.
Proof: The characteristic equation of the variation matrix (J E • ) is given by; Therefore, all the roots of equation (36) will have negative real roots if È = P Ð and by Hurwitz Routh criteria of sufficient condition for stability, E -(x -, y -, z -) is stable.

Numerical Simulations
Numerical examples or simulation study is used to confirm and support the results of mathematical or theoretical analysis.
That is, in order to verify the theoretical predictions of the model the numerical simulations are carried out. For that purpose a set of biologically valid parametric and initial values are used.
All parameter values chosen here are realistic and at the same time they obey the stability conditions.
The pictorial representation of the simulation study of population dynamics of host-parasite model in absence of commensal population is given in Figure 2. In Figure 2 the parametric values used are, h + 0.020, k 0.500, m ' 0.200, σ=0.9, ε 0.08, m + 0.60. Figure 2 shows the population dynamics of two dimensional host-parasite systems in absence of the commensal population. It is expected that the host population grows logistically but decays due to parasites. Similarly, the parasite population decays naturally but grows due to host. In the figure three phases of population dynamics are observed: (i) Phase 1 is the time interval from origin to the time of occurrence of maximum value for P (ii) Phase 2 is the time interval from the occurrence of maximum value for P to that for Qand (iii) Phase 3 is the time interval from the occurrence of maximum value for Q. In phase 1, both host and parasite populations increase. In phase 2, the host population decreases while the parasite population increases. In phase 3, both host and parasite populations decrease.
The following observations can be made from  In Figure 3    In Figure 5  Here as the time progresses population size of the parasite increases higher and the population sizes of both host and commensal converge to their respective carrying capacities. In this case the host population is dominated by both commensal and parasites. That is the dependents are more. Also, the disease is more in the system. However, this dynamics appears not interesting. Initially the population sizes fluctuate more but as the time progresses they converge to their optimum or stable values. Ultimately the population sizes maintain some balance among themselves. Hence, the dynamics shown in this simulation are highly accepted. The figure also supports that the interior equilibrium point is stable.

Conclusions
In the present study the population dynamics of a three dimensional commensal-host-parasite system have been taken up for investigation. A mathematical model for this three-dimensional system is developed by studying carefully and incorporating the futures of two-dimensional systems. The positivity and boundedness of the model variables are verified and hence shown that the newly developed threedimensional system is biologically valid. To reduce the model parameters and to highlight the more important parameters, the scaled parameters have been introduced and the scaled model equations are constructed. Six equilibrium points have been identified and following Jacobian matrix method their stability is tested. Though the calculation of interior equilibrium point involves very complex algebraic operations, it appeared that it is stable. However, the simulation study conducted promises that the interior equilibrium point is stable. In this paper, the model presented here only involves commensal, host and parasite, and only considers the case that the Monod type responses function.
It is interesting to consider the interaction between multicommensal, multi-hosts and multi-parasites and the other interaction of the parasites in the model, and study the dynamical behaviors of the multi-hosts and multi-parasites system. Beside this, in the model harvesting rates as well as recovery of the population is not assumed. These assumptions will be considered in a future study.
In the present model the commensal has all the three positive growth rates due to nature, host and parasite; the host grows naturally but decays due to parasite; the parasite grows due to host but decays naturally and due to Commensal. In other words it can be put as; Commensal has no negative growth rate; Host has one negative growth rate and the Parasite has two negative growth rates. Hence the following results: (i) The population size of Commensal always lies above that of both the host and parasite. This fact is well supported by the simulations of the figures (3) to (6).
(ii)Though the Commensal population lies higher, the host and parasite compete for the second position.