Hopf bifurcation and stability analysis of a delay differential equation model for biodegradation of a class of microcystins

1.   Introduction

With global warming and increasing eutrophication ^[1], the frequency and distribution of harmful cyanobacterial blooms in lake ecosystems are increasing worldwide ^[2,3,4]. When cyanobacteria, also known as "blue-green algae", congregate in large numbers in freshwater lakes, they form harmful algal blooms. These blooms produce a variety of secondary toxic metabolites called cyanotoxins ^[5]. Harmful algal blooms can have negative impacts on various aspects, including food security, tourism, local economies, human health, drinking water, and aquatic food ^[6,7].

Among the many different types of cyanotoxins produced by cyanobacteria, the most abundant, widespread, and harmful cyanotoxins are microcystins (MCs) ^[8,9]. MCs comprise a set of cyclic heptapeptide hepatotoxins produced by freshwater species of cyanobacteria, which share the common structure cyclo (-Adda-D-Glu-M-dha-D-Ala-L-X-D-MeAsp-L-Z), where X and Z represent variable amino acids ^[10,11]. MCs can form many different isomers due to the variation of amino acid species in the peptide composition. The most common variants detected in China are MC-LR, MC-RR, and MC-YR ^[12]. MCs can inhibit the activity of phosphoproteases, causing disruption of physiological and biochemical reactions in the human body and seriously endangering health. Long-term and frequent exposure to low-concentration MCs can lead to chronic apoptosis or uncontrolled cell proliferation of hepatocytes, which in turn promotes the occurrence and development of tumors, leading to primary hepatocellular carcinoma ^[13,14]. The International Agency for Research on Cancer (IARC) classified MC-LR as a class 2B carcinogen. The World Health Organization (WHO) recommends that the concentration of MC-LR in drinking water should not exceed 1 $\mu$g/L ^[15]. Therefore, effective removal of MCs is the key to drinking water treatment.

Since MCs have harmful effects on humans and the environment, research on their degradation has important practical significance. The degradation pathways of MCs mainly comprise physical, chemical, and biological methods. Physical degradation methods mainly include coagulation and precipitation, filtration, activated carbon adsorption, and membrane treatment. However, physical degradation methods are generally inefficient, costly, and difficult to recycle ^[16]. Chemical degradation methods, such as chemical reagents, ozone oxidation, and photodegradation, can effectively reduce MCs in water bodies. However, these methods can introduce toxic byproducts, causing secondary pollution. Furthermore, the photocatalytic degradation process has complex operational requirements ^[17]. Therefore, the limitations of physical and chemical degradation methods have restricted their further application in degrading MCs. However, biodegradation has proven to be highly efficient, cost-effective, and free from secondary pollution, making it the safest and most effective method for degrading MCs ^[18].

Microorganisms are capable of reducing or losing the toxicity of MCs by modifying the structure of the Adda active group in the side chain of MCs or by breaking down the ring structure. In 1994, Jones et al. first isolated a strain of Sphingomonas sp. ACM-3962 from natural water bodies, which could degrade MC-LR ^[19]. Subsequently, Bourne et al. showed for the first time that the microcystinases MlrA, MlrB, MlrC, and transport protein (MlrD) encoded by four genes, mlrA, mlrB, mlrC, and mlrD, respectively, were involved in the degradation of MC-LR in Sphingomonas sp. ACM-3962. The cyclic MC-LR was sequentially degraded to linear MC-LR (Adda-Glu-M-dha-Ala-Leu-MaspArg-OH), tetra-compound (Adda-Glu-M-dha-Ala-OH), poly-compound, and amino acid by the catalytic action of MlrA, MlrB, and MlrC. MlrA peptidase activity has been shown to be the most efficient enzymatic process and the most specific catalyst in all known detoxification pathways of MCs ^[20]. Since then, various strains capable of degrading MCs have been isolated, such as Pseudomonas aeruginosa ^[21], Paucibacter toxinivorans ^[22], and Brevibacterium sp. ^[23]. However, the degrading bacteria are still mainly concentrated in the genus Sphingomonas ^[24,25,26]. Yan et al. ^[27] successfully extracted the Sphingopyxis sp. USTB-05, which can effectively degrade MCs. They studied the biodegradation of MCs by USTB-05 at both cellular and enzymatic levels.

On the other hand, the studies of the dynamics of the chemostat models and their variants have achieved rich results ^[28,29]. Tai et al. ^[30] proposed a time-delayed microorganism flocculation model. They studied the existence and local stability of the equilibria of the presented model and found that the model can display forward or backward bifurcation. Guo et al. ^[31] further studied the uniform persistence of the model and global stability of the equilibria. Building upon the research in ^[31], Guo et al. ^[32] considered a time-delayed microorganism flocculation model with saturated functional responses. Using the generalized Lyapunov-LaSalle theorem, they established some conditions for the global stability of the equilibria. Song et al. ^[33] proposed a dynamic model for microbial flocculant with nutrient competition and metabolic products, and analyzed the global dynamic properties of the model.

