Mortality can produce limit cycles in density-dependent models with a predator-prey relationship

Tahani Mtar; Radhouane Fekih-Salem; Tewfik Sari

doi:10.3934/dcdsb.2022049

Article Contents

2022, Volume 27, Issue 12: 7445-7467. Doi: 10.3934/dcdsb.2022049

This issue Previous Article Mathematical analysis of a triple age dependent epidemiological model with including a protection strategy Next Article Factorization method for inverse time-harmonic elastic scattering with a single plane wave

Mortality can produce limit cycles in density-dependent models with a predator-prey relationship

a.
University of Tunis El Manar, National Engineering School of Tunis, LAMSIN, 1002, Tunis, Tunisia
b.
ITAP, Univ Montpellier, INRAE, Institut Agro, Montpellier, France
c.
University of Monastir, Higher Institute of Computer Science of Mahdia, 5111, Mahdia, Tunisia

^* Corresponding author: Radhouane Fekih-Salem
^* Corresponding author: Radhouane Fekih-Salem

Received: June 30, 2021

Revised: December 31, 2021

Early access: March 2022

Published: December 2022

Supported by the UNESCO ICIREWARD project ANUMAB and the Euro-Mediterranean research network TREASURE (http://www.inrae.fr/treasure).

Abstract Full Text(HTML) Figure(13) / Table(8) Related Papers Cited by

Abstract

We study an interspecific, density-dependent model of two species competing for a single nutrient in a chemostat, allowing for a predator-prey relationship between them. We have previously examined the system in the absence of species mortality, showing that multiple positive steady states can appear and disappear through a saddle-node or transcritical bifurcation. In this paper we include mortality. We give a complete analysis for the existence and local stability of all steady states of the three-dimensional system that cannot be reduced to two dimensional ones. Specializing the forms of the rate functions, we show how mortality destabilizes the positive steady state and that stable limit cycles emerge through supercritical Hopf bifurcations. To describe how the process behaves with respect to the choice of dilution rate and input concentration as control parameters, we determine the operating diagram theoretically and also numerically by using the software package MATCONT. The bifurcation diagram based on the input concentration shows various types of bifurcations of steady states, and coexistence either at a positive steady state or via sustained oscillations.

Keywords:

Mathematics Subject Classification: Primary: 34A34, 34D20; Secondary: 92B05, 92D25.

Citation:

Full Text(HTML)

Figure 1. Existence of solutions of equation $ \widetilde{f}_2(x_1,0) = 0 $. The color in this figure and all other Figures is available online only

Download: Full-size image PowerPoint slide

Figure 2. Number of positive steady states: (a) $ \mbox{Case 1} $: no positive steady state when $ (S_{in}, D) = (0.26,0.1) $, (b) $ \mbox{Case 2} $: an odd number when $ (S_{in}, D) = (0.35,0.1) $, $ {} (c) $ $ \mbox{Case 3} $: an even number when $ (S_{in}, D) = (4.5,0.8) $

Download: Full-size image PowerPoint slide

Figure 3. Curves of the function $ C_3(S) $ for different values of $ D $ when $ S>\max(\lambda_1(D),\bar{\lambda}_2(D)) $

Download: Full-size image PowerPoint slide

Figure 4. (a) Steady-state characteristics describing the local asymptotic behavior of the positive steady state $ \mathcal{E}^* $ when $ D = 0.25 $. Magnification for (b) $ S \leqslant 2 $ and (c) $ 0.32 \leqslant S \leqslant 0.36 $

Download: Full-size image PowerPoint slide

Figure 5. Case $ D = 0.25 \leqslant D^* \simeq 0.2648 $: (a) Change of sign of $ C_4 $ when $ S_1\simeq 0.3299 $ (or equivalently $ \sigma_3 \simeq 0.5255 $), $ S_2\simeq 0.3423 $ (or equivalently $ \sigma_4 \simeq 0.7159 $) and $ S_3\simeq 1.4365 $ (or equivalently $ \sigma_5 \simeq 12.4809 $). (b) Magnification for $ S \in [0.316,1.7] $ and (c) magnification for $ S \in [0.316,0.361] $

Download: Full-size image PowerPoint slide

Figure 6. Operating diagram of (1) in MAPLE. (b)-(e) Magnification when $ (S_{in},D)\in [0,2.6]\times[0,0.3] $. (c)-(f) Magnification when $ (S_{in},D)\in [0,0.6]\times[0,0.3] $

Download: Full-size image PowerPoint slide

Figure 7. Operating diagram of (1) in MATCONT. (b) Magnification when $ (S_{in},D)\in [0,2.6]\times[0,0.3] $. (c) Magnification when $ (S_{in},D)\in [0,0.6]\times[0,0.3] $

Download: Full-size image PowerPoint slide

Figure 8. Operating diagram of (1) showing the disappearance of the region $ \mathcal{J}_4 $ when $ a_1 $ and $ a_2 $ diminish

Download: Full-size image PowerPoint slide

Figure 9. Scilab simulation showing projections of the $ \omega $-limit set in variable $ S $ when $ D = 0.25 $: (a) emergence and the disappearance of limit cycle at $ \sigma_3 $ and $ \sigma_4 $ for $ S_{in}\in[0.3,0.8] $; (b) emergence of limit cycle at $ \sigma_5 $ for $ S_{in}\in[0.8,30] $

Download: Full-size image PowerPoint slide

Figure 10. (a) Variation of the pair of complex-conjugate eigenvalues (33) as $ S_{in} $ increases from 0 to 40 when $ D = 0.25 $. (b) Magnification on $ \lambda_{\pm} $ for $ S_{in}\in [0.4,0.8] $

Download: Full-size image PowerPoint slide

Figure 11. Case $ S_{in} = 0.5 \in (\sigma_2,\sigma_3) $ and $ D = 0.25 $: convergence to $ \mathcal{E}^* $

Download: Full-size image PowerPoint slide

Figure 12. Case $ S_{in} = 0.6 \in (\sigma_3,\sigma_4) $ and $ D = 0.25 $: convergence towards a stable limit cycle (in red)

Download: Full-size image PowerPoint slide

Figure 13. Case $ S_{in} = 15 >\sigma_5 $ and $ D = 0.25 $: convergence to a stable limit cycle showing the sustained oscillations

Download: Full-size image PowerPoint slide

Table 1. Necessary and sufficient conditions of existence and stability of steady states of system (1) where $ c_4 $ is defined by (20)

	Existence	Local stability
$ \mathcal{E}_0 $	always exists	$ S_{in}<\lambda_1(D) $
$ \mathcal{E}_1 $	$ S_{in}>\lambda_1(D) $	$ S_{in}< \varphi(D) $
$ \mathcal{E}^* $	(16) has a solution	$ F_2^{'}\left(x_1^\right)>F_1^{'}\left(x_1^\right) $ and $ c_4(S_{in},D)>0 $

| Show Table

DownLoad: CSV

Table 2. The set $ \Upsilon $ and the corresponding colors in Figs. 6 and 7 where $ \varphi(D) $ and $ c_4 $ are defined by (27) and (20), resp

$ \Upsilon $	Color
$ \Upsilon_1=\left\{ (S_{in},D): S_{in} = \lambda_1(D) \right\} $	Black
$ \Upsilon_2=\left\{ (S_{in},D): S_{in} = \varphi(D) \right\} $	Blue
$ \Upsilon_3:=\left\{(S_{in},D): S_{in}=S_{in}^{\rm SN}(D) \right\} $
$ \Upsilon_4=\left\{ (S_{in}, D): c_4(S_{in}, D)=0 \right\} $	Green

| Show Table

DownLoad: CSV

Table 3. Existence and stability of steady states according to regions in the operating diagrams of Figs. 6 and 7. The letter S (resp. U) means stable (resp. unstable) steady state. Absence of letter means that the corresponding steady state does not exist

Condition	Region	Color	$ \mathcal{E}_0 $	$ \mathcal{E}_1 $	$ \mathcal{E}^* $
$ S_{in}<\lambda_1(D) $	$ \mathcal{J}_1 $	Cyan	S
$ \lambda_1(D)<S_{in}<\varphi(D) $	$ \mathcal{J}_2 $	Pink	U	S
$ \varphi(D)<S_{in} $ and $ c_4(S_{in}, D)>0 $	$ \mathcal{J}_3 $	Grey	U	U	S
$ \varphi(D)<S_{in} $ and $ c_4(S_{in}, D)<0 $	$ \mathcal{J}_4 $	Yellow	U	U	U

| Show Table

DownLoad: CSV

Table 4. Nature of bifurcations of the steady states of (1) by crossing to the surfaces of $ \Upsilon $. The letter TB (resp. SHB) means a transcritical bifurcation (resp. Supercritical Hopf bifurcation)

Subset	Transition	Bifurcation
$ \Upsilon_1 $	$ \mathcal{J}_1 $ to $ \mathcal{J}_2 $	TB: $ \mathcal{E}_0=\mathcal{E}_1 $
$ \Upsilon_2 $	$ \mathcal{J}_2 $ to $ \mathcal{J}_3 $	TB: $ \mathcal{E}_1=\mathcal{E}^* $
$ \Upsilon_4 $	$ \mathcal{J}_3 $ to $ \mathcal{J}_4 $	SHB: $ \mathcal{E}^* $

| Show Table

DownLoad: CSV

Table 5. Definitions of the critical values $ \sigma_i $, $ i = 1,\ldots,5 $ of $ S_{in} $ and the corresponding nature of bifurcations when $ D = 0.25 $

Definition	Value	Bifurcation
$ \sigma_1=\lambda_1(D) $	0.31884	TB
$ \sigma_2=\varphi(D) $	0.35394	TB
$ \sigma_3 $ is the first solution of equation $ c_4(S_{in})=0 $	0.52555	SHB
$ \sigma_4 $ is the second solution of equation $ c_4(S_{in})=0 $	0.71593	SHB
$ \sigma_5 $ is the third solution of equation $ c_4(S_{in})=0 $	12.4809	SHB

| Show Table

DownLoad: CSV

Table 6. Existence and stability of steady states according to $ S_{in} $

Interval of $ S_{in} $	$ \mathcal{E}_0 $	$ \mathcal{E}_1 $	$ \mathcal{E}^* $
$ (0,\sigma_1) $	S
$ (\sigma_1,\sigma_2) $	U	S
$ (\sigma_2,\sigma_3) $	U	U	S
$ (\sigma_3,\sigma_4) $	U	U	U
$ (\sigma_4,\sigma_5) $	U	U	S
$ (\sigma_5,+\infty) $	U	U	U

| Show Table

DownLoad: CSV

Table 7. Parameter values used for model (1) when the growth rates $ f_1 $ and $ f_2 $ are given by (31).

Parameter	$ m_1 $	$ K_1 $	$ L_1 $	$ m_2 $	$ K_2 $	$ L_2 $	$ \alpha_1 $	$ \alpha_2 $	$ a_1 $	$ a_2 $
Fig. 2(c)	2.75	2	1.2	2.95	1.8	1.5	$ 10^{-3} $	0.1	0.95	0.7
Figs. 1, 2(a, b), 6, 7, 9–13									0.3	0.2
Fig. 8(a)	4	2	3	8	0.1	0.2	1	1	0.3	0.05
Fig. 8(b)	4	2	3	8	0.1	0.2	1	1	0.1	0.05
Fig. 8(c)									0	0

| Show Table

DownLoad: CSV

Table 8. Break-even concentration, notations and auxiliary functions

λ₁(D)	S = λ₁(D) is the solution of equation f₁(S, 0)=α₁D+a₁
λ₁(D)	It is defined for D < , (f₁(+∞, 0) − a₁)=α₁, see(4)
$\tilde{x}_{1} $	$\tilde{x}_{1}=\frac{D}{D_{1}}\left(S_{i n}-\lambda_{1}(D)\right), \text { see }(5)$
$\widetilde{f}_{1}\left(x_{1}, x_{2}\right)$	$\widetilde{f}_{1}\left(x_{1}, x_{2}\right)=f_{1}\left(S_{i n}-\frac{D_{1}}{D} x_{1}-\frac{D_{2}}{D} x_{2}, x_{2}\right)-D_{1}, \text { see }(11) $
$\widetilde{f}_{2}\left(x_{1}, x_{2}\right)$	$\widetilde{f}_{2}\left(x_{1}, x_{2}\right)=f_{2}\left(S_{i n}-\frac{D_{1}}{D} x_{1}-\frac{D_{2}}{D} x_{2}, x_{1}\right)-D_{2}, \text { see }(11)$
F₁(x₁)	x₂ = F₁(x₁) is the unique solution of equation $\widetilde{f}_{1}$(x₁, x₂) = 0
F₁(x₁)	It is defined for 0 ≤ x₁ ≤ $\tilde{x}_{1} $, see Lemma 1
F₂(x₁)	x₂ = F₂(x₁) is the unique solution of equation $\widetilde{f}_{2}$ (x₁, x₂) = 0
	It is defined for $x_{1}^{1}$ ≤ x₁ ≤ $x_{1}^{2}$, where $x_{1}^{1}$ and $x_{1}^{2}$ are
	the solutions of equation $\widetilde{f}_{2}$(x₁, 0) = 0, see Lemma 3
$\left(x_{1}^{}, x_{2}^{}\right)$	$\left(x_{1}^{}, x_{2}^{}\right)$ is a solution of x₂ = F₁(x₁) = F₂(x₁), see Prop. 3
X₁(S, D)	x₁ = X₁(S, D) is the solution of equation f₂(S, x₁) = α₂D + a₂
	It is defined for S > λ₂(D), where λ₂(D) is the unique solution,
	if it exists, of equation f₂(S, +∞) = α₂D + a₂, see (26)
X₂(S, D)	x₂ = X₂(S, D) is the solution of equation f₁(S, x₂) = α₁D + a₁
	It is defined for λ₁(D) ≤ S < λ₁(D), where λ₁(D) is the unique solution,
	if it exists, of equation f₁(S, +∞) = α₁D + a₁, see (25)
$\varphi$(D)	$\varphi(D)=\lambda_{1}(D)+\frac{D_{1}}{D} X_{1}\left(\lambda_{1}(D), D\right), \text { see }(27)$
K(S, D)	K(S, D) = D₁X₁(S, D) + D₂X₂(S, D), see (30)

| Show Table

DownLoad: CSV

Related Papers

Cited by

References

[1]	N. Abdellatif, R. Fekih-Salem and T. Sari, Competition for a single resource and coexistence of several species in the chemostat, Math. Biosci. Eng., 13 (2016), 631-652. doi: 10.3934/mbe.2016012.
[2]	M. Ballyk, R. Staffeldt and I. Jawarneh, A nutrient-prey-predator model: Stability and bifurcations, Discrete Contin. Dyn. Syst. Ser. S, 13 (2020), 2975-3004. doi: 10.3934/dcdss.2020192.
[3]	B. Bar and T. Sari, The operating diagram for a model of competition in a chemostat with an external lethal inhibitor, Discrete Contin. Dyn. Syst. Ser. B, 25 (2020), 2093-2120. doi: 10.3934/dcdsb.2019203.
[4]	B. Benyahia, T. Sari, B. Cherki and J. Harmand, Bifurcation and stability analysis of a two step model for monitoring anaerobic digestion processes, J. Proc. Control, 22 (2012), 1008-1019. doi: 10.1016/j.jprocont.2012.04.012.
[5]	O. Bernard, Z. Hadj-Sadok, D. Dochain, A. Genovesi and J.-P. Steyer, Dynamical model development and parameter identification for an anaerobic wastewater treatment process, Biotechnol. Bioeng., 75 (2001), 424-438. doi: 10.1002/bit.10036.
[6]	M. P. Boer, B. W. Kooi and S. A. L. M. Kooijman, Food chain dynamics in the chemostat, Math. Biosci., 150 (1998), 43-62. doi: 10.1016/S0025-5564(98)00010-8.
[7]	F. Borsali and K. Yadi, Contribution to the study of the effect of the interspecificity on a two nutrients competition model, Int. J. Biomath., 8 (2015), 1550008, 17 pp. doi: 10.1142/S1793524515500084.
[8]	M. Dali-Youcef, A. Rapaport and T. Sari, Study of performance criteria of serial configuration of two chemostats, Math. Biosci. Eng., 17 (2020), 6278-6309. doi: 10.3934/mbe.2020332.
[9]	Y. Daoud, N. Abdellatif, T. Sari and J. Harmand, Steady state analysis of a syntrophic model: The effect of a new input substrate concentration, Math. Model. Nat. Phenom., 13 (2018), Paper No. 31, 21 pp. doi: 10.1051/mmnp/2018037.
[10]	P. De Leenheer, D. Angeli and E. D. Sontag, Crowding effects promote coexistence in the chemostat, J. Math. Anal. Appl., 319 (2006), 48-60. doi: 10.1016/j.jmaa.2006.02.036.
[11]	M. Dellal and B. Bar, Global analysis of a model of competition in the chemostat with internal inhibitor, Discrete Contin. Dyn. Syst. Ser. B, 26 (2021), 1129-1148. doi: 10.3934/dcdsb.2020156.
[12]	M. Dellal, M. Lakrib and T. Sari, The operating diagram of a model of two competitors in a chemostat with an external inhibitor, Math. Biosci., 302 (2018), 27-45. doi: 10.1016/j.mbs.2018.05.004.
[13]	A. Dhooge, W. Govaerts, Y. A. Kuznetsov, H. G. E. Meijer and B. Sautois, New features of the software MatCont for bifurcation analysis of dynamical systems, Math. Comput. Model. Dyn. Syst., 14 (2008), 147-175. doi: 10.1080/13873950701742754.
[14]	M. El-Hajji, How can inter-specific interferences explain coexistence or confirm the competitive exclusion principle in a chemostat?, Int. J. Biomath., 11 (2018), 1850111, 20 pp. doi: 10.1142/S1793524518501115.
[15]	M. El-Hajji, F. Mazenc and J. Harmand, A mathematical study of a syntrophic relationship of a model of anaerobic digestion process, Math. Biosci. Eng., 7 (2010), 641-656. doi: 10.3934/mbe.2010.7.641.
[16]	R. Fekih-Salem, J. Harmand, C. Lobry, A. Rapaport and T. Sari, Extensions of the chemostat model with flocculation, J. Math. Anal. Appl., 397 (2013), 292-306. doi: 10.1016/j.jmaa.2012.07.055.
[17]	R. Fekih-Salem, C. Lobry and T. Sari, A density-dependent model of competition for one resource in the chemostat, Math. Biosc., 268 (2017), 104-122. doi: 10.1016/j.mbs.2017.02.007.
[18]	R. Fekih-Salem, A. Rapaport and T. Sari, Emergence of coexistence and limit cycles in the chemostat model with flocculation for a general class of functional responses, Appl. Math. Modell., 40 (2016), 7656-7677. doi: 10.1016/j.apm.2016.03.028.
[19]	R. Fekih-Salem and T. Sari, Properties of the chemostat model with aggregated biomass and distinct removal rates, SIAM J. Appl. Dyn. Syst. (SIADS), 18 (2019), 481-509. doi: 10.1137/18M1171801.
[20]	R. Fekih-Salem and T. Sari, Operating diagram of a flocculation model in the chemostat, ARIMA J., 31 (2020), 45-58. doi: 10.46298/arima.5593.
[21]	B. Haegeman and A. Rapaport, How flocculation can explain coexistence in the chemostat, J. Biol. Dyn., 2 (2008), 1-13. doi: 10.1080/17513750801942537.
[22]	M. Hanaki, J. Harmand, Z. Mghazli, A. Rapaport, T. Sari and P. Ugalde, Mathematical study of a two-stage anaerobic model when the hydrolysis is the limiting step, Processes, 9 (2021). doi: 10.3390/pr9112050.
[23]	S. R. Hansen and S. P. Hubbell, Single-nutrient microbial competition: Qualitative agreement between experimental and theoretically forecast outcomes, Science, 207 (1980), 1491-1493. doi: 10.1126/science.6767274.
[24]	J. Harmand, C. Lobry, A. Rapaport and T. Sari, The Chemostat: Mathematical Theory of Microorganism Cultures, Chemostat and bioprocesses set., Vol. 1. ISTE, London; John Wiley & Sons, Inc., Hoboken, NJ, 2017. doi: 10.1002/9781119437215.
[25]	J. Harmand, A. Rapaport, D. Dochain and C. Lobry, Microbial ecology and bioprocess control: Opportunities and challenges, Journal of Process Control, 18 (2008), 865-875. doi: 10.1016/j.jprocont.2008.06.017.
[26]	S.-B. Hsu, C. A. Klausmeier and C.-J. Lin, Analysis of a model of two parallel food chains, Discrete Contin. Dyn. Syst. Ser. B, 12 (2009), 337-359. doi: 10.3934/dcdsb.2009.12.337.
[27]	Z. Khedim, B. Benyahia, B. Cherki, T. Sari and J. Harmand, Effect of control parameters on biogas production during the anaerobic digestion of protein-rich substrates, Appl. Math. Model., 61 (2018), 351-376. doi: 10.1016/j.apm.2018.04.020.
[28]	B. W. Kooi and M. P. Boer, Chaotic behaviour of a predator-prey system in the chemostat, Dyn. Contin. Discrete Impuls. Syst. Ser. B Appl. Algorithms, 10 (2003), 259-272.
[29]	B. Li and Y. Kuang, Simple food chain in a chemostat with distinct removal rates, J. Math. Anal. Appl., 242 (2000), 75-92. doi: 10.1006/jmaa.1999.6655.
[30]	C. Lobry and J. Harmand, A new hypothesis to explain the coexistence of n species in the presence of a single resource, C. R. Biol., 329 (2006), 40-46. doi: 10.1016/j.crvi.2005.10.004.
[31]	C. Lobry and F. Mazenc, Effect on persistence of intra-specific competition in competition models, Electron. J. Diff. Equ., 125 (2007), 1-10.
[32]	C. Lobry, F. Mazenc and A. Rapaport, Persistence in ecological models of competition for a single resource, C. R. Acad. Sci. Paris, Ser. I, 340 (2005), 199-204. doi: 10.1016/j.crma.2004.12.021.
[33]	C. Lobry, A. Rapaport and F. Mazenc, Sur un modèle densité-dépendant de compétition pour une ressource, C. R. Biol., 329 (2006), 63-70. doi: 10.1016/j.crvi.2005.11.004.
[34]	MAPLE [Software], Version 13.0, Maplesoft, a division of Waterloo Maple Inc., Waterloo, Ontario, (2009).
[35]	MATCONT [Software], Matcontsoft, (2021).
[36]	T. Mtar, R. Fekih-Salem and T. Sari, Interspecific density-dependent model of predator-prey relationship in the chemostat, Int. J. Biomath., 14 (2021), 2050086, 22 pp. doi: 10.1142/S1793524520500862.
[37]	T. Mtar, R. Fekih-Salem and T. Sari, Effect of the mortality on a density-dependent model with a predator-prey relationship, CARI'2020, Proceedings of the 15th African Conference on Research in Computer Science and Applied Mathematics, (2020).
[38]	S. Nouaoura, N. Abdellatif, R. Fekih-Salem and T. Sari, Mathematical analysis of a three-tiered model of anaerobic digestion, SIAM J. Appl. Math., 81 (2021), 1264-1286. doi: 10.1137/20M1353897.
[39]	S. Nouaoura, R. Fekih-Salem, N. Abdellatif and T. Sari, Mathematical analysis of a three-tiered food-web in the chemostat, Discrete Contin. Dyn. Syst. Ser. B, 26 (2021), 5601-5625. doi: 10.3934/dcdsb.2020369.
[40]	S. Nouaoura, R. Fekih-Salem, N. Abdellatif and T. Sari, Operating diagrams for a three-tiered microbial food web in the chemostat, Preprint HAL, (2021).
[41]	A. Rapaport and M. Veruete, A new proof of the competitive exclusion principle in the chemostat, Discrete Contin. Dyn. Syst. Ser. B, 24 (2019), 3755-3764. doi: 10.3934/dcdsb.2018314.
[42]	T. Sari and B. Benyahia, The operating diagram for a two-step anaerobic digestion model, Nonlinear Dynam., 105 (2021), 2711-2737. doi: 10.1007/s11071-021-06722-7.
[43]	T. Sari and J. Harmand, A model of a syntrophic relationship between two microbial species in a chemostat including maintenance, Math. Biosci., 275 (2016), 1-9. doi: 10.1016/j.mbs.2016.02.008.
[44]	T. Sari and M. J. Wade, Generalised approach to modelling a three-tiered microbial food-web, Math. Biosci., 291 (2017), 21-37. doi: 10.1016/j.mbs.2017.07.005.
[45]	M. Sbarciog, M. Loccufier and E. Noldus, Determination of appropriate operating strategies for anaerobic digestion systems, Biochem. Eng. J., 51 (2010), 180-188. doi: 10.1016/j.bej.2010.06.016.
[46]	SCILAB [Software], version 6.0.1, Scilab, Enterprises SAS (2018).
[47]	S. Shen, G. C. Premier, A. Guwy and R. Dinsdale, Bifurcation and stability analysis of an anaerobic digestion model, Nonlinear Dynam., 48 (2007), 391-408. doi: 10.1007/s11071-006-9093-1.
[48]	H. L. Smith and P. Waltman, The Theory of the Chemostat, Dynamics of Microbial Competition, Cambridge University Press, 1995. doi: 10.1017/CBO9780511530043.
[49]	S. Sobieszek, M. J. Wade and G. S. K. Wolkowicz, Rich dynamics of a three-tiered anaerobic food-web in a chemostat with multiple substrate inflow, Math. Biosci. Eng., 17 (2020), 7045-7073. doi: 10.3934/mbe.2020363.
[50]	G. A. K. Van Voorn, B. W. Kooi and M. P. Boer, Ecological consequences of global bifurcations in some food chain models, Math. Biosci., 226 (2010), 120-133. doi: 10.1016/j.mbs.2010.04.005.
[51]	D. V. Vayenas and S. Pavlou, Chaotic dynamics of a food web in a chemostat, Math. Biosci., 162 (1999), 69-84. doi: 10.1016/S0025-5564(99)00044-9.
[52]	M. J. Wade, J. Oakley, S. Harbisher, N. G. Parker and J. Dolfing, MI-Sim: A MATLAB package for the numerical analysis of microbial ecological interactions, PLoS ONE, 12 (2017), 0173249. doi: 10.1371/journal.pone.0173249.
[53]	M. J. Wade, R. W. Pattinsonm, N. G. Parker and J. Dolfing, Emergent behaviour in a chlorophenol-mineralising three-tiered microbial 'food web', J. Theor. Biol., 389 (2016), 171-186. doi: 10.1016/j.jtbi.2015.10.032.
[54]	M. Weedermann, G. S. K. Wolkowicz and J. Sasara, Optimal biogas production in a model for anaerobic digestion, Nonlinear Dyn., 81 (2015), 1097-1112. doi: 10.1007/s11071-015-2051-z.
[55]	G. S. K. Wolkowicz, Successful invasion of a food web in a chemostat, Math. Biosci., 93 (1989), 249-268. doi: 10.1016/0025-5564(89)90025-4.
[56]	A. Xu, J. Dolfing, T. P. Curtis, G. Montague and E. Martin, Maintenance affects the stability of a two-tiered microbial 'food chain'?, J. Theor. Biol., 276 (2011), 35-41. doi: 10.1016/j.jtbi.2011.01.026.