\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2015 (2015), No. 02, pp. 1--20.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2015 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2015/02\hfil Stability of mutualisms] {Stability of mutualisms in a lattice gas system of two species} \author[Y. Wang, H. Wu \hfil EJDE-2015/02\hfilneg] {Yuanshi Wang, Hong Wu} % in alphabetical order \address{Yuanshi Wang \newline School of Mathematics and Computational Science, Sun Yat-sen University, Guangzhou 510275, China} \email{mcswys@mail.sysu.edu.cn} \address{Hong Wu (corresponding author) \newline School of Mathematics and Computational Science, Sun Yat-sen University, Guangzhou 510275, China} \email{wuhong@mail.sysu.edu.cn} \thanks{Submitted November 19, 2014. Published January 5, 2015.} \subjclass[2000]{34C37, 92D25, 37N25} \keywords{Stability; persistence; cooperation; saddle-node bifurcation; \hfill\break\indent Holling Type II functional response} \begin{abstract} This article considers mutualisms in a lattice gas system of two species. The species are mutualistic since each one can provide resources to the other. They are also competitive since they compete for empty sites on the same lattice. The mutualisms are assumed to have a saturated response, and the intraspecific competition is considered because of self-limitation. The mutualism system is characterized by differential equations, which are derived from reactions on lattice and are extension of a previous model. Global stability analysis demonstrates that (i) When neither species can survive alone, they can coexist if mutualisms between them are strong and population densities are large, which exhibits the Allee effect in obligate mutualism; (ii) When one species can survive alone but the other cannot, the latter one will survive if the mutualistic effect from the former is strong. Even if the effect is intermediate, the latter species can survive by strengthening its mutualistic effect on the former and enhancing its population density; (iii) When either species can survive alone, a weak mutualism will lead to extinction of one species. When in coexistence, intermediate strength of mutualism is shown to be beneficial under certain parameter range, while over- or under- mutualism is not good. Furthermore, extremely strong/weak mutualism is exhibited to result in extinction of one/both species. While seven typical dynamics are displayed by numerical simulation in a previous work, they are proved in this work and the eighth one is exhibited. Numerical simulations validate and extend our conclusions. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \allowdisplaybreaks \section{Introduction} Mutualistic interactions are ubiquitous in nature since most biomass survive by cooperating with the others \cite{Bego06}. For example, many microbial species are observed to play a role in the abundance of interrelated species \cite{Kell06, Madi00}, while various bacterial species coexist in syntrophic colonies, in which one species consumes resources produced by another (obligate mutualisms) \cite{Agui07,Shap95}. Many models have been presented to characterize mutualisms, among which Lotka-Volterra equations (LVEs) are the most famous \cite{Smith95,Smith11}. The LVEs can be modeled by \begin{equation} \label{eq11} \begin{gathered} \frac{dx_1}{dt} = x_1(r_1-d_1x_1 + e_1x_2) \\ \frac{dx_2}{dt} = x_2(r_2-d_2x_2 + e_2x_1), \end{gathered} \end{equation} where variable $x_i$ represents population density of species $i$, while parameters $r_i$ and $ d_i$ denote the intrinsic growth rate and self-competition degree in species $i$, respectively $(i=1,2)$. $e_i$ represents the mutualistic effect of species $j$ on $i$, $i \ne j$, $i,j=1,2$. It is known that the two species can coexist at a steady state if mutualistic effects between them are weak ($e_1e_20, r_2 > 0$, while cases of $r_1 >0, r_2=0$ and $r_1=r_2 = 0$ are considered in Appendices A and B respectively. \begin{figure}[htb] \begin{center} \includegraphics[width = .6\textwidth]{fig1} \end{center} \caption{Intersection points of parabolas $L_1$ and $L_2$ on the whole plane. Let $m_1=m_2=0.02$, $e_1=e_2=40$, $c_1=c_2=0.001$, $d_1=0.8$, $d_2=0.75$. Then $L_1$ and $L_2$ have four intersection points, while two of them are in the second and fourth quadrants, respectively.} \label{fig1} \end{figure} Assume $r_1>0, r_2 > 0$. Then each species can reproduce in the absence of the other as shown in \cite{Iwat11}. Denote \begin{equation} \label{eq20} m_1 := \frac{m_1}{r_1},\quad m_2 := \frac{m_2}{r_2},\quad e_1 := \frac{e_1}{r_1},\quad e_2:= \frac{e_2}{r_2}\,. \end{equation} Then model \eqref{eq11} can be rewritten as \begin{equation} \label{eq21} \begin{gathered} \frac{dx_1}{dt} =r_1 x_1 [-m_1 + (1+\frac{e_1 x_2}{1+c_1 x_2}) (1-x_1-x_2) -d_1x_1 ],\\ \frac{dx_2}{dt} =r_2 x_2 [-m_2 + (1+\frac{e_2 x_1}{1+c_2 x_1}) (1-x_1-x_2) -d_2x_2]. \end{gathered} \end{equation} Equilibria of \eqref{eq21} are determined by their relative positions of isoclines $L_i$, which can be described as \begin{equation} \label{eq210} \begin{gathered} L_1:-m_1 + (1+\frac{e_1 x_2}{1+c_1 x_2})(1-x_1-x_2)-d_1x_1=0\\ L_2:-m_2 + (1+\frac{e_2 x_1}{1+c_2 x_1})(1-x_1-x_2)-d_2x_2=0. \end{gathered} \end{equation} The expression of $L_1$ can be rewritten as \begin{equation} \label{eq22} \begin{gathered} (\alpha_1 -x_1- \beta_1 x_2) [1 +d_1+ (e_1 +c_1+ c_1 d_1) x_2] = \gamma_1,\\ x_1=x_1(x_2,e_1)= \alpha_1 - \beta_1 x_2 - \frac{\gamma_1}{1 +d_1+ (e_1 +c_1+ c_1 d_1) x_2}, \end{gathered} \end{equation} where \begin{gather*} \alpha_i = \frac{e_i d_i +(e_i+c_i- m_i c_i) (e_i+c_i+ c_i d_i)}{(e_i+c_i+ c_i d_i)^2}, \quad \beta_i = \frac{e_i+c_i}{e_i+c_i+ c_i d_i},\\ \gamma_i = \frac{e_im_i}{e_i+c_i+ c_i d_i} + \frac{e_i d_i (1+d_i + e_i+c_i+ c_i d_i)}{(e_i+c_i+ c_i d_i)^2},\quad i=1,2. \end{gather*} Thus, $P_1((1-m_1)/(1+d_1),0) \in L_1$ as shown in the first quadrant of Figure \ref{fig1}. $L_1$ is a parabola which is convex rightward and has a vertex $\bar{P}(\bar{x}_1,\bar{x}_2)$ with $$ \bar{x}_1 = x_1(\bar{x}_2,e_1),\quad \bar{x}_2 = \sqrt{\frac{\gamma_1}{\beta_1 (e_1 +c_1+ c_1 d_1)} } - \frac{1+d_1}{e_1 +c_1+ c_1 d_1}. $$ The asymptotes of $L_1$ are $$ L_{11}: \alpha_1 -x_1- \beta_1 x_2=0,\quad L_{12}: x_2 = - \frac{1+d_1}{e_1 +c_1+ c_1 d_1}. $$ Since \begin{equation} \label{eq23} \frac{\partial x_1}{\partial e_1} = \frac{x_2 (1-x_1-x_2)}{1+d_1+(e_1+c_1+ c_1 d_1)x_2}>0\quad\text{as }1-x_1-x_2>0 \end{equation} the function $x_1=x_1(x_2,e_1)$ increases monotonously as $e_1$ increases. From \eqref{eq22}, we have $$ \lim_{e_1 \to +\infty}x_1 (x_2,e_1) = 1-x_2,\quad \lim_{e_1 \to 0} x_1 (x_2,e_1) = (1-m_1-x_2)/(1+d_1). $$ If $\alpha_1 \le 0$, then $L_1 \cap\operatorname{int}R_+^2 = \emptyset$ by \eqref{eq22} and $dx_1/dt<0$ by \eqref{eq21}, which implies that species 1 goes to extinction. A similar discussion can be given for species 2. Thus, we assume the following hypothesis in this work $$ \alpha_i>0,\quad i=1,2. $$ When $1-m_1>0$, we have $L_1 \cap\operatorname{int}R_+^2 \ne \emptyset$ since $P_1 \in L_1$ and $L_1$ has an asymptote $\alpha_1 -x_1- \beta_1 x_2=0$. When $1-m_1\le 0$, we have $L_1 \cap\operatorname{int}R_+^2 \ne \emptyset$ if and only if $L_1$ and the positive $x_2$-axis have two intersection points. That is, the following equation has two positive roots $$ G(x_2) \equiv \tilde{a}x_2^2 +\tilde{b} x_2 +\tilde{c} =0, $$ where $$ \tilde{a}= e_1+c_1,~\tilde{b}=1-e_1+c_1(m_1-1),\quad \tilde{c}= m_1-1. $$ Let $\tilde{\Delta} = \tilde{b}^2 - 4\tilde{a} \tilde{c}$, then $$ \tilde{\Delta} = e_1^2 -2 e_1 [ c_1(m_1-1) + 2m_1-1] + [ c_1(m_1-1) -1]^2. $$ Thus the equation $G(x_2) =0$ has two positive roots if and only if $\tilde{b} <0$ and $\tilde{\Delta} >0$. From $\tilde{b} <0$ we have $e_1 >\tilde{e}_1 = c_1(m_1-1)+1$. From $\tilde{\Delta} >0$ we have $e_1>e_1^{+}$ or $e_10, \end{gather*} the equation $G(x_2) =0$ has two positive roots if and only if $e_1>e_1^{(1)}$ with \begin{equation} \label{eqz1} e_i^{(1)} = c_i(m_i-1) + 2m_i-1+2\sqrt{m_i(m_i-1)(1+c_i)},\quad i=1,2. \end{equation} Similarly, the expression of $L_2$ can be rewritten as \begin{equation} \label{eq24} \begin{gathered} (\alpha_2 -x_2- \beta_2 x_1) [1 +d_2+ (e_2 +c_2+ c_2 d_2) x_1] = \gamma_2\\ x_2=x_2(x_1,e_2)= \alpha_2 - \beta_2 x_1 - \frac{\gamma_2}{1 +d_2+ (e_2 +c_2+ c_2 d_2) x_1}. \end{gathered} \end{equation} Thus, $P_2(0,(1-m_2)/(1+d_2)) \in L_2$ and $L_2$ is a parabola and convex upward, which has a vertex $\hat{P}(\hat{x}_1,\hat{x}_2)$. The asymptotes of $L_2$ are $$ L_{21}: \alpha_2 -x_2- \beta_2 x_1=0,\quad L_{22}: x_1 = - \frac{1+d_2}{e_2 +c_2+ c_2 d_2}. $$ From equations \eqref{eq21}-\eqref{eq22} we have \begin{equation} \label{eq25} \begin{gathered} \frac{\partial x_2}{\partial e_2} >0\quad \text{as } 1-x_1-x_2>0 \\ \lim_{e_2 \to +\infty}x_2 (x_1,e_2) = 1-x_1,\quad \lim_{e_2 \to 0} x_2 (x_1,e_2) = (1-m_2-x_1)/(1+d_2). \end{gathered} \end{equation} When $1-m_2>0$, then $L_2 \cap\operatorname{int}R_+^2 \ne \emptyset$ since $P_2 \in L_2$ and $L_2$ has an asymptote $\alpha_2 -x_2- \beta_2 x_1=0$. When $1-m_2\le 0$, we have $L_2 \cap\operatorname{int}R_+^2 \ne \emptyset$ if and only if $e_2>e_2^{(1)}$. \begin{lemma} \label{le21} The parabolas $L_1$ and $L_2$ have at most four intersection points on the plane, while system \eqref{eq21} has at most two positive equilibria $P^+(x_1^+,x_2^+)$ and $P^-(x_1^-,x_2^-)$ with $x_1^+ >x_1^-$. \end{lemma} \begin{proof} Let $P(x_1,x_2)$ be the intersection point of $L_1$ and $L_2$. From \eqref{eq22} and \eqref{eq24}, a direct computation shows that $P(x_1,x_2)$ satisfies $$ a_0 x_1^4 +a_1 x_1^3 +a_2 x_1^2 +a_3 x_1 +a_4=0,\quad a_0= \beta_2 (1- \beta_1\beta_2)(e_1 +c_1+ c_1 d_1)(e_2 +c_2+ c_2 d_2)^2. $$ Thus, the parabolas $L_1$ and $L_2$ have at most four intersection points on the plane. Assume asymptotes $L_{11}$ and $L_{21}$ coincide. From \eqref{eq22} and \eqref{eq24}, $x_1$ is a linear function of $x_2$. By \eqref{eq22}, $x_1$ satisfies a quadratic equation, which implies that system \eqref{eq21} has at most two positive equilibria. Assume asymptotes $L_{11}$ and $L_{21}$ do not coincide but are parallel. Then we have $1- \beta_1\beta_2=0$ and $a_0=0$. Thus, $x_1$ satisfies a cubic equation, which implies that the parabolas $L_1$ and $L_2$ have at most three intersection points on the plane. If $L_{11}$ is above $L_{21}$. then $L_1$ and $L_2$ have an intersection point in the second quadrant. If $L_{11}$ is below $L_{21}$. then $L_1$ and $L_2$ have an intersection point in the fourth quadrant. Thus, system \eqref{eq21} has at most two positive equilibria. Assume asymptotes $L_{11}$ and $L_{21}$ intersect. Without loss of generality, we suppose $\beta_2 <1/\beta_1$ as shown in Figure \ref{fig1}. Let $L_1^+$ be the part of $L_1$ above $L_{11}$. Let $L_2^-$ be the part of $L_2$ below $L_{21}$. Then $L_1^+$ divide the fourth quadrant into two parts. If $L_2^-$ starts at a point above $L_1^+$, then $L_2^-$ and $L_1^+$ have an intersection point in the fourth quadrant since they have asymptotes $L_{21}$ and $L_{12}$, respectively. If $L_2^-$ starts at a point below $L_1^+$, then $L_2^-$ and $L_1^+$ have an intersection point in the fourth quadrant since they have asymptotes $L_{22}$ and $L_{11}$, respectively. Thus, $L_1$ and $L_2$ always have an intersection point in the fourth quadrant. A similar discussion could show that $L_1$ and $L_2$ always have an intersection point in the second quadrant. Therefore, system \eqref{eq21} has at most two positive equilibria. \end{proof} \begin{lemma} \label{le22} There is no periodic orbit of system \eqref{eq21} and solutions of \eqref{eq21} are bounded. \end{lemma} \begin{proof} Let $f_i(x_1,x_2)$ be the right-hand sides of the equalities in \eqref{eq21}, respectively. Denote $g(x_1,x_2) = 1/x_1x_2$. Then we obtain $$ \frac{\partial (gf_1)}{\partial x_1}+\frac{\partial (gf_2)}{\partial x_2} = - \frac{\gamma_1}{x_2}(1+\frac{e_1 x_2}{1+c_1 x_2}) - \frac{\gamma_2}{x_1}(1+\frac{e_2 x_1}{1+c_2 x_1})<0. $$ By Bendixson-Dulac Theorem \cite{Hofb98}, there is no periodic orbit of \eqref{eq21}. If $x_1 + x_2 \ge 1$, then $dx_1/dt <0$ and $dx_2/dt <0$ by \eqref{eq21}. Thus, all solutions of \eqref{eq21} satisfy $x_1 + x_2 \le 1$ as $t$ is sufficiently large, which implies that solutions of \eqref{eq21} are bounded. \end{proof} The equilibria of \eqref{eq21} on axes are analyzed as follows, while their local stability is determined by eigenvalues of Jacobian matrix of \eqref{eq21} at these equilibria. (a) The trivial equilibrium $O(0,0)$ always exists and has eigenvalues $r_1(1-m_1)$, $r_2(1-m_2)$. (b) The semi-trivial equilibrium $P_1((1-m_1)/(1+d_1),0)$ exists if $1-m_1 > 0$, while $P_2(0, (1-m_2)/(1+d_2))$ exists if $1-m_2 > 0$. The eigenvalues of $P_i$ are $- r_i(1-m_i), -r_j D_j$ with $$ D_j = m_j - \frac{d_i+m_i}{1+d_i} \big[1+ \frac{e_j(1-m_i)}{1+d_i +c_j(1-m_i)}\big], \quad i,j=1,2,\; i \neq j. $$ Stability analysis of system \eqref{eq21} is considered in three situations: (i) obligate mutualisms, i.e., $1-m_1 \le 0, 1-m_2 \le 0$; (ii) obligate-facultative mutualisms, i.e., $1-m_1 > 0, 1-m_2 \le 0$; (iii) facultative mutualisms, i.e., $1-m_1 > 0, 1-m_2 > 0$. %% \input 3-Dyna5 Section 3 \section{Obligate mutualisms} In this section, we consider the situation $1-m_1 \le 0, 1-m_2 \le 0$, which means that neither species can survive in the absence of the other. Denote $$ k_i= \frac{\gamma_i(e_i +c_i+ c_i d_i)}{ (1+ d_i)^2} -\beta_i,\quad i=1,2. $$ \begin{figure}[htb] \begin{center} \includegraphics[width = .9\textwidth]{fig2} \end{center} \caption{ Population dynamics of \eqref{eq21}. Red and green curves are isoclines $L_1$ and $L_2$, and black lines are separatrices of saddle points. Grey arrows denote the direction and strength of the vector fields on the portrait, while filled and open circles represent the stable and unstable equilibria, respectively. (a) Let $m_1= m_2=1, e_1=e_2=10, $ $c_1=c_2=0.001, d_1=d_2=0.8$. Then system \eqref{eq21} has a unique positive equilibrium, which is globally asymptotically stable. (b-c) Fix $m_1= m_2=1.2, $ $c_1=c_2=0.001, d_1=d_2=0.8$ and let $e_1, e_2$ vary. When $e_1=e_2=10$ in (b), there are two positive equilibria $P^-$ and $P^+$. $P^-$ is unstable while $P^+$ is asymptotically stable. When $e_1=e_2=4.75$ in (c), $P^-$ and $P^+$ coincide and form a saddle-node point $P^+$. Separatrices of the saddle points subdivide the first quadrant into two regions, which are basins of attraction of $O$ and $P^+$ respectively. When $e_1=e_2=4.5$ in (d), all solutions converge to $O$.} \label{fig2} \end{figure} \begin{theorem} \label{th31} Assume $1-m_1 = 0,1-m_2 = 0$. If $e_1>1, e_2>1$ and $k_1 k_2>1$, then system \eqref{eq21} has a unique positive equilibrium $P^+$, which is globally asymptotically stable as shown in Figure \ref{fig2}a. Otherwise, all positive solutions of \eqref{eq21} converge to $O$. \end{theorem} \begin{proof} Since $1-m_1 = 1-m_2 = 0$, we obtain $e_1^{(1)}=e_2^{(1)}=1$ and $O \in L_1 \cap L_2$. A direct computation shows that $1/k_1$ and $k_2$ are slopes of $L_1$ and $L_2$ at $O$, respectively. From $e_i>1$, we have $k_i>0$ and $L_i \cap$ int$R_+^2 \ne \emptyset, i=1,2$. If $1/k_1 < k_2$, it follows from the convexity of $L_1$ and $L_2$ that they have a unique intersection point $P^+$ in the first quadrant. By Lemma \ref{le22}, $P^+$ is globally asymptotically stable. In other situations, $L_1$ and $L_2$ have no intersection point in the first quadrant. Thus, system \eqref{eq21} has no positive equilibrium, which implies that all positive solutions of \eqref{eq21} converge to $O$ by Lemma \ref{le22}. \end{proof} In the following result, we focus on the case of $1- m_1 \le 0, 1- m_2 < 0$, while a similar one can be given for the case of $1- m_1 < 0, 1- m_2 \le 0$. \begin{theorem} \label{th32} Assume $1-m_1 \le 0,1-m_2 < 0$. {\rm (i)} Let $e_i > e_i^{(1)}, i=1,2$. There exist $e_1^{(2)}$ and $e_2^{(2)}$ such that when $e_1> e_1^{(2)}$ or $e_2> e_2^{(2)}$, there are two positive equilibria $P^-$ and $P^+$ of \eqref{eq21} as shown in Figure \ref{fig2}b. $P^-$ is a saddle point while $P^+$ is asymptotically stable. When $e_1= e_1^{(2)}$ or $e_2= e_2^{(2)}$, $P^-$ and $P^+$ coincide and form a saddle-node point. Separatrices of the saddle point subdivide the first quadrant into two regions, one is the basin of attraction of $O$ while the other is that of $P^+$, as shown in Figure \ref{fig2}b-c. In other situations, all positive solutions of \eqref{eq21} converge to $O$ as shown in Figure \ref{fig2}d. {\rm (ii)} If $e_1 \le e_1^{(1)}$ or $e_2 \le e_2^{(1)}$, then equilibrium $O$ is globally asymptotically stable. \end{theorem} \begin{proof} (i) When $e_i> e_i^{(1)}$, we obtain $L_i \cap$ int$R_+^2 \ne \emptyset, i=1,2$. It follows from the convexity of $L_2$ and $\lim_{e_2 \to +\infty}x_2 (x_1,e_2) = 1-x_1$ in \eqref{eq25} that there is $e_2^{(2)}>0$ such that when $e_2> e_2^{(2)}$, $L_1$ and $L_2$ have two intersection points $P^-$ and $P^+$; when $e_2= e_2^{(2)}$, $P^-$ and $P^+$ coincide and the isoclines are tangent. Thus, $P^-$ and $P^+$ are positive equilibria of \eqref{eq21} if $e_1> e_1^{(1)},e_2> e_2^{(1)}$ and $e_2> e_2^{(2)}$. The local stability of $P^+$ can be shown as follows. Let $r_ix_iF_i$ denote the righthand sides of \eqref{eq21}. Let $k_1^+$ (resp. $k_2^+$) denote the slope of $L_1$ (resp. $L_2$) at $P^+$. From the expression of $F_i$, we have $$ \frac{1}{k_1^+} = -\frac{\partial F_1}{\partial x_2} \big/ \frac{\partial F_1}{\partial x_1}|_{P^+},\quad k_2^+ = -\frac{\partial F_2}{\partial x_1} \big/ \frac{\partial F_2}{\partial x_2} |_{P^+} $$ where \begin{gather*} \frac{\partial F_1}{\partial x_1}|_{P^+} = -x_1[1+\frac{e_1 x_2}{1+c_1 x_2}+d_1]|_{P^+}<0,\\ \frac{\partial F_2}{\partial x_2}|_{P^+} = -x_2[1+\frac{e_2 x_1}{1+c_2 x_1}+d_2]|_{P^+}<0. \end{gather*} The Jacobian matrix of \eqref{eq21} at $P^+$ is \begin{equation} \label{eq42} J(P^+) = \begin{pmatrix} r_1 x_1 \frac{\partial F_1}{\partial x_1} &r_1 x_1 \frac{\partial F_1}{\partial x_2} \\ r_2 x_2 \frac{\partial F_2}{\partial x_1} &r_2 x_2 \frac{\partial F_2}{\partial x_2} \end{pmatrix}_{P^+}. \end{equation} Thus, $$ \operatorname{tr} J(P^+) = (r_1 x_1 \frac{\partial F_1}{\partial x_1}+ r_2 x_2 \frac{\partial F_2}{\partial x_2})|_{P^+} < 0. $$ Moreover, we have $$ \det J(P^+) = r_1 r_2 x_1 x_2 \big[\frac{\partial F_1}{\partial x_1} \frac{\partial F_2}{\partial x_2} - \frac{\partial F_1}{\partial x_2} \frac{\partial F_2}{\partial x_1} \big] = r_1 r_2 x_1 x_2 \frac{\partial F_1}{\partial x_1} \frac{\partial F_2}{\partial x_2}\big|_{P^+} (1 - \frac{k_2^+}{k_1^+}). $$ Assume $k_1^+ <0$ as shown in Figure \ref{fig2}b. Then $P^+$ is above the vertex of $L_1$. Notice that $L_1$ and the positive $x_2$-axis form a convex area $\Omega_1$. At the point $P^+$, $L_2$ intersects with $L_1$ from the inside of $\Omega_1$ to its outside. Thus, if $k_2^+ <0$, then we have $k_1^+ < k_2^+$ as shown in Figure \ref{fig2}b, which implies that $\det J(P^+)>0$. If $k_2^+ >0$, we can see that $\det J(P^+)>0$. Assume $k_1^+ >0$. Then $P^+$ is below the vertex of $L_1$. When $k_2^+ >0$, we have $k_1^+ >k_2^+ $, which implies that $\det J(P^+)>0$. When $k_2^+ \le 0$, we can see that $\det J(P^+)>0$. It follows from convexity of $L_1$ that $k_1^+ \ne 0$. Thus we have $\det J(P^+)>0$, which implies that $P^+$ is asymptotically stable. Similarly, let $k_1^-$ (resp. $k_2^-$) denote the slope of $L_1$ (resp. $L_2$) at $P^-$. Then we have $k_1^- < k_2^-$ and $\det J(P^-)<0$, which implies that $P^-$ is a saddle point. By Lemma \ref{le22}, the omega limit set of any interior point is an equilibrium since system \eqref{eq21} is analytic and has no graphic here. Therefore, the stable manifold of $P^-$ subdivide the first quadrant into two regions, one is the basin of attraction of $O$ while the other is that of $P^+$. When $e_2= e_2^{(2)}$, $P^-$ and $P^+$ coincide and form a saddle-node point \cite{Perk01}. Thus, separatrices of the saddle point subdivide the first quadrant into two regions, one is the basin of attraction of $O$ while the other is that of $P^+$. A similar discussion can be given for $L_1$ and $e_1^{(2)}$. {\rm (ii)} If $e_1 \le e_1^{(1)}$ or $e_2 \le e_2^{(1)}$, there is no positive equilibrium of \eqref{eq21}, which implies that all solutions of \eqref{eq21} converge to $O$. \end{proof} Theorems \ref{th31}-\ref{th32} demonstrate essential features of obligate mutualisms. \textit{First,} mutualisms between species can lead to survival of obligate species. Recall that the mutualistic effect consists of two factors, e.g., $e_1 = \mu_1 \nu_2$, where $\nu_2$ denotes the quantity of resources that an individual of species 2 produces for species 1, and $\mu_1$ is the efficiency of species 1 in converting the resources into fitness. In this section, we focus on the role of resources such as food, while we focus on the role of efficiency in section 6. In the situation considered by Theorem \ref{th32}, neither species can survive in the absence of the other. Theorem \ref{th32}(i) exhibits that they can coexist if (a) mutualistic effects between them are above a threshold ($e_i >e_i^{(1)}$), (b) one of the effects is sufficiently strong (e.g., $e_2 > e_2^{(2)}$), and (c) population densities of the species are large. The reason is that under these conditions, each species can produce abundant food for the other, which results in persistence of both species. Ecologically, this result can provide an explanation for the reason why obligate bacteria can coexist by consuming products of the other. When one of the above conditions is not satisfies, at least one of the species cannot obtain sufficient food from the other, which eventually leads to extinction of both species. Since neither species can survive alone, it is the mutualism that leads to their persistence. A similar discussion can be given for Theorem \ref{th31}. \textit{Second,} intermediate mutualism is beneficial in certain parameter ranges. In the situation considered by Theorems \ref{th31}-\ref{th32}, we have $1- m_i \le 0, i=1,2$. $L_1 \cap\operatorname{int}R_+^2 \ne \emptyset$ if $e_1>e_1^{(1)}$. It follows from the convexity of $L_1$ that its vertex $\bar{P}(\bar{x}_1,\bar{x}_2)$ is in the first quadrant. For a fix $e_1(>e_1^{(1)})$, it follows from the monotonicity of $L_2$ that there is $\bar{e}_2>0$ such that when $e_2=\bar{e}_2$, we have $P^+=\bar{P}$ and $x_1^+=\bar{x}_1$. Thus, when $e_2=\bar{e}_2$, species 1 approaches a maximal density $\bar{x}_1$, as shown in Figure \ref{fig2}b. An over-mutualism ($e_2>\bar{e}_2$) or under-mutualism ($e_2<\bar{e}_2$) means that $x_1^+<\bar{x}_1$, which is not the best for species 1. The reason is as follows: (I) When $e_2>\bar{e}_2$, the over-mutualism leads to the increase of species 2 who will occupy more sites and will decrease the density of species 1; (II) When $e_2<\bar{e}_2$, the under-mutualism leads to the decrease of species 2. Thus species 2 cannot produce abundant food for species 1 to approach its maximum. Therefore, \textit{only} the intermediate mutualistic effect $e_2=\bar{e}_2$ is the best for species 1. A similar discussion can be given for species 2. \textit{Third,} extremely strong mutualism will lead to extinction of both species. Indeed, Theorem \ref{th32}(i) exhibits that there is a stable positive equilibrium $P^+(x_1^+, x_2^+)$ if mutualistic effects ($e_i$) are large. From \eqref{eq25}, we obtain $\lim_{e_2 \to +\infty}x_2 (x_1,e_2) = 1-x_1$. Since $P_2 \in L_2$, it follows from the convexity of $L_1$ and $L_2$ that when $e_2 \to +\infty$, $P^+$ tends to the $x_2$-axis with $x_1^+ \to 0$. That is, $\lim_{e_2 \to +\infty} x_1^+ = 0 $. Similarly, we have $\lim_{e_1 \to +\infty} x_2^+ = 0 $. Therefore, we conclude the following result. \begin{lemma} \label{le41} If $P^+(x_1^+, x_2^+)$ is a positive equilibrium of \eqref{eq21}, then $$ \lim_{e_2 \to +\infty} x_1^+ = 0,\quad \lim_{e_1 \to +\infty} x_2^+ = 0. $$ \end{lemma} When the mutualistic effect of species 1 on 2 is extremely strong ($e_2 \to +\infty$), Lemma \ref{le41} shows that species 1 goes to extinction, which implies the extinction of species 2 who cannot survive alone. The reason is that the extremely strong mutualism leads to an explosive growth of species 2 who will occupy most of the sites and drive species 1 into extinction. Thus, both species will go to extinction since neither species can survive alone. Therefore, extremely strong mutualism will lead to extinction of both species. \textit{Finally,} extremely weak mutualism will result in extinction of both species. Indeed, when the mutualism is extremely weak (e.g., $e_2 \to 0$), Theorems \ref{th31}-\ref{th32} demonstrate that both species go to extinction. This is because one of the species will get little food from the other and will go to extinction, which implies extinction of both species since neither one can survive alone. Thus, extremely weak mutualism will result in extinction of both species. \section{Obligate-facultative mutualisms} In this section, we consider the situation where one species can persist in the absence of the other but the other cannot survive alone. We focus on the case of $1-m_1 > 0, 1-m_2 \le 0$, while a similar discussion can be given for $1-m_1 \le 0,1-m_2 > 0$. When $m_1\ne 1, m_2\ne 1$, we denote \begin{equation} \label{eq31} \begin{gathered} e_1^{(3)} =(c_1+\frac{1+d_2 }{1-m_2}) ( \frac{m_1+d_2m_1}{d_2+m_2}-1),\\ e_2^{(3)} =(c_2+\frac{1+d_1 }{1-m_1}) ( \frac{m_2+d_1m_2}{d_1+m_1}-1). \end{gathered} \end{equation} It follows from \eqref{eq22} and \eqref{eq24} that if $e_1 > e_1^{(3)}$, $P_2$ is at the right of $L_1$; if $e_2 > e_2^{(3)}$, $P_1$ is below $L_2$. \begin{theorem} \label{th33} Assume $1-m_1 > 0, 1-m_2 = 0$. If $e_2 > e_2^{(3)}$, system \eqref{eq21} has a unique positive equilibrium $P^+$, which is globally asymptotically stable. Otherwise, all positive solutions of \eqref{eq21} converge to $P_1$. \end{theorem} \begin{proof} From $m_2 = 1$, we have $O \in L_2$. If $e_2 > e_2^{(3)}$, then $P_1$ is below $L_2$ and is a saddle point. It follows from the convexity of $L_1$ and $L_2$ that they have a unique intersection point $P^+$ in the first quadrant. Thus system \eqref{eq21} has a unique positive equilibrium $P^+$. By Lemma \ref{le22}, $P^+$ is globally asymptotically stable. If $e_2 \le e_2^{(3)}$, system \eqref{eq21} has no positive equilibrium, which implies that all positive solutions of \eqref{eq21} converge to $P_1$. \end{proof} Assume $1-m_1 > 0, 1-m_2 < 0$. When $e_2 > e_2^{(3)}$, equilibrium $P_1$ is below $L_2$ and is a saddle point. By a proof similar to that of Theorem \ref{th33}, system \eqref{eq21} has a unique positive equilibrium $P^+$, which is globally asymptotically stable. Suppose $e_2 < e_2^{(3)}$. Then $P_1$ is above $L_2$. When $e_1 \le e_1^{(3)}$, $P_2$ is at the left of $L_1$ and system \eqref{eq21} has no positive equilibrium. Thus $P_1$ is globally asymptotically stable. Let $e_1 > e_1^{(3)}$. It follows from the convexity of $L_2$ and $\lim_{e_2 \to +\infty}x_2 (x_1,e_2) = 1-x_1$ in \eqref{eq25} that there is $e_2^{(4)}>0$ such that when $e_2 > e_2^{(4)}$, there is at least one intersection point of $L_1$ and $L_2$ in the first quadrant. Thus, when $e_2 > \max\{e_2^{(2)},e_2^{(4)} \}$, $L_1$ and $L_2$ have two positive intersection points $P^-$ and $P^+$. By Lemma \ref{le22}, phase-portrait analysis shows that equilibrium $P^-$ is a saddle point and $P^+$ is asymptotically stable. When $e_2 =e_2^{(2)}>e_2^{(4)}$, $P^-$ and $P^+$ coincide and form a saddle-node point. In other situations, there is no positive equilibrium of \eqref{eq21}. \begin{figure}[htb] \begin{center} \includegraphics[width = 0.9\textwidth]{fig3} \end{center} \caption{ (a) Population dynamics of \eqref{eq21}. Let $m_1=0.5, m_2=1.5$, $c_i=0.001, d_i=0.8, e_i=8, i=1,2$. Then system \eqref{eq21} has a unique positive equilibrium, which is globally asymptotically stable. (b) Fix $m_1=0.95, m_2=1.15, e_2=5, $ $c_1=c_2=0.001, d_1=d_2=0.8$ and Let $e_1$ vary. When $e_1$ increases from 2, 3.5 to 6.5, isocline $L_1$ increases monotonically while species 2 could persist at $e_1 \ge 3.5$. (c) Fix $m_1=0.95, m_2=1.15, e_1=4.5, e_2=5, $ $c_1=0.001, d_1=0.8, d_2=0.02$ and let $c_2$ vary. When $c_2$ decreases from 5.5, 2.5 to 0.02, isocline $L_2$ increases monotonically while species 2 could persist at $c_2 =2.5, 0.02$. (d) Fix $m_1=0.95, m_2=1.15, e_1=4.5, e_2=5, $ $c_1=c_2=0.001, d_1=0.8$ and let $d_2$ vary. When $d_2$ decreases from 5.8, 0.9 to 0.02, isocline $L_2$ increases monotonically while species 2 could persist at $d_2 = 0.9, 0.02$. } \label{fig3} \end{figure} Suppose $e_2 = e_2^{(3)}$. Then $P_1 \in l_2$. When $e_1 > e_1^{(2)}$, there is a unique positive equilibrium $P^+$, which is globally asymptotically stable. Otherwise, there is no positive equilibrium and all positive solutions of \eqref{eq21} converge to $P_1$. Therefore, we conclude the following result. \begin{theorem} \label{th34} Assume $1-m_1 > 0$, $1-m_2 < 0$. {\rm (i)} If $e_2 > e_2^{(3)}$, system \eqref{eq21} has a unique positive equilibrium $P^+$, which is globally asymptotically stable as shown in Figure \ref{fig3}a. {\rm (ii)} If $e_2 < e_2^{(3)}$ and $e_1 \le e_1^{(3)}$, then all positive solutions of \eqref{eq21} converge to $P_1$. {\rm (iii)} Let $e_2 < e_2^{(3)}$ and $e_1 > e_1^{(3)}$. \begin{itemize} \item[(a)] If $e_2 > \max\{e_2^{(2)},e_2^{(4)} \}$, system \eqref{eq21} has two positive equilibria $P^-$ and $P^+$. $P^-$ is a saddle point and $P^+$ is asymptotically stable. If $e_2 =e_2^{(2)}>e_2^{(4)}$, $P^-$ and $P^+$ coincide and form a saddle-node point. The separatrices of the saddle point subdivide the first quadrant into two regions, one is the basin of attraction of $P_1$ while the other is that of $P^+$, as shown in Figure \ref{fig3}b. \item[(b)] In other situations, system \eqref{eq21} has no positive equilibrium and $P_1$ is globally asymptotically stableas shown in Figure \ref{fig3}b. \end{itemize} {\rm (iv)} Let $e_2 =e_2^{(3)}$. If $e_1 > e_1^{(2)}$, system \eqref{eq21} has a unique positive equilibrium $P^+$, which is globally asymptotically stable. Otherwise, all positive solutions of \eqref{eq21} converge to $P_1$. \end{theorem} Theorems \ref{th33}-\ref{th34} demonstrate essential features of obligate-facultative mutualisms. \textit{First,} an obligate species can survive by cooperating with a facultative one. In the situation considered by Theorem \ref{th34}, species 1 can persist in the absence of species 2, while species 2 cannot survive alone. If the mutualistic effect of species 1 on 2 is strong ($e_2 > e_2^{(3)}$), Theorem \ref{th34}(i) exhibits that species 2 can survive. The reason is that in this case, species 1 can provide abundant food for species 2, which leads to its survival. Moreover, even when the mutualistic effect of species 1 on 2 is intermediate ($e_2^{(4)} < e_2 < e_2^{(3)}$), Theorem \ref{th34}(iia) demonstrates that species 2 can survive if the effect of species 2 on 1 is strong ($e_1 > \max\{e_1^{(2)}, e_1^{(3)} \}$) and their population densities are large (Here, $e_1 > e_1^{(2)}$ is equivalent to $e_2 > e_2^{(2)}$ in guaranteeing that $L_1$ and $L_2$ can intersect in the first quadrant). The reason is that in this case, species 2 can provide abundant food for species 1, which leads to the increase of species 1. Thus, the amount of food produced by species 1 for 2 is enhanced, which leads to the survival of species 2 in return. Hence, Theorem \ref{th34}(iia) exhibits a strategy for obligate species when cooperating with facultative ones: if the mutualistic effect from facultative species is intermediate, the obligate species can survive by strengthening its mutualistic effect on the facultative species (e.g., producing more food or providing better service for the obligate species) and enhancing its population density. A similar discussion can be given for Theorem \ref{th33}. \textit{Second,} intermediate mutualism is beneficial in certain parameter ranges. By \eqref{eq25}, there exists $\check{e}_1$ such that when $e_1>\check{e}_1$, the vertex $\bar{P}(\bar{x}_1,\bar{x}_2)$ of $L_1$ is in the first quadrant, as shown in Figure \ref{fig3}a. Similar to the discussion for obligate mutualisms, there is $\bar{e}_2>0$ such that when $e_2=\bar{e}_2$, species 1 approaches a maximal density $\bar{x}_1$, as shown in Figure \ref{fig3}a. An over-mutualism ($e_2>\bar{e}_2$) or under-mutualism ($e_2<\bar{e}_2$) means that $x_1^+<\bar{x}_1$, which is not the best for species 1. Therefore, \textit{only} the intermediate mutualistic effect $e_2=\bar{e}_2$ is the best for species 1. A similar discussion can be given for species 2. \textit{Third,} extremely strong mutualism will result in extinction of one/both species. (a) If the mutualistic effect of species 1 on 2 is extremely strong, Lemma \ref{le41} exhibits that species 1 goes to extinction, which implies extinction of species 2 who cannot survive alone. Thus, extremely strong mutualism from facultative species to obligate species will result in extinction of both species. (b) If the mutualistic effect of species 2 on 1 is extremely strong, Lemma \ref{le41} exhibits that species 2 goes to extinction, while species 1 approaches its carrying capacity. Thus, extremely strong mutualism from obligate species to facultative species will result in extinction of the obligate species itself. \textit{Finally,} extremely weak mutualism can result in extinction of obligate species. (a) If the mutualistic effect of species 1 on 2 is extremely weak, Theorems \ref{th33}-\ref{th34} show that species 2 goes to extinction because it depends upon the strong mutualism of species 1 for survival. (b) If the mutualistic effect of species 2 on 1 is extremely weak, species 2 goes to extinction when the mutualism from species 1 is not strong. An interesting phenomenon is: when the mutualism from species 1 on 2 is intermediate, Theorem \ref{th34}(iia) shows that species 2 can survive if its mutualism on species 1 is strong. However, if its mutualism on species 1 is weak, species 2 goes to extinction. What a pity! \section{Facultative mutualisms} \begin{figure}[htb] \begin{center} \includegraphics[width = 0.6\textwidth]{fig4} \end{center} \caption{ Fix $m_1=0.2, m_2=0.8, e_1=5, $ $c_1=c_2=0.001, d_1=d_2=0.8$ and let $e_2$ vary. When $e_2=0.8$, species 2 goes to extinction. When $e_2=5$, the two species coexist and form a win-win situation. When $e_2 =500 ($ i.e. $ e_2 \to +\infty)$, species 1 goes to extinction.} \label{fig4} \end{figure} In this section, we consider the situation of $1-m_1 > 0, 1-m_2 > 0$, where either species can survive in the absence of the other. \begin{theorem} \label{th35} Assume $1-m_1 > 0$, $1-m_2 > 0$. {\rm (i)} If $e_1 > e_1^{(3)}$ and $e_2 > e_2^{(3)}$, system \eqref{eq21} has a unique positive equilibrium $P^+$, which is globally asymptotically stable as shown in Figure \ref{fig4}. {\rm (ii)} If $e_1 \le e_1^{(3)}$, then $P_2$ is globally asymptotically stable. If $e_2 \le e_2^{(3)}$, then $P_1$ is globally asymptotically stable as shown in Figure \ref{fig4}. \end{theorem} \begin{proof} (i) Since $e_1 > e_1^{(3)}$, $P_2$ is at the left of $L_1$ and is a saddle point. Since $e_2 > e_2^{(3)}$, $P_1$ is below $L_2$ and is a saddle point. Thus, $L_1$ and $L_2$ have intersection points in the first quadrant. Moreover, since $L_1$ and $L_2$ have asymptotes $x_2 = - \frac{1+d_1}{e_1 +c_1+ c_1 d_1}$, $x_1 = - \frac{1+d_2}{e_2 +c_2+ c_2 d_2}$ respectively, they have an intersection point in the third quadrant. Thus, $L_1$ and $L_2$ have a unique intersection point in the first quadrant by Lemma \ref{le21}. Thus system \eqref{eq21} has a unique positive equilibrium $P^+$. By Lemma \ref{le22}, $P^+$ is globally asymptotically stable. (ii) Assume $e_2 \le e_2^{(3)}$. Then $P_1$ is above $L_2$ and is asymptotically stable. Thus $L_1$ and $L_2$ have an intersection point in the fourth quadrant. Since they have an intersection point in the third quadrant as shown in (i), $L_1$ and $L_2$ have no intersection point in the first quadrant by Lemma \ref{le21}. Thus system \eqref{eq21} has no positive equilibrium. By Lemma \ref{le22}, $P_1$ is globally asymptotically stable. A similar discussion can be given for the case $e_1 \le e_1^{(3)}$. \end{proof} Theorem \ref{th35} demonstrates essential features of facultative mutualisms. \textit{First,} mutualisms can lead to interaction outcomes $(+~+)$, where each species can approach a density larger than its carrying capacity in the absence of the other. Indeed, when $(1-m_1)/(1+d_1) + (1-m_2)/(1+d_2) <1$, it follows from the monotonicity of $L_i$ with $e_i$ that there is a region $R_{12}$ such that when $(e_1,e_2) \in R_{12}$, $P^+$ is at the right of $P_1$ and above $P_2$, which implies that either species approaches a density larger than its carrying capacity in the absence of the other. Thus the interaction outcomes are $(+~+)$. \textit{Second,} intermediate mutualism is beneficial in certain parameter ranges. By \eqref{eq25}, there exists $\check{e}_1$ such that when $e_1>\check{e}_1$, the vertex $\bar{P}(\bar{x}_1,\bar{x}_2)$ of $L_1$ is in the first quadrant. Similar to the discussion for obligate mutualisms, there is $\bar{e}_2>0$ such that when $e_2=\bar{e}_2$, species 1 approaches a maximal density $\bar{x}_1$, as shown in Figure \ref{fig4}. An over-mutualism ($e_2>\bar{e}_2$) or under-mutualism ($e_2<\bar{e}_2$) means that $x_1^+<\bar{x}_1$, which is not the best for species 1. Therefore, \textit{only} the intermediate mutualistic effect $e_2=\bar{e}_2$ is the best for species 1. A similar discussion can be given for species 2. \textit{Third,} extremely strong mutualism will result in extinction of species. If the mutualistic effect of species 1 on 2 is extremely strong ($e_2 \to +\infty$), Lemma \ref{le41} exhibits that species 1 goes to extinction while species 2 persists. A similar discussion can be given for $e_1 \to +\infty$. \textit{Finally,} extremely weak mutualism can result in extinction of species. If the mutualistic effect of species 1 on 2 is extremely weak, Theorem \ref{th35}(ii) exhibits that species 2 goes to extinction while species 1 persists. The reason is that species 2 obtains little food from species 1, which leads to its failure in spatial competition with species 1. A similar discussion can be given for species 1. \section{Discussion} In this article, we extend a lattice gas model of mutualisms in \cite{Iwat11} by considering saturated response and self-competition. Global dynamics of the extended model demonstrate some basic properties of mutualisms. \textit{First,} mutualisms (i.e. $e_i$) can lead to survival of mutualists. As discussed in Section 2, the mutualistic effect $e_i$ consists of two factors: one is the quantity of resources provided by collaborators, while the other is the efficiency of the species in converting the resources into fitness. While we focus on resources (food) in Sections 3-5, we consider the efficiency in this section. In the situation considered by Theorem \ref{th34}, species 1 can persist in the absence of species 2 while species 2 cannot survive alone. Theorem \ref{th34} demonstrate that species 2 can survive if its efficiency in converting the food (provided by species 1) into fitness is high. Even when the efficiency is intermediate, species 2 can survive if it has a large population density and the efficiency of species 1 is high. In numerical simulations of Figure \ref{fig3}b, we fix other parameters and let $e_1$ vary. Here, the mutualistic effect of species 1 on 2 is intermediate. When $e_1$ increases from 2 to 6.5, isocline $L_1$ increases monotonically and species 2 with large initial density can survive if $e_1 \ge 3.5$, which confirms that species 2 with large population density can survive if the efficiency of species 1 is high. A similar discussion can be given for Theorems \ref{th33}, \ref{th51}, \ref{th52}, \ref{th61}. \textit{Second,} intermediate mutualisms are favorable under certain parameter conditions, while extremely weak/strong mutualisms will lead to extinction of species. For example, when the mutualistic effect of species 2 on 1 is strong such that $\bar{x}_2>0$, there exists an interval of $e_2$ (intermediate mutualisms), in which species 1 can approach a density larger than its carrying capacity in the absence of species 2. On the other hand, as discussed in section 6, extremely weak mutualisms lead to little benefit to mutualists, while extremely strong mutualisms results in dramatic growth of one species, which implies extinction of the other. Therefore, extremely weak/strong mutualisms will lead to extinction of species. A similar discussion can be given for Theorems \ref{th51},\ref{th52},\ref{th61}. \textit{Third,} parameters $c_i$ and $d_i$ play an important role in survival of species. By \eqref{eq21}, $e_i/c_i$ represents the saturation level while $1/c_i$ is the half-saturation density. Thus, the decrease of $c_i$ promotes persistence of species $i$ by enlarging the functional response $e_ix_j/(1+c_ix_j)$ in \eqref{eq21}. Similarly, $d_i$ represents the degree of intraspecific competition in species $i$. The decrease of $d_i$ promotes the growth of species $i$ by enlarging the function $x_i=x_i(x_j, d_i)$ in \eqref{eq210}. In Figure \ref{fig3}c, we fix other parameters but let $c_2$ vary. When $c_2$ decreases from 5.5 to 0.02, isocline $L_2$ increases monotonically while species 2 could persist at $c_2 =0.02$. In Figure \ref{fig3}d, we fix $m_1=0.95$, $m_2=1.15$, $e_1=4.5$, $e_2=5$, $c_1=c_2=0.001$, $d_1=0.8$ and let $d_2$ vary. When $d_2$ decreases from 5.8 to 0.02, isocline $L_2$ increases monotonically while species 2 could persist at $d_2 =0.02$. Therefore, the decrease of $c_i$ and/or $d_i$ promotes persistence of the species. \textit{Finally,} population densities are crucial to survival of species. When neither species can survive in the absence of the other and the mutualistic effects are strong, Theorem \ref{th32} exhibits that the species can persist only if their initial densities are large. Otherwise, both species go to extinction. Therefore, Theorem \ref{th32} predicts the Allee effect in obligate mutualisms of two species. Similar discussions can be given for cases considered by Theorems \ref{th31},\ref{th33}-\ref{th34}. Iwata et al \cite{Iwat11} displayed seven typical types of dynamics of a lattice gas model in their Figure \ref{fig4} and section 5.2. However, they did not obviously exhibit the type of coexistence in Theorems \ref{th34}(i) as shown in Figure \ref{fig3}a. Indeed, if $c_i=d_i=0$, $i=1,2$, then system \eqref{eq21} becomes the model in \cite{Iwat11}. Let $c_i=d_i=0$, $i=1,2$, $m_1=0.5$, $m_2=1.5$, $e_1=e_2=8$. Numerical simulations could show that species 2 could survive by the mutualism of species 1. Thus Theorem \ref{th34}(i) extends the result by Iwata et al \cite{Iwat11}. Although there is no real data to verify dynamics of the model, the model demonstrates mechanisms which seems to be consistent with ecological situations. For example, as shown in Theorems \ref{th33}-\ref{th34} and Figure \ref{fig3}, different mutualistic effects (or exclusive competitions) can lead to different components of syntrophic colonies, which is crucial to the development of colony architectures. Thus, this model may be useful in the study of cooperative association like bacterial species. Although the model is simple, its global dynamics demonstrate some essential features of mutualisms, which may be helpful for understanding complexity of mutualisms in real situations. \section{Appendix: The case of $r_1>0, r_2 = 0$} In this case, species 1 can reproduce alone but species 2 cannot, while a similar discussion can be given for the case $r_1=0, r_2 > 0$, as shown by Iwata et al \cite{Iwat11}. Denote $$ m_1 := \frac{m_1}{r_1},\quad e_1 := \frac{e_1}{r_1},\quad d_1 := \frac{d_1}{r_1} $$ then model \eqref{eq13} can be rewritten as \begin{equation} \label{eq51} \begin{gathered} \frac{dx_1}{dt} =r_1 x_1 \big[-m_1 + (1+\frac{e_1 x_2}{1+c_1 x_2}) (1-x_1-x_2) -d_1x_1 \big],\\ \frac{dx_2}{dt} =x_2 \big[-m_2 + \frac{e_2 x_1}{1+c_2 x_1} (1-x_1-x_2) -d_2x_2\big]. \end{gathered} \end{equation} The isoclines $L_2$ of \eqref{eq51} can be rewritten as $$ (\bar{\alpha}_1 -x_1- x_2) (d_2+ e_2 x_1) = \bar{\gamma}_1 $$ where $$ \bar{\alpha}_1 = \frac{e_2+ d_2- m_2 c_2}{e_2},\quad \bar{\gamma}_1 = m_1 + \frac{d_2(e_2+ d_2- m_2 c_2)}{e_2}. $$ Thus, $L_2$ is a parabola with asymptotes $$ L_{21}: \bar{\alpha}_1 -x_1- x_2=0,\quad L_{22}: x_1 = - \frac{d_2}{e_2}. $$ Denote $$ \bar{e}_1^{(1)} = e_1^{(1)} , ~~ \bar{e}_2^{(1)} = c_2 m_2 + 2m_2 +2\sqrt{m_2(1+c_2)}. $$ When $e_i>\bar{e}_i^{(1)}$, we have $L_i \cap\operatorname{int}R_+^2 \ne \emptyset$, $i=1,2$. The equilibria of \eqref{eq51} are as follows, while their local stability is determined by eigenvalues of Jacobian matrix of \eqref{eq51} at these equilibria. (a) Equilibrium $O(0,0)$ always exists and has eigenvalues $r_1(1-m_1), -m_2$. (b) Equilibrium $P_1((1-m_1)/(1+d_1),0)$ exists if $1-m_1 > 0$. The eigenvalues of $P_1$ are $- r_1(1-m_1), - \bar{D}_2$ with $$ \bar{D}_2 = m_2 - \frac{e_2 (1-m_1) (d_1+m_1)}{(1+d_1)[1+d_1 +c_2(1-m_1)]}. $$ (c) There are at most two positive equilibria $P^-$ and $P^+$ of \eqref{eq51} by a proof similar to that in section 2. When they exist, $P^-$ is a saddle point and $P^+$ is asymptotically stable. Assume $1- m_1 \le 0$ and $e_i>\bar{e}_i^{(1)}, i=1,2$. Then $L_i \cap$ int$R_+^2 \ne \emptyset, i=1,2$. It follows from the convexity of $L_2$ that there is $\bar{e}_2^{(2)}>0$ such that when $e_2> \bar{e}_2^{(2)}$, $L_1$ and $L_2$ have two intersection points $P^-$ and $P^+$; when $e_2= \bar{e}_2^{(2)}$, $P^-$ and $P^+$ coincide and the isoclines are tangent. Thus, $P^-$ and $P^+$ are positive equilibria of \eqref{eq21} if $e_1> \bar{e}_1^{(1)},e_2> \bar{e}_2^{(1)}$ and $e_2 \ge \bar{e}_2^{(2)}$. A similar discussion can be given for $L_1$ and $\bar{e}_1^{(2)}$. If $e_1< \bar{e}_1^{(2)}$ or $e_2 < \bar{e}_2^{(2)}$, there is no positive equilibrium of \eqref{eq51}, which implies that all solutions of \eqref{eq51} converge to $O$. Therefore, we conclude the following result. \begin{theorem} \label{th51} Assume $1-m_1 \le 0$. {\rm (i)} Let $e_i > \bar{e}_i^{(1)}$, $i=1,2$. If $e_1> \bar{e}_1^{(2)}$ or $e_2> \bar{e}_2^{(2)}$, there are two positive equilibria $P^-$ and $P^+$ of \eqref{eq51}. $P^-$ is a saddle point while $P^+$ is asymptotically stable. If $e_1= \bar{e}_1^{(2)}$ or $e_2= \bar{e}_2^{(2)}$, $P^-$ and $P^+$ coincide and form a saddle-node point. Separatrices of the saddle point subdivide the first quadrant into two regions, one is the basin of attraction of $O$ while the other is that of $P^+$. In other situations, all positive solutions of \eqref{eq51} converge to $O$. {\rm (ii)} If $e_1 \le \bar{e}_1^{(1)}$ or $e_2 \le \bar{e}_2^{(1)}$, then equilibrium $O$ is globally asymptotically stable. \end{theorem} Assume $1-m_1 > 0$. Denote $$ \bar{e}_1^{(3)} = (c_1 - \frac{d_2}{m_2}) ( \frac{m_1d_2}{m_2+d_2} -1 ),\quad \bar{e}_2^{(3)} = \frac{m_2 (1+d_1) [1+d_1 +c_2(1-m_1)] }{(1-m_1) (d_1+m_1)}. $$ When $e_2 > \bar{e}_2^{(3)}$, $P_1$ is below $L_2$ and $\bar{D}_2 <0$; when $e_1 > \bar{e}_1^{(3)}$, $P_2(0,-m_2/d_2)$ is at the right of $L_1$. It follows from the convexity of $L_2$ that there is $\bar{e}_2^{(4)}>0$ such that when $e_2 > \bar{e}_2^{(4)}$, there is at least one positive intersection point of $L_1$ and $L_2$. Thus, when $e_1 > \bar{e}_1^{(3)}$ and $ e_2 > \max\{\bar{e}_2^{(2)},\bar{e}_2^{(4)} \}$, $L_1$ and $L_2$ have two intersection points. Similar to the discussion in section 2, we conclude the following result. \begin{theorem} \label{th52} Assume $1-m_1 > 0$. {\rm (i)} If $e_2 > \bar{e}_2^{(3)}$, system \eqref{eq51} has a unique positive equilibrium $P^+$, which is globally asymptotically stable. {\rm (ii)} If $e_2 < \bar{e}_2^{(3)}$ and $e_1 \le \bar{e}_1^{(3)}$, then all positive solutions of \eqref{eq51} converge to $P_1$. {\rm (iii)} Let $e_2 < \bar{e}_2^{(3)}$ and $e_1 > \bar{e}_1^{(3)}$. \begin{itemize} \item[(a)] If $e_2 > \max\{\bar{e}_2^{(2)},\bar{e}_2^{(4)} \}$, system \eqref{eq51} has two positive equilibria $P^-$ and $P^+$. $P^-$ is a saddle point and $P^+$ is asymptotically stable. If $e_2 =\bar{e}_2^{(2)}>\bar{e}_2^{(4)}$, $P^-$ and $P^+$ coincide and form a saddle-node point. The separatrices of the saddle point subdivide the first quadrant into two regions, one is the basin of attraction of $P_1$ while the other is that of $P^+$. \item[(b)] In other situations, system \eqref{eq51} has no positive equilibrium and $P_1$ is globally asymptotically stable. \end{itemize} {\rm (iv)} Let $e_2 =\bar{e}_2^{(3)}$. If $e_1 > \bar{e}_1^{(2)}$, system \eqref{eq51} has a unique positive equilibrium $P^+$, which is globally asymptotically stable. Otherwise, all positive solutions of \eqref{eq51} converge to $P_1$. \end{theorem} \section{Appendix: The case of $r_1 = r_2 = 0$} In this case, neither species can reproduce alone, which is the same as reactions in a male-female system as shown by Iwata et al \cite{Iwat11}. Let $r_1 = r_2 = 0$. Model \eqref{eq13} becomes \begin{equation} \label{eq61} \begin{gathered} \frac{dx_1}{dt} =x_1 [-m_1 + \frac{e_1 x_2}{1+c_1 x_2} (1-x_1-x_2) -d_1x_1 ]\\ \frac{dx_2}{dt} =x_2 [-m_2 + \frac{e_2 x_1}{1+c_2 x_1} (1-x_1-x_2) -d_2x_2]. \end{gathered} \end{equation} The isoclines $L_i$ of \eqref{eq61} can be rewritten as \begin{gather*} L_1: (\hat{\alpha}_1 -x_1- x_2) (d_1+ e_1 x_2) = \hat{\gamma}_1,\\ L_2: (\hat{\alpha}_2 -x_1- x_2) (d_2+ e_2 x_1) = \hat{\gamma}_2 \end{gather*} where $$ \hat{\alpha}_i = \frac{e_i+ d_i- m_i c_i}{e_i},\quad \hat{\gamma}_i = m_i + \frac{d_i(e_i+ d_i- m_i c_i)}{e_i},~ i=1,2. $$ Thus, $L_i$ is a parabola with asymptotes $$ L_{i1}: \hat{\alpha}_i -x_1- x_2=0,\quad L_{i2}: x_j = - \frac{d_i}{e_i},\quad i\ne j, \; i,j=1,2. $$ Denote $$ \hat{e}_1^{(1)} = c_1 m_1 + 2m_1 +2\sqrt{m_1(1+c_1)} , \quad \hat{e}_2^{(1)} = c_2 m_2 + 2m_2 +2\sqrt{m_2(1+c_2)}. $$ When $e_i>\hat{e}_i^{(1)}$, we have $L_i \cap\operatorname{int}R_+^2 \ne \emptyset, i=1,2$. The equilibria of \eqref{eq61} are as follows, while their local stability is determined by eigenvalues of Jacobian matrix of \eqref{eq61} at these equilibria. (a) Equilibrium $O(0,0)$ always exists and has eigenvalues $-m_1, -m_2$. (b) There are at most two positive equilibria $P^-$ and $P^+$ of \eqref{eq61} by a proof similar to that in section 2. When they exist, $P^-$ is a saddle point and $P^+$ is asymptotically stable. Assume $e_i>\hat{e}_i^{(1)}$, $i=1,2$. Then $L_i \cap\operatorname{int}R_+^2 \ne \emptyset$, $i=1,2$. It follows from the convexity of $L_2$ that there is $\hat{e}_2^{(2)}>0$ such that when $e_2> \hat{e}_2^{(2)}$, $L_1$ and $L_2$ have two intersection points $P^-$ and $P^+$; when $e_2= \hat{e}_2^{(2)}$, $P^-$ and $P^+$ coincide and the isoclines are tangent. Thus, $P^-$ and $P^+$ are positive equilibria of \eqref{eq61} if $e_1> \hat{e}_1^{(1)},e_2> \hat{e}_2^{(1)}$ and $e_2 \ge \hat{e}_2^{(2)}$. A similar discussion can be given for $L_1$ and $\hat{e}_1^{(2)}$. If $e_1< \hat{e}_1^{(2)}$ or $e_2 < \hat{e}_2^{(2)}$, there is no positive equilibrium of \eqref{eq61}, which implies that all solutions of \eqref{eq61} converge to $O$. Therefore, we conclude the following result. \begin{theorem} \label{th61} {\rm (i)} Let $e_i > \hat{e}_i^{(1)}, i=1,2$. If $e_1> \hat{e}_1^{(2)}$ or $e_2> \hat{e}_2^{(2)}$, there are two positive equilibria $P^-$ and $P^+$ of \eqref{eq61}. $P^-$ is a saddle point while $P^+$ is asymptotically stable. If $e_1= \hat{e}_1^{(2)}$ or $e_2= \hat{e}_2^{(2)}$, $P^-$ and $P^+$ coincide and form a saddle-node point. Separatrices of the saddle point subdivide the first quadrant into two regions, one is the basin of attraction of $O$ while the other is that of $P^+$. In other situations, all positive solutions of \eqref{eq61} converge to $O$. {\rm (ii)} If $e_1 \le \hat{e}_1^{(1)}$ or $e_2 \le \hat{e}_2^{(1)}$, then equilibrium $O$ is globally asymptotically stable. \end{theorem} \subsection*{Acknowledgements} This work was supported by NSF of Guangdong grants S2012010010320 and 1414050000636, and by STF of Guangzhou grant 1563000413. \begin{thebibliography}{99} \bibitem{Agra07} A. A. Agrawal; \emph{Filling key gaps in population and community ecology}, Front. Ecol. Environ., 5(2007): 145-152. \bibitem{Agui07} C. Aguilar, H. Vlamakis, R. Losick, R. Kolter; \emph{Thinking about bacillus subtilis as a multicellular organism}, Curr Opin Microbiol. 10(2007): 638-643. \bibitem{Alle30} W. C. Allee; \emph{Animal Aggregations}, University of Chicago Press, 1930. \bibitem{Bego06} M. Begon, C. R. Townsend, J. L. Harper; \emph{Ecology: From Individuals to Ecosystems}. Wiley, New York, 2006. \bibitem{Harr74} T. E. Harris; \emph{Contact interaction on a lattice}, Annals of Probability 2(1974): 969-988. \bibitem{Hern98} M. J. Hernandez; \emph{Dynamics of transitions between population interactions: a non-linear interaction $\alpha$-function defined}, Proc. R. Soc. London B 265(1998): 1433-1440. \bibitem{Hofb98} J. Hofbauer, K. Sigmund; \emph{Evolutionary Games and Population Dynamics}, Cambridge University Press, Cambridge, UK, 1998. \bibitem{Holl10} J. N. Holland, D. L. DeAngelis; \emph{A consumer-resource approach to the density-dependent population dynamics of mutualism}, Ecology, 91(2010): 1286-1295. \bibitem{Iwat11} S. Iwata, K. Kobayashi, S, Higa, J. Yoshimura, K. Tainaka; \emph{A simple population theory for mutualism by the use of lattice gas model}, Ecological Modelling, 222(2011): 2042-2048. \bibitem{Kell06} L. Keller, M. G. Surette; \emph{Communication in bacteria: an ecological and evolutionary perspective}. Nature Reviews Microbiology 4( 2006): 249-258. \bibitem{Lara12} T. Lara, J. Rebaza; \emph{Dynamics of transition in population interaction}, Nonlinear Analysis: Real World Applications, 13(2012): 1268-1277. \bibitem{Madi00} M. T. Madigan, J. M. Martinko, J. Parker; \emph{Biology of Microorganisms}. Prentice Hall Inc., New Jersey, 2000. \bibitem{Mats92} H. Matsuda, N. Ogita, A. Sasaki, K. Sato; \emph{Statistical mechanics of population: the lattice Lotka-Volterra model}, Progress of Theoretical Physics 88(1992): 1035-1049. \bibitem{Murray} J. D. Murray; \emph{Mathematical Biology}, Springer-Verlag, New York, 1998. \bibitem{Perk01} L. Perko; \emph{Differential Equations and Dynamical Systems}, Springer-Verlag, New York, 2001. \bibitem{Seo11} G. Seo, D. L. DeAngelis; \emph{A predator-prey model with a Holling type I functional response including a predator mutual interference}, J. Nonlinear Science. 21(2011): 811-833. \bibitem{Shap95} J. A. Shapiro; \emph{The significance of bacterial coloby patterns}, BioEssays. 17(1995): 597-607. \bibitem{Smith95} H. L. Smith; \emph{Monotone Dynamical Systems}, Amer Mathematical Society. 1995. \bibitem{Smith11} H. L. Smith, H. R.Thieme; \emph{Dynamical Systems and Population Persistence}, Amer Mathematical Society, 2011. \bibitem{Tain04} K. Tainaka, M. Kushida, Y. Ito, J. Yoshimura; \emph{Phase interspecific segregation in a lattice ecosystem with intraspecific competition}, Journal of the Physical Society of Japan 73(2004): 2914-2915. \bibitem{Wang13} Y. Wang, H. Wu; \emph{Invasibility of nectarless flowers in plant-pollinator systems}, Bull. Math. Biol. 75(2013): 1138-1156. \bibitem{Wang15} Y. Wang, D. L. DeAngelis, J. N. Holland; \emph{Dynamics of an ant-plant-pollinator model}, Commun. Non. Sci. Nume. Simu. 20 (2015): 950-964. \bibitem{Zhang03} Z. Zhang; \emph{Mutualism or cooperation among competitors promotes coexistence and competitive ability}, Ecological Modelling. 164(2003): 271-282. \end{thebibliography} \end{document} \textbf{\LARGE Figure Captions} \\\\ \textbf{Captions of Figure 1.} Intersection points of parabolas $L_1$ and $L_2$ on the whole plane. Let $m_1=m_2=0.02, e_1=e_2=40, $ $c_1=c_2=0.001, d_1=0.8, d_2=0.75$. Then $L_1$ and $L_2$ have four intersection points, while two of them are in the second and fourth quadrants, respectively. \\\\ \textbf{Captions of Figure 2.} Population dynamics of \eqref{eq21}. Red and green curves are isoclines $L_1$ and $L_2$, and black lines are separatrices of saddle points. Grey arrows denote the direction and strength of the vector fields on the portrait, while filled and open circles represent the stable and unstable equilibria, respectively. (a) Let $m_1= m_2=1, e_1=e_2=10, $ $c_1=c_2=0.001, d_1=d_2=0.8$. Then system \eqref{eq21} has a unique positive equilibrium, which is globally asymptotically stable. (b-c) Fix $m_1= m_2=1, $ $c_1=c_2=0.001, d_1=d_2=0.8$ and let $e_1, e_2$ vary. When $e_1=e_2=10$ in (b), there are two positive equilibria $P^-(0.093,0.093)$ and $P^+(0.266,0.266)$. $P^-$ is unstable while $P^+$ is asymptotically stable. When $e_1=e_2=8.7$ in (c), $P^-$ and $P^+$ coincide and form a saddle-node point $P^+( 0.17,0.17)$. Separatrices of the saddle points subdivide the first quadrant into two regions, which are basins of attraction of $O$ and $P^+$ respectively. When $e_1=e_2=8$ in (d), all solutions converge to $O$. \\\\ \textbf{Captions of Figure 3.} (a) Population dynamics of \eqref{eq21}. Let $m_1=0.5, m_2=1.5$, $c_i=0.001, d_i=0.8, e_i=8, i=1,2$. Then system \eqref{eq21} has a unique positive equilibrium, which is globally asymptotically stable. (b) Fix $m_1=0.95, m_2=1.15, e_2=5, $ $c_1=c_2=0.001, d_1=d_2=0.8$ and Let $e_1$ vary. When $e_1$ increases from 2, 3.5 to 6.5, isocline $L_1$ increases monotonically while species 2 could persist at $e_1\ge 3.5$. (c) Fix $m_1=0.95, m_2=1.15, e_1=4.5, e_2=5, $ $c_1=0.001, d_1=0.8, d_2=0.02$ and let $c_2$ vary. When $c_2$ decreases from 5.5, 2.5 to 0.02, isocline $L_2$ increases monotonically while species 2 could persist at $c_2 =2.5, 0.02$. (d) Fix $m_1=0.95, m_2=1.15, e_1=4.5, e_2=5, $ $c_1=c_2=0.001, d_1=0.8$ and let $d_2$ vary. When $d_2$ decreases from 5.8, 0.9 to 0.02, isocline $L_2$ increases monotonically while species 2 could persist at $d_2 = 0.9, 0.02$. \\\\ \textbf{Captions of Figure 4.} Fix $m_1=0.2, m_2=0.8, e_1=5, $ $c_1=c_2=0.001, d_1=d_2=0.8$ and let $e_2$ vary. When $e_2=0.8$, species 2 goes to extinction. When $e_2=5$, the two species coexist and form a win-win situation. When $e_2 =500 ($ i.e. $ e_2 \to +\infty)$, species 1 goes to extinction. %------------------------figure 2--------------------------------------- \begin{figure}[htb] \begin{center} % \includegraphics[width = .9\textwidth]{fig2} \end{center} \caption{ Population dynamics of \eqref{eq21}. Red and green curves are isoclines $L_1$ and $L_2$, and black lines are separatrices of saddle points. Grey arrows denote the direction and strength of the vector fields on the portrait, while filled and open circles represent the stable and unstable equilibria, respectively. (a) Let $m_1= m_2=1, e_1=e_2=10, $ $c_1=c_2=0.001, d_1=d_2=0.8$. Then system \eqref{eq21} has a unique positive equilibrium, which is globally asymptotically stable. (b-c) Fix $m_1= m_2=1.2, $ $c_1=c_2=0.001, d_1=d_2=0.8$ and let $e_1, e_2$ vary. When $e_1=e_2=10$ in (b), there are two positive equilibria $P^-$ and $P^+$. $P^-$ is unstable while $P^+$ is asymptotically stable. When $e_1=e_2=4.75$ in (c), $P^-$ and $P^+$ coincide and form a saddle-node point $P^+$. Separatrices of the saddle points subdivide the first quadrant into two regions, which are basins of attraction of $O$ and $P^+$ respectively. When $e_1=e_2=4.5$ in (d), all solutions converge to $O$.} \label{fig2} \end{figure} %------------------------figure 2--------------------------------------- %------------------------figure 3--------------------------------------- \begin{figure}[htb] \begin{center} % \includegraphics[width = 0.9\textwidth]{fig3} \end{center} \caption{ (a) Population dynamics of \eqref{eq21}. Let $m_1=0.5, m_2=1.5$, $c_i=0.001, d_i=0.8, e_i=8, i=1,2$. Then system \eqref{eq21} has a unique positive equilibrium, which is globally asymptotically stable. (b) Fix $m_1=0.95, m_2=1.15, e_2=5, $ $c_1=c_2=0.001, d_1=d_2=0.8$ and Let $e_1$ vary. When $e_1$ increases from 2, 3.5 to 6.5, isocline $L_1$ increases monotonically while species 2 could persist at $e_1 \ge 3.5$. (c) Fix $m_1=0.95, m_2=1.15, e_1=4.5, e_2=5, $ $c_1=0.001, d_1=0.8, d_2=0.02$ and let $c_2$ vary. When $c_2$ decreases from 5.5, 2.5 to 0.02, isocline $L_2$ increases monotonically while species 2 could persist at $c_2 =2.5, 0.02$. (d) Fix $m_1=0.95, m_2=1.15, e_1=4.5, e_2=5, $ $c_1=c_2=0.001, d_1=0.8$ and let $d_2$ vary. When $d_2$ decreases from 5.8, 0.9 to 0.02, isocline $L_2$ increases monotonically while species 2 could persist at $d_2 = 0.9, 0.02$. } \label{fig3} \end{figure} %------------------------figure 3--------------------------------------- %------------------------figure 4--------------------------------------- \begin{figure}[htb] \begin{center} % \includegraphics[width = 0.6\textwidth]{fig4} \end{center} \caption{ Fix $m_1=0.2, m_2=0.8, e_1=5, $ $c_1=c_2=0.001, d_1=d_2=0.8$ and let $e_2$ vary. When $e_2=0.8$, species 2 goes to extinction. When $e_2=5$, the two species coexist and form a win-win situation. When $e_2 =500 ($ i.e. $ e_2 \to +\infty)$, species 1 goes to extinction.} \label{fig4} \end{figure} %------------------------figure 4--------------------------------------- \end{document}