For the flocculation and microbial degradation of MCs, Yang et al. ^[34] proposed the following model:

Here, $x_{1}(t)$, $x_{2}(t)$, and $x_{3}(t)$ represent the concentrations of MCs, Sphingomonas sp., and the degrading enzymes produced by Sphingomonas sp. at time $t$. The constant $a_{10} > 0$ is the input concentration of MCs. The constant $D > 0$ is the continuous input rate of MCs, Sphingomonas sp. and the degrading enzymes into the chemostat. The constant $a_{12}\ge0$ is the consumption rate of MCs. The constant $a_{21}\ge0$ is the maximum growth rate of Sphingomonas sp. The constant $a_{20}\ge0$ is the consumption rate of Sphingomonas sp. The constant $a_{30}\ge0$ is the rate at which Sphingomonas sp. produce the degrading enzymes that can be used for the degradation of MCs. The constant $a_{13}\ge0$ is the degradation rate of MCs. The constant $a_{31}\ge0$ is the consumption rate of the degrading enzymes. The constants $d_{1}$, $d_{2}$, and $d_{3}$ represent the death rates of MCs, Sphingomonas sp., and the degrading enzymes, respectively. Yang et al. ^[34] studied the local stability of the equilibria and the global stability of the boundary equilibrium of model (1.1), as well as the uniform persistence of model (1.1). Furthermore, Song et al. ^[35] studied the global stability of the positive equilibrium of model (1.1) by constructing suitable Lyapunov functions. They also found that, under certain conditions, the parameter $a_{13}/D$ can cause Hopf bifurcation.

Model (1.1) is a set of ordinary differential equations that is valid under the assumption that the process of nutrient (MCs) storage by the organism (Sphingomonas sp.) is instantaneous. However, during the continuous cultivation of microorganisms, the growth of microorganisms and the consumption of nutrients often show a time delay. Time delay is caused by factors such as the storage of nutrients by microorganisms and their metabolic processes. Therefore, time delay becomes an important factor in accurately characterizing microbial culture processes ^{[36,37,38,39]}. In reference ^[40], Song et al. rewrote the second equation of model (1.1) into the following form:

where the constant $\tau_{1}\ge0$ represents the time that the organism (Sphingomonas sp.) stores the nutrient (MCs). The term $e^{-\delta_{1}\tau_{1}}$ represents the approximate proportion of individuals remaining in the chemostat during the conversion process. Considering that there may be time delay in the production of degrading enzymes by microorganisms ^[41,42], and under the experimental conditions of model (1.1), in order to obtain the degrading enzymes produced by Sphingomonas sp., it is necessary to separate the degrading enzymes from the chemostat by centrifugation and sonicate the Sphingomonas sp. to release the degrading enzymes. This process is not instantaneous. Therefore, we construct the following time-delayed model:

Here, the constant $\tau_{2}\ge 0$ represents the time required for Sphingomonas sp. to produce degrading enzymes that can be used for the degradation of MCs. The term $e^{-\delta_{2}\tau_{2}}$ represents the approximate proportion of the degrading enzymes produced by Sphingomonas sp. that remains in the chemostat during the conversion process. For clarity, the biological meanings of the parameters of model (1.2) are summarized in Table 1.

In model (1.2), for $\tau_{1} > 0$ and $\tau_{2} = 0$, Song et al. ^[40] studied the local and global stability of the boundary equilibrium, local stability of the positive equilibrium, and the existence of Hopf bifurcations caused by the time delay $\tau_{1}$. However, they did not obtain global stability results for the positive equilibrium of model (1.2) for $\tau_{1} > 0$ and $\tau_{2} = 0$. The global stability or attractivity of the positive equilibrium of model (1.2) deserves further study. The main purpose of this paper is to study the effect of the time delay $\tau_2$ on the dynamical properties of model (1.2) and to obtain sufficient conditions for the global attractivity of the positive equilibrium of model (1.2).

The rest of this paper is organized as follows. In Section 2, the dimensionless model for model (1.2) is first obtained. Further, the classification of equilibria and the global stability of the boundary equilibrium are obtained for model (2.1). In Section 3, if $\tau_{1} = 0, \tau_{2} > 0$, or $\tau_{1} = \tau_{2} > 0$, we study the local stability of the positive equilibrium and the existence of Hopf bifurcations of model (2.1). In Section 4, by constructing suitable Lyapunov functionals and applying some inequality analysis techniques, we establish some sufficient conditions for the global attractivity of the positive equilibrium. In Section 5, we provide several specific examples and perform numerical simulations to validate our conclusions. Finally, using experimental data from reference ^[43], we estimate the specific values of model parameters through the least squares method and perform numerical simulations.

2.   Preliminaties

First, we make the following dimensionless transformations of model (1.2):

For convenience, we remove the superscript and obtain the following time-delayed model:

Let $\mathbb C = \mathbb C([-\hat{\tau}, 0], \mathbb R^{3})$ be the Banach space of continuous functions mapping $[-\hat{\tau}, 0]$ to $\mathbb R^{3}$, equipped with the sup-norm, where $\hat{\tau} = \max \left\{ \tau_{1}, \tau_{2} \right\}$. We assume that model (2.1) always satisfies the initial condition

where

Then, we have the following lemma.

Lemma 2.1. The solution $(x(t), y(t), z(t))^{T}$ of model (2.1) with initial condition (2.2) is existent, unique, and nonnegative on $\left[0, +\infty\right)$, which satisfies

where $l = \min \left\{ D_{1}, D_{2} \right\}$.

Clearly, model (2.1) always has a boundary equilibrium $E_{0} = (x_{0}, 0, 0)$. For convenience, let us define

When $R_{0} > 1$ ($0\le \tau_{1} < \tau_{1max}$), model (2.1) has a unique positive equilibrium $E^{*} = (x^{*}, y^{*}, z^{*})$, where

By a similar argument as in the proof of Theorem 2.1 in ^[40], we obtain the following result.

Theorem 2.1. (i) If $R_{0} < 1$, then the boundary equilibrium $E_{0}$ is globally asymptotically stable.

(ii) If $R_{0} = 1$, then the boundary equilibrium $E_{0}$ is linearly stable.

(iii)If $R_{0} > 1$, then the boundary equilibrium $E_{0}$ is unstable.

3.   Local stability of the positive equilibrium and Hopf bifurcation analysis

In this section, we assume that $R_{0} > 1$. The method provided in references ^[44,45] for studying the stability of the models with coefficient-dependent time delays is employed here. We utilize the stability criteria discussed in Section 2.4 of reference ^[45] to analyze the local stability of the positive equilibrium.

The characteristic equation of model (2.1) at the positive equilibrium $E^{*}$ can be written as

where

When $\tau _{1} = \tau _{2} = 0$, it is easy to obtain that $A_{2}(0, 0)+B_{2} > 0$, $A_{1}(0, 0)+B_{1}(0, 0) > 0$, and $A_{0}(0, 0)+B_{0}(0, 0)+C_{0}(0, 0) > 0$. Define the following condition:

According to the Routh-Hurwitz criterion, if condition $(H_{1})$ holds, then all the roots of Eq (3.1) have negative real parts and the positive equilibrium $E^{*}$ is locally asymptotically stable if $\tau _{1} = \tau _{2} = 0$.

Define

Remark 3.1. From reference ^[35], we state the following two facts.

(i) If $\Pi \ge 0$, then condition $(H_{1})$ holds.

(ii) If $\Pi < 0$, then there exists a positive constant $a_{2}^{*}$ such that if $a_{2} < a_{2}^{*}$, then condition $(H_{1})$ holds, and if $a_{2}\ge a_{2}^{*}$, then condition $(H_{1})$ does not hold.

Case Ⅰ. $\tau_{1} = 0$, $\tau_{2}\ge 0$.

For the convenience of analysis, we first assume $\tau_{1} = 0$ and choose $\tau_{2}$ as the bifucation parameter.

In this case, the characteristic equation of model (2.1) at the positive equilibrium $E^{*}$ can be written as

where

If $\tau_{2} = 0$ and condition $(H_{1})$ holds, then all the roots of Eq (3.2) have negative real parts, and the positive equilibrium $E^{*}$ is locally asymptotically stable.

When $\tau_{2} > 0$, stability switching may occur when the existence of a pair of pure imaginary roots $\lambda = \pm i \omega$ exist in Eq (3.2) that crosses the imaginary axis as the value of $\tau_{2}$ increases. It can be easily proven that $P(\lambda, \tau _{2})$ and $Q(\lambda, \tau _{2})$ in Eq (3.2) satisfy the following properties:

(ⅰ) $P(0, \tau _{2}) +Q(0, \tau _{2})\ne 0$;

(ⅱ) $P(i\omega, \tau _{2}) +Q(i\omega, \tau _{2})\ne 0$;

(ⅲ) $\lim \limits _{ \left| \lambda \right| \rightarrow \infty} \left| \frac { Q(\lambda, \tau _{2}) } {P(\lambda, \tau _{2})} \right| = 0$;

(ⅳ) $F(\omega, \tau_{2}) = |P(\omega, \tau_{2})|^{2}-|Q(\omega, \tau_{2}|^{2}$ has a finite number of zero roots;

(ⅴ) The equation $F(\omega, \tau_{2}) = 0$ has positive roots $\omega(\tau_{2})$, and each positive root is continuously differentiable.

Let $\lambda = i\omega$ be a purely imaginary root of Eq (3.2). Then we have

Further, we have

By squaring both sides of the above two equations and adding them together, we have

where

Letting $v = {\omega}^{2}$, we have

Letting $\Delta(\tau_{2}) = P_{2}^{2}(\tau_{2})-3P_{1}(\tau_{2})$. If $\Delta(\tau_{2})\geq 0$, we define $v_{+} = \frac{-P_{2}(\tau_{2})+\sqrt{\Delta(\tau_{2})}}{3}$. According to reference ^[46], the following conclusion holds.

Lemma 3.1. (i) When $P_{0}(\tau_{2}) < 0$, Eq (3.4) has at least one positive real root;

(i) When $P_{0}(\tau_{2})\ge 0$, Eq (3.4) has positive roots if and only if $v_{+} = \frac{-P_{2}(\tau_{2})+\sqrt{\Delta(\tau_{2})}}{3} > 0$ and $h(v_{+})\le 0$;

(i) When conditions (i) or (ii) are not satisfied, Eq (3.4) has no positive real roots.

For convenience, in our subsequent discussions, we assume that Eq (3.4) has only one positive root. There exists a set $I$ such that, when $\tau_{2}\in I$, the equation

has a positive real root $\omega = \omega(\tau_{2}) > 0$. Furthermore, for $\tau_{2}\in I$, we can define an angle $\theta(\tau_{2})\in \left[ 0, 2\pi\right)$ as the solution of (3.3):

For $\tau_{2}\in I$, the angle $\theta$ in the above equation must satisfy the following relationship with $\omega \tau_{2}$ in Eq (3.3):

Define

where $S_{n}(\tau_{2})$ is continuously differentiable with respect to $\tau_{2}$.

Theorem 3.1. Suppose that $\tau _{1} = 0$, $R_{0} > 1$, and condition $(H_{1})$ hold. When $\tau_{2} = \tau_{2}^{*}\in I$ such that $S_{n}(\tau_{2}^{*}) = 0$ for some $n\in \mathbb{N}_{0}$ hold, the characteristic equation (3.2) has a pair of purely imaginary roots $\pm i\omega(\tau_{2}^{*})$. If $\delta(\tau_{2}^{*}) > 0$ ($< 0$), then the roots $\pm i\omega(\tau_{2}^{*})$ cross the imaginary axis from left to right (from right to left), where

Theorem 3.2. Suppose that $\tau_{1} = 0$ and $R_{0} > 1$. If condition $(H_{1})$ holds, the following conclusions hold.

(i) When the function $S_{0}(\tau_{2})$ has no positive roots for $\tau_{2}\in I$, the positive equilibrium $E^{*}$ is locally asymptotically stable for any $\tau_{2}\ge 0$;

(ii) When the function $S_{n}(\tau_{2})$ has roots such that $\tau_{2}^{0} < \tau_{2}^{1} < \cdots < \tau_{2}^{m}\in I$ and $S^{'}_{n}(\tau_{2}^{j})\ne 0 (j = 0, 1, 2, \cdots, m)$, then the positive equilibrium $E^{*}$ is locally asymptotically stable on $\left[0, \tau_{2}^{0}\right)\cup \left(\tau_{2}^{m}, +\infty\right)\cap I$.

Case Ⅱ. $\tau_{1} = \tau_{2} = \tau$.

We assume $\tau_{1} = \tau_{2} = \tau$, and we increase the value of $\tau$ from 0 to observe possible bifurcations. In this situation, the characteristic equation becomes

where

Let $\lambda = i\omega$ be a purely imaginary root of Eq (3.6). Then we have

From the two equations above, we have

where $v = \omega ^{2}, \bar P_{i} = \bar P_{i}(\tau)(i = 0, 1, 2, 3, 4, 5),$

We choose $\theta(\tau)\in [0, 2\pi)$, and $\sin\theta(\tau)$ and $\cos\theta(\tau)$ are given by Eq (3.7). Similar to Case Ⅰ, there exists a set $\bar I \subset (0, \tau_{1\text{max}})$ such that, when $\tau \in \bar I$, Eq (3.8) has a positive real root. Then we have $S_{n}(\tau) = \tau-\frac{\theta(\tau)+2n\pi}{\omega(\tau)}$, where $n\in\mathbb{N}_{0}$ and $\tau\in \bar I$.

Theorem 3.3. Suppose that $\tau_{1} = \tau_{2} = \tau$, $R_{0} > 1$, and condition $(H_{1})$ hold. When $\tau = \tau^{*}\in \bar I$ such that $S_{n}(\tau^{*}) = 0$ for some $n\in \mathbb{N}_{0}$ hold, the characteristic Eq (3.6) has a pair of purely imaginary roots $\pm i\omega (\tau^{*})$. If $\Delta(\tau^{*}) > 0$ ($< 0$), then the roots $\pm i\omega (\tau^{*})$ cross the imaginary axis from left to right (from right to left), where

Theorem 3.4. Suppose that $\tau_{1} = \tau_{2} = \tau$ and $R_{0} > 1$. If condition $(H_{1})$ holds, the following conclusions hold.

(i) When the function $S_{0}(\tau)$ has no positive roots for $\tau\in \bar I$, the positive equilibrium $E^{*}$ is locally asymptotically stable for $\tau\in \bar I$;

(ii) When the function $S_{n}(\tau)$ has roots such that $\tau^{0} < \tau^{1} < \cdots < \tau^{m}\in \bar I$ and $S^{'}_{n}(\tau^{j})\ne 0 (j = 0, 1, 2, \cdots, m)$, then the positive equilibrium $E^{*}$ is locally asymptotically stable on $\left[0, \tau^{0}\right)\cup \left(\tau^{m}, \tau_{1max} \right) \cap \bar I$.

4.   Global attractivity of the positive equilibrium

Similar to the method used to prove the uniform persistence in reference ^[40], we have the following result.

Theorem 4.1. If $R_{0} > 1$, then model (2.1) is uniformly persistent, and the solution $(x(t), y(t), z(t))^{T}$ of model (2.1) with initial condition (2.2) satisfies

where $\rho > 0$ and $d > 0$ satisfy $q \equiv \frac {1}{a_{1}\rho y^{*}+\frac{a_{2}c_{1}e^{-\delta _{2}\tau _{2}}\rho y^{*}}{D_{3}}+D_{1}} > x^{*}\; and\; x^{\Delta}\equiv q(1-e^{-\frac {d}{q}}) > x^{*}.$

Next, similar to the method used to prove the global stability of equilibria using Barbalat's lemma in references ^[47,48], we consider the global stability of the positive equilibrium $E^*$ when $R_{0} > 1$.

For convenience of description, let $\epsilon$ be any sufficiently small positive number such that $0 < \epsilon < \min\{\nu_{1}, \nu_{2}, \nu_{3}\}$. We define

We define a real symmetric matrix

where

Then we have the following result.

Theorem 4.2. If $R_{0} > 1$ and the symmetric matrix $J(0)$ is positive definite, then the positive equilibrium $E^*$ is globally attractive in $\mathbb C_{1}: = \{\phi \in \mathbb C^{+}\; :\; \phi(0) > 0\}$.

Proof. If $R_{0} > 1$, then there exists a sufficiently small $\epsilon$ such that $0 < \epsilon < \min\{\nu_{1}, \nu_{2}, \nu_{3}\}$ and $J(\epsilon)$ is a positive definite matrix. For this chosen $\epsilon$, there exists a sufficiently large $T(\epsilon) > \hat{\tau}$ such that when $t\ge T(\epsilon)$,

We define

When $t\ge T(\epsilon)$, we have

where

Since $E^{*}$ is an equilibrium of model (2.1), $\dot x(t)$ and $\dot y(t)$ can be rewritten as

When $t\ge T(\epsilon)+\hat\tau$, we have

and

Thus, when $t\ge T(\epsilon)+\hat\tau$, we have

When $t\ge T(\epsilon)+\hat\tau$, we have

and

Thus, when $t\ge T(\epsilon)+\hat\tau$, we have

and

Let $g(x) = x - 1 - \ln(x)$. It is evident that $g(x)\ge 0(x > 0)$, and $g(x) = 0$ if and only if $x = 1$. For $t\ge T(\epsilon)+\hat{\tau}$, we define the function

When $t\ge T(\epsilon)$, we have

where

When $t\ge T(\epsilon)+\hat\tau$, we have

and

Thus, when $t\ge T(\epsilon)+\hat\tau$, we have

When $t\ge T(\epsilon)+\hat\tau$, we have

and

Thus, when $t\ge T(\epsilon)+\hat\tau$, we have

and

We define

When $t\ge T(\epsilon)+\hat\tau$, we have

where

and

When $t\ge T(\epsilon)+\hat\tau$, we define

and

Thus, when $t\ge T(\epsilon)+\hat\tau$, we have

Finally, we define

When $t\ge T(\epsilon)+\hat\tau$, we define

Since $J(\epsilon)$ is positive definite, by applying the Barbalat's lemma ^[49], we have

Therefore, the positive equilibrium $E^{*}$ is globally attractive. □

Remark 4.1. Suppose that $\tau_{1} = \tau_{2} = 0$ and $R_{0} > 1$. If we choose $\theta_{1} = \frac{a_{1}D_{2}+a_{1}D_{1}}{b_{1}^{2}}$ and $\theta_{2} = \frac{a_{1}a_{2}x^{*}}{b_{1}c_{1}}$, then matrix $J(0)$ in our Theorem 4.2 becomes $J$, where $J$ can be found in the proof of Theorem 3.1 in reference ^[35]. Therefore, our Theorem 4.2 extends Theorem 3.1 in reference ^[35].

5.   Numerical simulations

In this section, we run simulations in two parts.

Part Ⅰ. In this part, we give some numerical examples to summarize the applications of the main results. These include the global stability of the boundary equilibrium $E_{0}$, the local stability, and global attractivity of the positive equilibrium $E^{*}$, and the existence of a Hopf bifurcation.

First, we determine the values of the parameters $a_{1}$, $a_{2}$, $D_{1}$, and $D_{3}$ to be $a_{1} = 1$, $a_{2} = 2$, $D_{1} = 1.01$, and $D_{3} = 1.02$.

Then, we choose $b_{1} = 5$, $c_{1} = 5$, $c_{2} = 15$, $D_{2} = 8.01$, $\delta_{1} = 2$, $\tau_{1} = 2$, $\delta_{2} = 2$, and $\tau_{2} = 2$. These values satisfy the condition $R_{0} = 0.0113 < 1$. According to Theorem 2.1 and Figure 1, the boundary equilibrium $E_{0}(0.9901, 0, 0)$ is globally asymptotically stable.

Next, we choose $b_{1} = 5$, $c_{1} = 15$, $c_{2} = 5$, $D_{2} = 4.01$, and $\delta_{2} = 0.1$. Under these parameters, it can be checked that Eq (3.3) has only one positive real root. Based on Theorem 3.1, the picture of $S_{0}(\tau_{2})$ can be drawn clearly. It follows from Figure 2 that there exists 2 roots denoted by $\tau_{2}^{0} = 1.5704$ and $\tau_{2}^{1} = 17.6442$. Then we have $\omega(\tau_{2}^{0}) = 0.6661$ and $\omega(\tau_{2}^{1}) = 0.1547$, and the transversality condition is determined by the sign of $S_{0}^{'}(\tau_{2})$.

Additionally, when $\tau_{2}\in \left[ 0, \tau_{2}^{0}\right)$, the positive equilibrium $E^{*}$ is locally asymptotically stable. For example (see Figure 3), when $\tau_{2} = 1\in \left[ 0, \tau_{2}^{0}\right)$, we have the positive equilibrium $E^{*} = (0.8020, 0.0370, 0.0999)$. As the time delay increases, a Hopf bifurcation occurs at $\tau_{2} = \tau_{2}^{0}$, and the positive equilibrium $E^{*}$ loses stability for $\tau_{2}\in (\tau_{2}^{0}, \tau_{2}^{1})$. For example (see Figure 4), when $\tau_{2} = 14\in (\tau_{2}^{0}, \tau_{2}^{1})$. As $\tau_{2}$ continues to increase, when $\tau_{2}\in(\tau_{2}^{1}, +\infty)\cap I$, the positive equilibrium $E^{*}$ is locally asymptotically stable. For example (see Figure 5), when $\tau_{2} = 17.79\in(\tau_{2}^{1}, +\infty)\cap I$, we have the positive equilibrium $E^{*} = (0.8020, 0.1179, 0.0595)$.

Next, we choose $b_{1} = 16.1$, $c_{1} = 3.1$, $c_{2} = 14.1$, $D_{2} = 7.1$, and $\delta_{2} = 0.1$. Under these parameters, it can be checked that Eq (3.3) has two positive real roots. The picture of $S_{0}(\tau_{2})$ follows from Figure 6, and there exist 2 roots denoted by $\tau_{22}^{0} = 0.6617$ and $\tau_{21}^{0} = 1.2539$. Then we have $\omega(\tau_{22}^{0}) = 2.0852$ and $\omega(\tau_{21}^{0}) = 1.6120$.

Additionally, when $\tau_{2}\in \left[ 0, \tau_{22}^{0}\right)$, the positive equilibrium $E^{*}$ is locally asymptotically stable. For example (see Figure 7), when $\tau_{2} = 0.5\in \left[ 0, \tau_{22}^{0}\right)$, we have the positive equilibrium $E^{*} = (0.4410, 0.6930, 0.2823)$. As the time delay increases, a Hopf bifurcation occurs at $\tau_{2} = \tau_{22}^{0}$, and the positive equilibrium $E^{*}$ loses stability for $\tau_{2}\in (\tau_{22}^{0}, \tau_{21}^{0})$. For example (see Figure 8), when $\tau_{2} = 1.2\in (\tau_{22}^{0}, \tau_{21}^{0})$. As $\tau_{2}$ continues to increase, when $\tau_{2}\in(\tau_{21}^{0}, +\infty)\cap I$, the positive equilibrium $E^{*}$ is locally asymptotically stable. For example (see Figure 9), when $\tau_{2} = 1.35\in(\tau_{21}^{0}, +\infty)\cap I$, we have the positive equilibrium $E^{*} = (0.4410, 0.7193, 0.2692)$.

Lastly, based on reference ^[48], Theorem 4.2 gives some sufficient conditions for the global attractivity of the positive equilibrium $E^{*}$ in model (2.1) with time delays. For example (see Figure 10), let us choose the parameter values $a_{1} = 0.8890$, $a_{2} = 0.0970$, $b_{1} = 37.2064$, $c_{1} = 6.1465$, $c_{2} = 1.0287$, $D_{1} = 1.0300$, $D_{2} = 1.0100$, $D_{3} = 19.9259$, $\delta_{1} = 1$, $\tau_{1} = {3.3007e-6}$, $\delta_{2} = 1$, and $\tau_{2} = {4.0184e-6}$. Then we have $R_{0} = 36.4839 > 1$, and model (2.1) has a positive equilibrium $E^{*} = (0.0271, 40.2768, 0.0011)$. By selecting $\theta_{1} = 0.9601$ and $\theta_{2} = 0.0018$, it can be easily verified that $J(0)$ is positive definite. Therefore, the conditions of Theorem 4.2 are satisfied, and the positive equilibrium $E^{*}$ is globally attractive. Without changing any other parameters, when we increase $\tau_{1}$ to $\tau_{1} = {6.3000e-6}$ and perform calculations, it does not satisfy the conditions of Theorem 4.2. However, through numerical simulations, it is observed that under the set of parameters, the positive equilibrium $E^{*}$ is globally attractive. Theorem 4.2 only provides a sufficient condition for the global attractivity of the positive equilibrium $E^{*}$ and is somewhat conservative.

Furthermore, we choose the parameter values $a_{1} = 1.8883$, $a_{2} = 0.6644$, $b_{1} = 40.0752$, $c_{1} = {6.1465e-4}$, $c_{2} = 98.1660$, $D_{1} = 1.0300$, $D_{2} = 1.0100$, $D_{3} = 14.1512$, $\delta_{1} = 1$, $\delta_{2} = 1$, $\theta_{1} = 0.2429$, and $\theta_{2} = {7.4998e-4}$. When $\tau_{1} = 0$, for any $\tau_{2} > 0$, the positive equilibrium $E^{*}$ is locally asymptotically stable. Using a genetic algorithm and MATLAB simulation, we determine the global attractivity region of model (2.1) with the specified parameter values, as shown by the blue dashed region in Figure 11(a). Further, we can observe that under the parameter values, $\tau_{2}$ has a relatively small impact on the global attractivity compared to $\tau_{1}$. To ensure that the equilibrium is globally attractive, we can observe from Figure 11 that the range of variation in $\tau_{1}$ is much smaller compared to the range of variation in $\tau_{2}$. Therefore, even small changes in $\tau_{1}$ can lead to instability of the equilibrium, while changes in $\tau_{2}$ seem to have a minor impact. From a biological perspective, the time delay ($\tau_{1}$) that the organism (Sphingomonas sp.) stores the nutrient (MCs) is more significant than the time delay ($\tau_{2}$) it takes for Sphingomonas sp. to produce degrading enzymes for the sustained degradation of MCs. When $\tau_{1} = \tau_{2} = \tau$, within the range $\tau\in\left[0, 3.65\right)$, the positive equilibrium $E^{*}$ is locally asymptotically stable, as represented by the red line in Figure 11(b). Therefore, the red dashed line falling within the blue region indicates that the positive equilibrium $E^{*}$ is globally asymptotically stable.

Part Ⅱ. In this part, based on the MCs degradation experiments by Sphingopyxis sp. USTB-05 and enzymes of USTB-05, as described in reference ^[43], we modify model (1.2) into the following two models:

The biological meanings of the parameters in models (5.1) and (5.2) are provided in Tables 2–4. The experimental data for Figures 2 and 4 in reference ^[43] are provided in Table 5. However, since reference ^[43] has not provided the accurate experimental data for 48 hours in Table 5(a) and 10 hours in Table 5(b), and from a mathematical point of view, it is unlikely that the experimental values at these two-time points are zero. Therefore, for simulations, we only use the first four experimental data sets as the parameter values.

Firstly, according to the experimental procedure described in reference ^[43], we set the initial concentration of USTB-05 to be $10 ml/L$ and the initial concentration of enzymes of USTB-05 to be $100 ml/L$. Furthermore, using the experimental data of Table 5 and the least squares method (LSM), the remaining parameter values in models (5.1) and (5.2) are determined as shown in Tables 2–4. Since $\delta_{1}$ is related to $d_{2}$, for simplicity let us assume $\delta_{1} = d_{2}$. Based on the parameter values provided in Tables 2–5, we have Figures 12 and 13. These figures demonstrate that under the condition where MCs are the sole carbon and nitrogen source, models (5.1) and (5.2) fit well with the MCs degradation experiments by USTB-05 and enzymes of USTB-05.

Furthermore, Table 6 provides the fitted values of MCs in Figures 12 and 13. We use mean squared error (MSE) and root mean squared error (RMSE) to evaluate the reliability of the data fitting for models (5.1) and (5.2) (for example, see ^[50,51]). For convenience, let us introduce the definitions of MSE and RMSE as follows:

Here, $y_{i}$ and $y'_{i}$ represent the experimental data and fitted data, respectively. The positive integer $n$ represents the number of fitted data points. Using the data from Table 6, we obtain Table 7. Since the initial values of the experimental data in Table 5 and fitted data are consistent, they are not included in the evaluation. According to Table 7, we can observe that for model (5.1), the fitted results of model (5.1) with time delay are always better than those of model (5.1) without time delay. From a biological point of view, model (5.1) with time delay is more reasonable.

Based on the simulation results for the experiments described in reference ^[43], it can be observed that when the degradation of MCs by USTB-05 is at 48 hours, the remaining amounts of MC-RR, MC-LR, and MC-YR are approximately $6.93 mg/L$, $1.74 mg/L$, and $2.18 mg/L$, respectively. Furthermore, USTB-05 degrades MC-RR, MC-LR, and MC-YR to $1\%$ of their initial values at approximately 57.5 hours, 54.2 hours, and 78.6 hours, respectively.

6.   Conclusions

This paper is mainly based on references ^[34,40]. We proposed and analyzed the time-delayed model (2.1). Model (2.1) describes the biodegradation process of MCs by Sphingomonas sp. and its degrading enzyme. Theorem 2.1 provides the condition ($R_{0} < 1$) for global asymptotic stability of the boundary equilibrium $E_{0}$ in model (2.1). It implies that, in the chemostat, Sphingomonas sp. cannot sustainably degrade MCs, when $R_{0} < 1$.

Theorems 3.1–3.4 provide some conditions for local asymptotic stability of the positive equilibrium $E^{*}$ and the existence of Hopf bifurcation in model (2.1). According to Theorem 3.1 and the results of simulations, when $\tau_{2} < \tau_{2}^{0}$, the positive equilibrium $E^{*}$ is locally asymptotically stable. When $\tau_{2} \in(\tau_{2}^{0}, \tau_{2}^{m})$, the positive equilibrium $E^{*}$ may becomes unstable. When $\tau_{2} > \tau_{2}^{m}$, the positive equilibrium $E^{*}$ becomes stable again. This indicates that time delay $\tau_{2}$ has a significant impact on the stability of the positive equilibrium $E^{*}$. From a biological point of view, it means time delay ($\tau_{2}$) in the production of degrading enzymes by Sphingomonas sp. makes it more difficult to stably degrade MCs.

Theorem 4.1 provides the condition ($R_{0} < 1$) for uniform persistence in model (2.1). From a biological point of view, under some conditions, the biodegradation of MCs is sustainable. Time delays $\tau_{1}$ and $\tau_{2}$ do not affect the sustained degradation of MCs. Theorem 4.2 provides some sufficient conditions for global attractivity of the positive equilibrium $E^{*}$ in model (2.1) by constructing appropriate Lyapunov functionals and analyzing inequalities. Considering numerical simulations, from a biological perspective, the time delay ($\tau_{1}$) that the organism (Sphingomonas sp.) stores the nutrient (MCs) is more significant than the time delay ($\tau_{2}$) it takes for Sphingomonas sp. to produce degrading enzymes for the sustained degradation of MCs.

Finally, based on experimental data from reference ^[43], a least squares fitting is performed, and the fitting results demonstrate the rationality of models (5.1) and (5.2) to some extent. Furthermore, as future work to make the model (2.1) more biologically plausible, the addition of some diffusion terms may have significant theoretical and practical implications.

Author contributions

L.Z.: Writing–original draft, Methodology, Software, Writing–review and editing; M.L.: Writing–original draft, Methodology, Software, Writing–review and editing; W.M.: Writing–original draft, Methodology, Writing–review and editing, Funding acquisition, Supervision. All authors have read and agreed to the published version of the manuscript.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors would like to thank the reviewers and the editors for their careful reading, helpful comments, and suggestions that greatly improved the paper. This paper has been supported by the National Natural Science Foundation of China (No. 12371481) and the Beijing Natural Science Fundation (No. 1202019).

Conflict of interest

The authors declare there is no conflict of interest.