Loading [MathJax]/jax/output/SVG/jax.js
Research article

Assessing the therapeutic response of tumors to hypoxia-targeted prodrugs with an in silico approach


  • Received: 06 May 2022 Revised: 16 June 2022 Accepted: 29 June 2022 Published: 01 August 2022
  • Tumor hypoxia is commonly recognized as a condition stimulating the progress of the aggressive phenotype of tumor cells. Hypoxic tumor cells inhibit the delivery of cytotoxic drugs, causing hypoxic areas to receive insufficient amounts of anticancer agents, which results in adverse treatment responses. Being such an obstruction to conventional therapies for cancer, hypoxia might be considered a target to facilitate the efficacy of treatments in the resistive environment of tumor sites. In this regard, benefiting from prodrugs that selectively target hypoxic regions remains an effective approach. Additionally, combining hypoxia-activated prodrugs (HAPs) with conventional chemotherapeutic drugs has been used as a promising strategy to eradicate hypoxic cells. However, determining the appropriate sequencing and scheduling of the combination therapy is also of great importance in obtaining favorable results in anticancer therapy. Here, benefiting from a modeling approach, we study the efficacy of HAPs in combination with chemotherapeutic drugs on tumor growth and the treatment response. Different treatment schedules have been investigated to see the importance of determining the optimal schedule in combination therapy. The effectiveness of HAPs in varying hypoxic conditions has also been explored in the study. The model provides qualitative conclusions about the treatment response, as the maximal benefit is obtained from combination therapy with greater cell death for highly hypoxic tumors. It has also been observed that the antitumor effects of HAPs show a hypoxia-dependent profile.

    Citation: Defne Yilmaz, Mert Tuzer, Mehmet Burcin Unlu. Assessing the therapeutic response of tumors to hypoxia-targeted prodrugs with an in silico approach[J]. Mathematical Biosciences and Engineering, 2022, 19(11): 10941-10962. doi: 10.3934/mbe.2022511

    Related Papers:

    [1] Xinyou Meng, Jie Li . Stability and Hopf bifurcation analysis of a delayed phytoplankton-zooplankton model with Allee effect and linear harvesting. Mathematical Biosciences and Engineering, 2020, 17(3): 1973-2002. doi: 10.3934/mbe.2020105
    [2] Zhichao Jiang, Xiaohua Bi, Tongqian Zhang, B.G. Sampath Aruna Pradeep . Global Hopf bifurcation of a delayed phytoplankton-zooplankton system considering toxin producing effect and delay dependent coefficient. Mathematical Biosciences and Engineering, 2019, 16(5): 3807-3829. doi: 10.3934/mbe.2019188
    [3] Juan Li, Yongzhong Song, Hui Wan, Huaiping Zhu . Dynamical analysis of a toxin-producing phytoplankton-zooplankton model with refuge. Mathematical Biosciences and Engineering, 2017, 14(2): 529-557. doi: 10.3934/mbe.2017032
    [4] Zhenyao Sun, Da Song, Meng Fan . Dynamics of a stoichiometric phytoplankton-zooplankton model with season-driven light intensity. Mathematical Biosciences and Engineering, 2024, 21(8): 6870-6897. doi: 10.3934/mbe.2024301
    [5] Elvira Barbera, Giancarlo Consolo, Giovanna Valenti . A two or three compartments hyperbolic reaction-diffusion model for the aquatic food chain. Mathematical Biosciences and Engineering, 2015, 12(3): 451-472. doi: 10.3934/mbe.2015.12.451
    [6] Xin Zhao, Lijun Wang, Pankaj Kumar Tiwari, He Liu, Yi Wang, Jianbing Li, Min Zhao, Chuanjun Dai, Qing Guo . Investigation of a nutrient-plankton model with stochastic fluctuation and impulsive control. Mathematical Biosciences and Engineering, 2023, 20(8): 15496-15523. doi: 10.3934/mbe.2023692
    [7] Xiaomeng Ma, Zhanbing Bai, Sujing Sun . Stability and bifurcation control for a fractional-order chemostat model with time delays and incommensurate orders. Mathematical Biosciences and Engineering, 2023, 20(1): 437-455. doi: 10.3934/mbe.2023020
    [8] Yuan Ma, Yunxian Dai . Stability and Hopf bifurcation analysis of a fractional-order ring-hub structure neural network with delays under parameters delay feedback control. Mathematical Biosciences and Engineering, 2023, 20(11): 20093-20115. doi: 10.3934/mbe.2023890
    [9] Saswati Biswas, Pankaj Kumar Tiwari, Yun Kang, Samares Pal . Effects of zooplankton selectivity on phytoplankton in an ecosystem affected by free-viruses and environmental toxins. Mathematical Biosciences and Engineering, 2020, 17(2): 1272-1317. doi: 10.3934/mbe.2020065
    [10] Wenjie Yang, Qianqian Zheng, Jianwei Shen, Linan Guan . Bifurcation and pattern dynamics in the nutrient-plankton network. Mathematical Biosciences and Engineering, 2023, 20(12): 21337-21358. doi: 10.3934/mbe.2023944
  • Tumor hypoxia is commonly recognized as a condition stimulating the progress of the aggressive phenotype of tumor cells. Hypoxic tumor cells inhibit the delivery of cytotoxic drugs, causing hypoxic areas to receive insufficient amounts of anticancer agents, which results in adverse treatment responses. Being such an obstruction to conventional therapies for cancer, hypoxia might be considered a target to facilitate the efficacy of treatments in the resistive environment of tumor sites. In this regard, benefiting from prodrugs that selectively target hypoxic regions remains an effective approach. Additionally, combining hypoxia-activated prodrugs (HAPs) with conventional chemotherapeutic drugs has been used as a promising strategy to eradicate hypoxic cells. However, determining the appropriate sequencing and scheduling of the combination therapy is also of great importance in obtaining favorable results in anticancer therapy. Here, benefiting from a modeling approach, we study the efficacy of HAPs in combination with chemotherapeutic drugs on tumor growth and the treatment response. Different treatment schedules have been investigated to see the importance of determining the optimal schedule in combination therapy. The effectiveness of HAPs in varying hypoxic conditions has also been explored in the study. The model provides qualitative conclusions about the treatment response, as the maximal benefit is obtained from combination therapy with greater cell death for highly hypoxic tumors. It has also been observed that the antitumor effects of HAPs show a hypoxia-dependent profile.



    Plankton includes plants and animals that float freely in some fresh water bodies, and almost all aquatic life is based on plankton [1]. Aquatic ecosystems are affected by many factors, including physical and chemical signals in the environment, plankton predation and competition [2,3]. Many scholars have carried out analysis concerning the impact of environment on the ecosystem and the treatment of sewage [4,5,6]. We also know how important plankton itself is to the wealth of marine ecosystems and ultimately to the planet itself. On the one hand, plankton species have positive effects on the environment, such as providing food for marine life, oxygen for animal life; on the other hand, it have harmful effects, such as economic losses to fisheries and tourism due to algae blooms [7,8].

    In recent years, different models of plankton have been established and studied, for example, model with two harmful phytoplankton [9], models with time delays [10,11] and stochastic models [12,13,14]. Toxins produced by harmful phytoplankton tend to be concentrated at higher levels in the food web, as they can spread through the marine food web, affecting herbivores at higher nutrient levels, reaching fish, and through them eventually reaching marine mammals, even in seabirds [9,11]. There is also some evidence that the occurrence of toxin-producing phytoplankton is not necessarily harmful, but rather helps maintain a stable balance of nutrient dynamics through the coexistence of all species. These results suggest that toxin-producing phytoplankton (TPP) play an important role in the growth of zooplankton populations [15].

    It was shown that aquatic plant systems not only have extraordinary memories of climatic events, but also exhibit phenomenological responses based on memory [16,17,18]. The authors noted that environmental factors often alter the expression of chromatin in multiple responsive genes in [19]. Environment-induced chromatin markers are at certain sites and is transmitted by cell division, allowing plants to acquire memories of environmental experiences. This ensures that the plant can adapt to changes in its environment or perform better if the event occurs again. In some cases, it is passed on to the next generation, namely, epigenetic mechanisms. This mechanism is crucial for plants' stress memory and adaptation to the environment, suggesting that plants do form memory and defense mechanisms in certain environments. In addition, a large amount of zooplankton chemical signal learning and corresponding reactions have been documented for aquatic systems [20,21]. In summary, such memory and genetic characteristics can not be neglected for plankton systems.

    As we all know, fractional order derivatives are a good tool for describing the memory and genetic properties of various materials and processes. In other words, the application of fractional order dynamical systems can fully reflect some long-term memory and non-local effects. That is, fractional differential equations have an advantage over classical integer differential equations for describing such systems. In recent years, more and more researchers began to study the qualitative theory and numerical solution of fractional order biological model [22,23,24]. The main reason is that fractional order equations are naturally related to memory systems that exist in most biological systems [25,26]. In addition, fractional-order derivative has also been widely studied and applied in physics[27], engineering [28], biology [29] and many other fields [30,31,32]. At present, there are more than six definitions of fractional derivative, among which Riemann-Liouville and Caputo derivatives are the most commonly used [33]. In the case of fractional Caputo derivative, the initial conditions are expressed by the values of the unknown function and its integer derivative with clear physical meaning [34]. So we will adapt the Caputo's definition in our paper.

    The interactions between phytoplankton and zooplankton do not occur instantaneously in real ecosystems. Instead, the response of zooplankton to contacts with phytoplankton is likely to be delayed due to gestation. For example, in [35], the authors discussed Hopf bifurcation in the presence of time delay required for toxin-phytoplankton maturation. The universality of time-delay coupled system indicates its importance, applicability and practicability in a wide range of biological systems [36,37]. In fact, time delay may change the qualitative behavior of dynamic system [38,39].

    In [40], the authors considered a fractional nutrient-phytoplankton-zooplankton system as follows

    {DαX(t)=x0aXb1XY+c1Y+c2Z,DαY(t)=b2XYd1YZe+Yc3Y,DαZ(t)=d2YZe+YfYZc4Z. (1.1)

    Based on the above model, we classify phytoplankton into non-toxic phytoplankton and toxic phytoplankton and put forward an improved fractional order four-dimensional ecological epidemiological model with delay. The system is established as follows:

    {DαX(t)=ΛμX(t)b1X(t)Y1(t)b2X(t)Y2(t)+c1Y1(t)+c2Y2(t)+c3Z(t),DαY1(t)=k1b1X(t)Y1(t)η1Y1(t)Z(t)h1Y1(t)Y2(t)μ1Y1(t),DαY2(t)=k2b2X(t)Y2(t)η2Y2(t)Z(t)h2Y1(t)Y2(t)μ2Y2(t),DαZ(t)=θ1η1Y1(tτ)Z(tτ)+θ2η2Y2(tτ)Z(tτ)δY2(t)Z(t)μ3Z(t), (1.2)

    subjected to the biologically feasible initial condition:

    X(0)0, Y1(t)=ϕ(t)0, Y2(t)=ψ(t)0, Z(t)=ζ(t)0,t[τ,0], (1.3)

    where ϕ(t), ψ(t) and ζ(t) are continuous function defined on t[τ,0].

    The meaning of state variables and parameters are listed in Table 1; Dα (0<α<1) denotes Caputo fractional differential operator, and the model are based on the following scenarios:

    Table 1.  Description of state variables and parameters in the system (1.2).
    variables Descriptions
    X(t) concentration of nutrient population at time t
    Y1(t) concentration of phytoplankton population at time t
    Y2(t) concentration of toxic phytoplankton population at time t
    Z(t) concentration of zooplankton population at time t
    Parameters Descriptions Default value
    Λ Constant input of nutrient [0.5, 3]
    b1 Nutrient uptake rate for the phytoplankton population [0, 3]
    b2 Nutrient uptake rate for the toxic phytoplankton population [0, 3]
    k1 Nutrient-phytoplankton conversion rate (0, 1)
    k2 Nutrient-toxic phytoplankton conversion rate (0, 1)
    c1 Nutrient recycling rate after the death of phytoplankton (0, 0.1)
    c2 Nutrient recycling rate after the death of toxic phytoplankton (0, 0.1)
    c3 Nutrient recycling rate after the death of zooplankton (0, 0.7)
    η1 Maximal zooplankton ingestion rate (0, 3)
    η2 Maximal zooplankton ingestion rate (0, 2.5)
    θ1 Maximal phytoplankton-zooplankton conversion rate (0, 0.7)
    θ2 Maximal toxic phytoplankton-zooplankton conversion rate (0, 0.8)
    μ Rate of nutrient loss (0, 1.5)
    μ1 Phytoplankton mortality rate (0, 1)
    μ2 Toxic phytoplankton mortality rate (0, 1)
    μ3 Zooplankton death rate (0, 0.6)
    δ Rate of zooplankton decay due to toxin producing phytoplankton (0, 0.2)
    h1 competition effect for phytoplankton (0, 0.3)
    h2 competition effect for toxic phytoplankton (0, 0.3)

     | Show Table
    DownLoad: CSV

    (H1) X(t), Y1(t), Y2(t) and Z(t) represent nutrient population, phytoplankton population, toxic phytoplankton population and zooplankton population, respectively.

    (H2) In real ecosystems, phytoplankton compete with each other for essential resources: nutrients and light. So as the model in the [41], we assume that, for nutrient X(t), phytoplankton population Y1(t) is in competition with toxic phytoplankton population Y2(t), h1 and h2 represent the influence on Y1(t) and Y2(t) in the competition, respectively.

    (H3) Zooplankton do not grow instantaneously after consuming phytoplankton, and pregnancy of predators requires a discrete time delay τ.

    (H4) Zooplankton populations feed only on phytoplankton, and only some of the dead phytoplankton and zooplankton are recycled into nutrients.

    (H5) The toxic phytoplankton has both positive and negative effect on zooplankton, corresponding to the term θ2η2Y2(tτ)Z(tτ) and δY2(t)Z(t) in the last equation of the system(1.2).

    The above assumption (H5) is based on the result in [42]. In fact, the authors concluded that the food selection mechanism of plankton may not yet mature. When different algae of the same species (toxic and non-toxic) coexist, some zooplankton may have poor ability to select between toxic and non-toxic algae, and even show a slight preference for toxic strains.

    The present paper is organized as follows. In section 2, some preliminaries are presented. In section 3, qualitative analysis of the system is performed. In section 4, some numerical examples and simulations are exploited to verify the theoretical results. In the last section, some conclusions and discussions are provided.

    For convenience, we list some of the basic definitions and lemmas of the fractional calculus. In fractional-order calculus, there are many fractional-order integration and fractional-order differentiation that have been defined, for example, the Grunwald-Letnikov (GL) definition, the Riemann-Liouville (RL) definition and the Caputo definition. Since the initial condition is the same as the form of integral differential equation, we will adopt the definition of Caputo in this paper.

    Definition 2.1. [34] The Riemann-Liouville fractional integral of order α>0 for a function f:R+R is defined by

    0Dαtf(t)=1Γ(α)t0(ts)α1f(s)ds, t0.

    Based on this definition of Riemann-Liouville fractional integral, the fractional-order derivative in Riemann-Liouville sense and Caputo sense are given.

    Definition 2.2. [34] The Riemann-Liouville fractional derivative of order α>0 for a function f:R+R is defined by

    RL0Dαtf(t)=dkdtk(0D(kα)tf(t))=1Γ(kα)dkdtkt0(ts)kα1f(s)ds, t0,

    where k1α<k, kN and Γ() is the Gamma function, Γ(α)=+0tα1etdt.

    In particular, when 0<α<1, we have

    RL0Dαtf(t)=1Γ(1α)ddtt0(ts)αf(s)ds.

    Definition 2.3. [34] The Caputo fractional derivative of order α>0 for a function f:R+R is defined by

    C0Dαtf(t)=0D(kα)tf(k)(t)=1Γ(kα)t0(ts)kα1f(k)(s)ds,t0,

    where k1α<k, kN and f(m)(t) is the m-order derivative of f(t). In particular, when 0<α<1, we have

    C0Dαtf(t)=1Γ(1α)t0f(s)(ts)αds.

    Definition 2.4. [34] The two-parameter Mittag-Leffler function is defined by

    Eα,β(z)=+i=0ziΓ(αi+β),α>0,β>0.

    When β=1, the two-parameter Mittag-Leffler function becomes to the one-parameter Mittag- Leffler function, i.e.

    Eα(z)=Eα,1(z)=+i=0ziΓ(αi+1),α>0.

    Theorem 2.5. [43] Consider the following commensurate fractional-order system:

    dαxdtα=f(x),x(0)=x0,

    with 0<α<1 and xRn. The equilibrium points of the above system are calculated by solving the equation: f(x)=0. These points are locally asymptotically stable if all eigenvalues λi of the Jacobian matrix evaluated at the equilibrium points satisfy the inequality: |arg(λi)|>απ2.

    Since the proof of the positivity and boundedness of the solution of the system(1.2) is similar to Theorem 2 and Theorem 3 in the Ref.[38], we will not prove it here.

    The equilibriums of model (1.2) are obtained by solving the following algebraic system

    {ΛμXb1XY1b2XY2+c1Y1+c2Y2+c3Z=0,k1b1XY1η1Y1Zh1Y1Y2μ1Y1=0,k2b2XY2η2Y2Zh2Y1Y2μ2Y2=0,θ1η1Y1Z+θ2η2Y2ZδY2Zμ3Z=0. (3.1)

    By simple calculation, we obtain seven equilibriums of system (1.2), namely:

    (1) E0=(X(0), 0, 0, 0) with X(0)=Λμ.

    (2) E1=(X(1), Y(1)1, 0, 0) with X(1)=μ1k1b1, Y(1)1=μ(1R1)b1(c1R1b1X(0)1), where R1=X(0)X(1). And the feasibility conditions for E1 are simplified as:

    X(0)<X(1)<c1b1orc1b1<X(1)<X(0).

    (3) E2=(X(2), 0, Y(2)2, 0) with X(2)=μ2k2b2, Y(2)2=μ(1R2)b2(c2R2b2X(0)1), where R2=X(0)X(2). And the feasibility conditions for E2 are simplified as:

    X(0)<X(2)<c2b2orc2b2<X(2)<X(0).

    (4) E3=(X(3), Y(3)1, Y(3)2, 0) with Y(3)1=k2b2X(3)μ2h2, Y(3)2=k1b1X(3)μ1h1, and X(3) is uniquely determined by the following equation:

    a1X2+a2X+a3=0, (3.2)

    where

    a1=b1b2(h1k2+h2k1)<0,a2=μh1h2+μ2b1h1+μ1b2h2+c1h1b2k2+c2h2k1b1,a3=Λh1h2μ2c1h1μ1c2h2.

    If Λ>μ2c1h2+μ1c2h1, then Descartes rule of sign ensures that the above Eq.(3.2) possesses a uniquely positive root. And the feasibility conditions for E3 are simplified as:

    Λ>μ2c1h2+μ1c2h1,X(3)>X(1)andX(3)>X(2).

    (5) E4=(X(4), Y(4)1, 0, Z(4)) with

    X(4)=Λη1+c1η1Y(4)1μ1c3μη1+b1η1Y(4)1k1b1c3,Y(4)1=μ3θ1η1,Z(4)=k1b1X(4)μ1η1.

    The feasibility conditions of E4 are simplified as:

    R3=X(4)X(1)>1.

    (6) E5=(X(5), 0, Y(5)2, Z(5)) with

    X(5)=Λη2+c2η2Y(5)2μ2c3μη2+b2η2Y(5)2k2b2c3,Y(5)2=μ3θ2η2δ,Z(5)=k2b2X(5)μ2η2.

    Considering the biological background, we assume θ2η2>δ is reasonable, and the feasibility conditions for E5 are simplified as:

    R4=X(5)X(2)>1.

    (7) E6=(X(6), Y(6)1, Y(6)2, Z(6)) with

    Y(6)1=k2b2X(6)η2Z(6)μ2h2,Y(6)2=k1b1X(6)η1Z(6)μ1h1,

    Z(6)=(θ1η1h1k2b2+θ2η2h2k1b1δh2k1b1)X(6)+δh2μ1θ1η1h1μ2θ2η2h2μ1μ3h1h2η1(θ1η2h1+θ2η2h2δh2). And X(6) is uniquely determined by the following equation:

    b1X2+b2X+b3=0, (3.3)

    where

    b1=(b1h1η2+b2h2η1)(θ1η1h1k2b2+θ2η2h2k1b1δh2k1b1)b1b2η1(h1k2+h2k1)(θ1η2h1+θ2η2h2δh2),b2=(b1h1η2+b2h2η1)(δh2μ1θ1η1h1μ2θ2η2h2μ1μ3h1h2)+(c3h1h2c1h1η2c2h2η1)(θ1η1h1k2b2+θ2η2h2k1b1δh2k1b1)+η1(θ1η2h1+θ2η2h2δh2)(b1h1μ2+b2h2μ1+c1h1k2b2+c2h2k1b1μh1h2),b3=(c3h1h2c1h1η2c2h2η1)(δh2μ1θ1η1h1μ2θ2η2h2μ1μ3h1h2)+(θ1η2h1+θ2η2h2δh2)(Λh1h2c1h1μ2c2h2μ1),

    If b1b3<0, then Descartes rule of sign ensures that the above Eq.(3.3) possesses a uniquely positive root. The feasibility conditions for E6 are simplified as: Y(6)1,Y(6)2,Z(6)>0,b1b3<0.

    Remark 3.1. (1) The necessary conditions for the existence of E3 are R1>1 and R2>1.

    (2) Because of the complexity of computation, we have not obtain the exact formula of positive equilibrium E6.

    In this subsection we discuss the stability of each equilibrium when τ = 0.

    Obviously, the eigenvalues of the Jacobian matrix of system (1.2) at equilibrium E0 are λ1=μ<0, λ2=μ3<0, λ3=k1b1X(0)μ1, λ4=k2b2X(0)μ2, so we get the following result.

    Theorem 3.2. If R1<1 and R2<1, then the disease-free equilibrium E0 is locally asymptotically stable and it is unstable if R1>1 or R2>1.

    The Jacobian matrix of system (1.2) at equilibrium E1 is

    (μb1Y(1)1b1X(1)+c1b2X(1)+c2c3k1b1Y(1)10h1Y(1)1η1Y(1)100k2b2X(1)h2Y(1)1μ20000θ1η1Y(1)1μ3).

    The characteristic equation at the equilibrium E1 is

    [λ(k2b2X(1)h2Y(1)1μ2)][λ+μ3θ1η1Y(1)1][λ2+(μ+b1Y(1)1)λ+k1b1Y(1)1(b1X(1)c1)]=0. (3.4)

    Thus, we get the following result.

    Theorem 3.3. If 1<R1<X(0)b1c1, R2<R1 and μ3>θ1η1Y(1)1, then the equilibrium E1 is locally asymptotically stable.

    The Jacobian matrix of system (1.2) at equilibrium E2 is

    (μb2Y(2)2b1X(2)+c1b2X(2)+c2c30k1b1X(2)h1Y(2)2μ100k2b2Y(2)2h1Y(2)20η2Y(2)2000θ2η2Y(2)2δY(2)2μ3).

    The characteristic equation at the equilibrium E2 is

    [λ(k1b1X(2)h1Y(2)2μ1)][λ+δY(2)2+μ3θ2η2Y(2)2][λ2+(μ+b2Y(2)2)λ+k2b2Y(2)2(b2X(2)c2)]=0. (3.5)

    So we get the following result.

    Theorem 3.4. If 1<R2<X(0)b2c2, R2>R1 and μ3>(θ2η2f)Y(2)2, then the equilibrium E2 is locally asymptotically stable.

    The Jacobian matrix of system (1.2) at equilibrium E3 is

    (μb1Y(3)1b2Y(3)2b1X(3)+c1b2X(3)+c2c3k1b1Y(3)10h1Y(3)1η1Y(3)1k2b2Y(3)2h2Y(3)20η2Y(3)2000θ1η1Y(3)1+θ2η2Y(3)2δY(3)2μ3),

    with the characteristic equation

    [λ(θ1η1Y(3)1+θ2η2Y(3)2δY(3)2μ3)](λ3+A1λ2+A2λ+A3)=0, (3.6)

    where

    A1=μ+b1Y(3)1+b2Y(3)2,A2=k1b1Y(3)1(b1X(3)c1)+k2b2Y(3)2(b2X(3)c2)h1h2Y(3)1Y(3)2,A3=h1h2Y(3)1Y(3)2(μ+b1Y(3)1+b2Y(3)2)h1k2b2Y(3)1Y(3)2(b1X(3)c1)h2k1b1Y(3)1Y(3)2(b1X(3)c1).

    Denote D(P) denote the discriminant of a polynomial

    Q(λ)=λ3+A1λ2+A2λ+A3. Then

    D(Q)=18A1A2A3+(A1A2)24A3A314A3227A33.

    If μ3>θ1η1Y(3)1+θ2η2Y(3)2δY(3)2, in order to discuss the stability of the equilibrium E3, we get the following result by use of the same method as in Ref [44].

    Proposition 3.5. The equilibrium E3 is asymptotically stable if one of the following conditions holds for polynomial Q and D(Q):

    (1)D(Q)>0,A1>0,A3>0 and A1A2>A3.

    (2)D(Q)<0,A10,A20,A30 and α<23.

    The Jacobian matrix of system (1.2) at equilibrium E4 is

    (m11m12m13m14m210m23m2400m3300m42m43m44),

    where

    m11=μb1Y(4)1,m12=b1X(4)+c1,m13=b2X(4)+c2,m14=c3,m21=k1b1Y(4)1,m23=h1Y(4)1m24=η1Y(4)1,m33=k2b2X(4)η2Z(4)h2Y(4)1μ2,m42=θ1η1Z(4),m43=θ2η2Z(4)δZ(4),m44=θ1η1Y(4)1μ3.

    The characteristic equation at the equilibrium E4 is

    [λ(k2b2X(4)η2Z(4)h2Y(4)1μ2)][λ3+B1λ2+B2λ+B3]=0, (3.7)

    where

    B1=m11m44,B2=m11m44m12m21m24m42,B3=m11m24m42+m12m21m44m14m21m42.

    Denote D(P) denote the discriminant of a polynomial

    Q(λ)=λ3+B1λ2+B2λ+B3. Then

    D(Q)=18B1B2B3+(B1B2)24B3B314B3227B33.

    If μ2>k2b2X(4)η2Z(4)h2Y(4)1, in order to discuss the stability of the equilibrium E4, we get the following result by use of the same method as in Ref [44].

    Proposition 3.6. The equilibrium E4 is asymptotically stable if one of the following conditions holds for polynomial Q and D(Q):

    (1)D(Q)>0,B1>0,B3>0 and B1B2>B3.

    (2)D(Q)<0,B10,B20,B30 and α<23.

    The Jacobian matrix of system (1.2) at equilibrium E5 is

    (ˆm11ˆm12ˆm13ˆm140ˆm2200ˆm31ˆm320ˆm340ˆm42ˆm43ˆm44),

    where

    ˆm11=μb2Y(5)2,ˆm12=b1X(5)+c1,ˆm13=b2X(5)+c2,ˆm14=c3,ˆm22=k1b1X(5)η1Z(5)h1Y(5)2μ1,ˆm31=k2b2Y(5)2ˆm32=h2Y(5)2,ˆm34=η2Y(5)2,ˆm42=θ1η1Z(5),ˆm43=θ2η2Z(5)δZ(5),ˆm44=θ2η2Y(5)2δY(5)2μ3.

    The characteristic equation at the equilibrium E5 is

    [λ(k1b1X(5)η1Z(5)h1Y(5)2μ1)][λ3+C1λ2+C2λ+C3]=0, (3.8)

    where

    C1=ˆm11ˆm44,C2=ˆm11ˆm44ˆm13ˆm31ˆm34ˆm43,C3=ˆm11ˆm34ˆm43+ˆm13ˆm31ˆm44ˆm14ˆm31ˆm43.

    Denote D(P) denote the discriminant of a polynomial

    Q(λ)=λ3+C1λ2+C2λ+C3. Then

    D(Q)=18C1C2C3+(C1C2)24C3C314C3227C33.

    If μ1>k1b1X(5)η1Z(5)h1Y(5)2, in order to discuss the stability of the equilibrium E5, we get the following result by use of the same method as in Ref [44].

    Proposition 3.7. The equilibrium E5 is asymptotically stable if one of the following conditions holds for polynomial Q and D(Q):

    (1)D(Q)>0,C1>0,C3>0 and C1C2>C3.

    (2)D(Q)<0,C10,C20,C30 and α<23.

    Theorem 3.8. The equilibrium E6 is locally asymptotically stable if the following conditions hold:

    (H5) X(6)>max{c1b1, c2b2},

    (H6) k1b1η1<k2b2η2,

    (H7) k3>η2(μ+b1Y(6)1+b2Y(6)2)c3b2,

    (H8) (e1e2e3)e3e4e21>0.

    Proof. The Jacobian matrix of system (1.2) at equilibrium E6 is

    (a11a12a13a14a210a23a24a31a320a340b42a43+b43a44+b44)

    where

    a11=μb1Y(6)1b2Y(6)2,a12=b1X(6)+c1,a13=b2X(6)+c2,a14=c3,a21=k1b1Y(6)1,a23=h1Y(6)1a24=η1Y(6)1,a31=k2b2Y(6)2,a32=h2Y(6)2a34=η2Y(6)2,a43=fZ(6),a44=fY(6)2μ3,b42=θ1η1Z(6),b43=θ2η2Z(6),b44=θ1η1Y(6)1+θ2η2Y(6)2.

    The characteristic equation at the equilibrium E6 is

    λ4+e1λ3+e2λ2+e3λ+e4=0, (3.9)

    where

    e1=a11>0,e2=a13a31a12a21a23a32a34a43a24b42a34b43,e3=a11a23a32a12a23a31a13a21a32+a11a34a43a14a31a43a24a32a43+a11a24b42a14a21b42+a11a34b43a14a31b43a23a34b42a24a32b43,e4=a11a24a32a43+a12a21a34a43a12a24a31a43a14a21a32a43+a11a23a34b42+a11a24a32b43+a12a21a34b43a12a24a31b43a13a21a34b42+a13a24a31b42a14a21a32b43a14a23a31b42.

    By simple calculation, if X(6)>max{c1b1, c2b2}, then e1e2>e3; if k3>η2(μ+b1Y(6)1+b2Y(6)2)c3b2 and k1b1η1<k2b2η2, then e4>0. In summary, the condition of the Routh-Hurwitz criterion above is satisfied for Eq.(3.9), that is,

    e1>0, e1e2>e3, (e1e2e3)e3e4e21>0, e4>0,

    hold. So, all the roots of this Eq.(3.9) have negative real part. This ends the proof.

    Remark 3.9 In Theorem (3.8), (H5)-(H8) is a sufficient condition for equilibrium E6 to be stable, and the necessary and sufficient condition for E6 to be stable is all roots of Eq.(3.9) satisfy |arg(λi)|>απ2.

    In this subsection, according to the research methods in literature [22,38,45], we study the Hopf bifurcation with time delay as the parameter.

    we will analyze the Hopf bifurcation of E6 when τ>0, and the characteristic equation at the equilibrium E6 is

    s4α+p1s3α+p2s2α+p3sα+p4+(q1s3α+q2s2α+q3sα+q4)esτ=0, (3.10)

    where

    p1=a11a44,p2=a11a44a13a31a12a21a23a32a34a43,p3=a11a23a32a12a23a31a13a21a32+a12a21a44+a11a34a43+a13a31a44a14a31a43+a23a32a44a24a32a43,p4=a11a23a32a44+a11a24a32a43+a12a21a34a43+a12a23a31a44a12a24a31a43+a13a21a32a44a14a21a32a43,q1=b44,q2=a11b44a24b42a34b43,q3=a11a24b42+a12a21b44a14a21b42+a11a34b43+a13a31b44a14a31b43+a23a32b44a23a34b42a24a32b43,q4=a11a23a32b44+a11a23a34b42+a11a24a32b43+a12a21a34b43+a12a23a31b44a12a24a31b43+a13a21a32b44a13a21a34b42+a13a24a31b42a14a21a32b43a14a23a31b42.

    Assume that s=iω=ω(cosπ2+isinπ2), ω>0 is a root of Eq.(3.10).

    Substituting s=iω into Eq.(3.10), one gets

    ω4α(cos2απ+isin2απ)+p1ω3α(cos3απ2+isin3απ2)+p2ω2α(cosαπ+isinαπ)+p3ωα(cosαπ2+isinαπ2)+p4+[q1ω3α(cos3απ2+isin3απ2)+q2ω2α(cosαπ+isinαπ)+q3ωα(cosαπ2+isinαπ2)+q4](cosωτisinωτ)=0. (3.11)

    and separating the real and imaginary parts of it, it results in

    {R2cos(ωτ)+I2sin(ωτ)=R1,I2cos(ωτ)R2sin(ωτ)=I1, (3.12)

    Ri,Ii are defined as follows:

    R1=ω4αcos2απ+p1ω3αcos3απ2+p2ω2αcosαπ+p3ωαcosαπ2+p4,R2=q1ω3αcos3απ2+q2ω2αcosαπ+q3ωαcosαπ2+q4,I1=ω4αsin2απ+p1ω3αsin3απ2+p2ω2αsinαπ+p3ωαsinαπ2,I2=q1ω3αsin3απ2+q2ω2αsinαπ+q3ωαsinαπ2.

    It can be acquired from Eq. (3.12) that

    {cos(ωτ)=R1R2+I1I2R22+I22=F(ω),sin(ωτ)=R2I1R1I2R22+I22=G(ω). (3.13)

    Adding the squares of the two equations of Eq.(3.12), we obtain

    ω8α+M+N=0. (3.14)

    where M is a polynomial containing ω7α, ω6α, ω5α, ω4α, ω3α, ω2α, ωα, and N is a constant.

    Define

    h(ω)=ω8α+M+N. (3.15)

    Suppose that N<0. Thus, h(ω) has at least one positive root. The delay τ is regarded as a bifurcation parameter. Let s(ω)=ξ(τ)+iω(τ) be the Eq.(3.10) such that for some initial value of the bifurcation parameter τ0 we have ξ(τ0)=0, ω(τ0)=ω0. Without loss of generality, we assume ω(0)>0. From Eq.(3.13), one can conclude

    τj=1ω0[arccosF(ω)+2jπ],j=0, 1, 2. (3.16)

    where

    τ0=min{τj},j=0, 1, 2.

    To derive the condition of the occurrence for Hopf bifurcation, we have the following Lemma.

    Lemma 3.10. Assume that N<0, then Hopf bifurcation occurs provided h(ω0)0.

    Proof. Differentiating both sides of Eq.(3.10) with respect to τ, it can be obtained that

    (4αs4α1+3αp1s3α1+2αp2s2α1+αp3sα1)dsdτ+(3αq1s3α1+2αq2s2α1+αq3sα1)esτdsdτ+(q1s3α+q2s2α+q3sα+q4)esτ(τdsdτs)=0.

    Hence, one gets

    (dsdτ)1=(4αs4α1+3αp1s3α1+2αp2s2α1+αp3sα1)+(3αq1s3α1+2αq2s2α1+αq3sα1)esτs(q1s3α+q2s2α+q3sα+q4)esττs(4αs4α1+3αp1s3α1+2αp2s2α1+αp3sα1)s(s4α+p1s3α+p2s2α+p3sα+p4)+(3αq1s3α1+2αq2s2α1+αq3sα1)s(q1s3α+q2s2α+q3sα+q4)τs. (3.17)

    Substitute s=iω0 into Eq.(3.17), we have

    Re[(dsdτ)1|τ=τ0]=Re[(4α(iω0)4α1+3αp1(iω0)3α1+2αp2(iω0)2α1+αp3(iω0)α1)(iω0)((iω0)4α+p1(iω0)3α+p2(iω0)2α+p3(iω0)α+p4)+(3αq1(iω0)3α1+2αq2(iω0)2α1+αq3(iω0)α1)(iω0)(q1(iω0)3α+q2(iω0)2α+q3(iω0)α+q4)]=Re[(4α(iω0)4α1+3αp1(iω0)3α1+2αp2(iω0)2α1+αp3(iω0)α1)(iω0)((iω0)4α+p1(iω0)3α+p2(iω0)2α+p3(iω0)α+p4)+(3αq1(iω0)3α1+2αq2(iω0)2α1+αq3(iω0)α1)(iω0)(q1(iω0)3α+q2(iω0)2α+q3(iω0)α+q4)]=h(ω0)2ω0G,

    where

    G=(q1ω3α0cos(3α+1)π2+q2ω2α0cos(2α+1)π2+q3ωα0cos(α+1)π2)2+(q1ω3α0sin(3α+1)π2+q2ω2α0sin(2α+1)π2+q3ωα0sin(α+1)π2+q4)2.

    then sign{dRe(λ)dτ|τ=τ0}=sign{Re[(dλdτ)1|τ=τ0]}=sign{h(ω0)}.

    Obviously, if h(ω0)0 the transversality condition holds, and Hopf bifurcation occurs at τ=τ0.

    Theorem 3.11. Suppose that (H5)-(H8) and N<0 hold, then the positive equilibrium E6 of system (1.2) is asymptotically stable when τ[0,τ0), h(ω0)<0 and unstable when τ>τ0, h(ω0)>0. When τ=τ0, h(ω0)0 a Hopf bifurcation occurs, that is a family of periodic solutions bifurcates from E6 as τ passes through the critical value τ0.

    In this section, some numerical examples are presented to verify the theoretical results. The simulation are based on Adama-Bashforth-Moulton predictor-corrector scheme [46].

    Example 1. For the following set of parameters: Λ= 1, b1= 0.3, b2= 0.25, c1= 0.06, c2= 0.06, c3= 0.06, k1= 0.7, k2= 0.7, η1= 2.1, η2= 0.2, μ= 1, μ1= 0.5, μ2= 0.3, μ3= 0.3, θ1= 0.6, θ2= 0.7, δ= 0.1, h1= 0.2, h2= 0.1.

    In this case R1=0.42<1, R2=0.5833<1. From Figure 1, we can see that the equilibrium E0=(1, 0, 0, 0) is stable for different values of α and different sets of initial values: [X(0), Y1(0), Y2(0), Z(0)]=[0.9, 0.2, 0.2, 0.2], [1.2, 0.5, 0.3, 0.6].

    Figure 1.  (a)-(d) are the time series of the system (1.2), which show that the equilibrium E0 is stable for different values of α (α = 1, 0.9) when τ = 0; the red and and the blue and .. line represents dynamics with initial value [0.9, 0.2, 0.2, 0.2]; the yellow and ... and the green and line represents dynamics with initial value [1.2, 0.5, 0.3, 0.6].

    Example 2. For the following set of parameters: Λ= 0.5, b1= 2.5, b2= 0.3, c1= 0.01, c2= 0.01, c3= 0.06, k1= 0.5, k2= 0.7, η1= 0.2, η2= 2.1, μ= 1, μ1= 0.3, μ2= 0.5, μ3= 0.3, θ1= 0.5, θ2= 0.1, δ= 0.01, h1= 0.2, h2= 0.1.

    In this case, 1<R1=2.0833<X(0)b1c1, R2=0.21<R1 and μ3>θ1η1Y(1)1. From Figure 2, we can see that the equilibrium E1=(0.24, 0.4407, 0, 0) is stable for different values of α and different sets of initial values: [X(0), Y1(0), Y2(0), Z(0)]=[0.1,0.6, 0.5, 0.2], [0.2,0.5, 0.8, 0.6].

    Figure 2.  (a)-(d) are the time series of the system (1.2), which show that the equilibrium E1 is stable for different values of α (α = 1, 0.9) when τ = 0; the red and and the blue and .. line represents dynamics with initial value [0.1, 0.6, 0.5, 0.2]; the yellow and ... and the green and line represents dynamics with initial value [0.2, 0.5, 0.8, 0.6].

    Example 3. For the following set of parameters: Λ= 0.5, b1= 0.3, b2= 2.5, c1= 0.01, c2= 0.01, c3= 0.06, k1= 0.5, k2= 0.6, η1= 2.1, η2= 0.2, μ= 1, μ1= 0.3, μ2= 0.5, μ3= 0.5, θ1= 0.1, θ2= 0.1, δ= 0.01, h1= 0.2, h2= 0.1.

    In this case, 1<R2=1.5<X(0)b2c2, R2>R1 and μ3>(θ2η2f)Y(2)2. From Figure 3, we can see that the equilibrium E2=(0.3333, 0, 0.2024, 0) is stable for different values of α and different sets of initial values: [X(0), Y1(0), Y2(0), Z(0)]=[0.2, 0.1, 0.5, 0.2], [0.1, 0.05, 0.3, 0.6].

    Figure 3.  (a)-(d) are the time series of the system (1.2), which show that the equilibrium E2 is stable for different values of α (α = 1, 0.9) when τ = 0; the red and and the blue and .. line represents dynamics with initial value [0.2, 0.1, 0.5, 0.2]; the yellow and ... and the green and line represents dynamics with initial value [0.1, 0.05, 0.3, 0.6].

    Remark 4.1. The above three examples corresponding to the following ecological interpretation.

    (1) Figure 1 indicates that if R1<1, R2<1, then the phytoplankton can not survive and the zooplankton will also die out. However, this phenomenon usually does not happen in the real world.

    (2) Figure 2 indicates that if R1>1, R1>R2, then the non-toxic phytoplankton will win the competition between phytoplankton and toxic phytoplankton, while the zooplankton will die out due to excessive mortality.

    (3) Figure 3 indicates that if R2>1, R2>R1, then the toxic phytoplankton will win the competition between phytoplankton and toxic phytoplankton, while the zooplankton will die out due to excessive mortality.

    Example 4. For the following set of parameters: Λ= 1.4, b1= 2.8, b2= 2.3, c1= 0.01, c2= 0.01, c3= 0.06, k1= 0.95, k2= 0.1, η1= 2.1, η2= 0.26, μ= 0.2, μ1= 0.2, μ2= 0.8, μ3= 0.5, θ1= 0.2, θ2= 0.5, δ= 0.01, h1= 0.1, h2= 0.1.

    In this case, simple calculation indicates that the sufficient condition (2) of Proposition 3.6 is satisfied. From Figure 4 we can see that the equilibrium E4=(0.4067, 1.1905, 0, 0.4199) is stable for different values of α and different sets of initial values: [X(0), Y1(0), Y2(0), Z(0)]=[0.2, 0.3, 0.5, 0.2], [0.3, 0.4, 0.8, 0.6].

    Figure 4.  (a)-(d) are the time series of the system (1.2), which show that the equilibrium E4 is stable for different values of α (α = 1, 0.95) when τ = 0; the red and and the blue and .. line represents dynamics with initial value [0.2, 0.3, 0.5, 0.2]; the yellow and ... and the green and line represents dynamics with initial value [0.3, 0.4, 0.8, 0.6].

    Example 5. For the following set of parameters: Λ= 1.4, b1= 2.3, b2= 2.5, c1= 0.01, c2= 0.01, c3= 0.06, k1= 0.1, k2= 0.95, η1= 0.26, η2= 2.1, μ= 0.2, μ1= 0.8, μ2= 0.2, μ3= 0.5, θ1= 0.5, θ2= 0.2, δ= 0.01, h1= 0.1, h2= 0.1.

    In this case, simple calculation indicates that the sufficient condition (2) of Proposition 3.7 is satisfied. From Figure 5 we can see that the equilibrium E5=(0.4422, 0, 1.2195, 0.4048) is stable for different values of α and different sets of initial values: [X(0), Y1(0), Y2(0), Z(0)]=[0.2, 0.5, 0.3, 0.2], [0.3, 0.8, 0.4, 0.6].

    Figure 5.  (a)-(d) are the time series of the system (1.2), which show that the equilibrium E5 is stable for different values of α (α = 1, 0.8) when τ = 0; the red and and the blue and .. line represents dynamics with initial value [0.2, 0.5, 0.3, 0.2]; the yellow and ... and the green and line represents dynamics with initial value [0.3, 0.8, 0.4, 0.6].

    Remark 4.2. The above examples corresponding to the following ecological interpretation.

    (1) Figure 4 indicates that if the toxic phytoplankton is less competitive than the non-toxic phytoplankton, then the toxic phytoplankton will die out.

    (2) Figure 5 indicates that if the toxic phytoplankton win the competition between non-toxic phytoplankton and toxic phytoplankton, then the zooplankton may still survive under certain conditions, that is, nutrients, toxic phytoplankton and zooplankton may theoretically coexist. However, this phenomenon usually does not appear in real world.

    (3) From Figures 15 we can see that as the value of α decreases, the steady speed becomes slow for each equilibrium. This indicates that the value of α has obvious effects on the dynamical behaviors of the system.

    Example 6. For the following set of parameters: α= 0.8, Λ= 1.4, b1= 0.32, b2= 0.54, c1= 0.06, c2= 0.08, c3= 0.6, k1= 0.7, k2= 0.6, η1= 1.8, η2= 0.6, μ= 0.2, μ1= 0.4, μ2= 0.9, μ3= 0.5, θ2= 0.5, δ= 0.1, h1= 0.1, h2= 0.1.

    In this example, we will consider the influence of toxic, i.e., θ1. Here, we choose θ1=0,0.6, with the initial conditions: [X(0), Y1(0), Y2(0), Z(0)] = [0.2, 0.2, 0.3, 0.2], [0.8, 0.8, 0.6, 0.8].

    Example 7. In this example, we will consider the influence of α.

    (1) For the following set of parameters: Λ= 1.4, b1= 2.3, b2= 2.5, c1= 0.01, c2= 0.01, c3= 0.06, k1= 0.1, k2= 0.95, η1= 0.26, η2= 2.1, μ= 0.2, μ1= 0.8, μ2= 0.2, μ3= 0.5, θ1= 0.5, θ2= 0.2, δ= 0.01, h1= 0.1, h2= 0.1, with the initial conditions: [X(0), Y1(0), Y2(0), Z(0)] = [0.2, 0.5, 0.3, 0.2].

    (2) For the following set of parameters: Λ= 1.4, b1= 0.32, b2= 0.54, c1= 0.06, c2= 0.08, c3= 0.6, k1= 0.7, k2= 0.6, η1= 1.8, η2= 0.6, μ= 0.2, μ1= 0.4, μ2= 0.9, μ3= 0.5, θ1= 0.6, θ2= 0.5, δ= 0.1, h1= 0.1, h2= 0.1, with the initial conditions: [X(0), Y1(0), Y2(0), Z(0)] = [0.5, 0.5, 0.5, 0.2].

    Remark 4.3. (1) From Figure 7, we can see that if the value of α is relatively big(i.e., α = 1, 0.8), then the equilibrium E5 is locally stable; if the value of α is relatively small(i.e., α = 0.4), then the equilibrium is unstable, and oscillation may occur.

    Figure 7.  (a)-(d) are the time series of the system (1.2), which show the influence of α. (The initial conditions: [0.2, 0.5, 0.3, 0.2]).

    (2) From Figure 8, we can see that if the value of α is relatively big(i.e., α = 1, 0.7), then the equilibrium E6 is locally stable; if the value of α is relatively small(i.e., α = 0.1), then the equilibrium E0 is locally stable.

    Figure 8.  (a)-(d) are the time series of the system (1.2), which show the influence of α. (The initial conditions: [0.5, 0.5, 0.5, 0.2]).

    (3) Figure 7 and Figure 8 indicate that if the value of α is relatively small, the system will be destabilized.

    Example 8. For the following set of parameters: α= 0.89, Λ= 1.4, b1= 0.32, b2= 0.54, c1= 0.06, c2= 0.08, c3= 0.6, k1= 0.7, k2= 0.6, η1= 1.8, η2= 0.6, μ= 0.2, μ1= 0.4, μ2= 0.9, μ3= 0.5, θ1= 0.6, θ2= 0.5, δ= 0.1, h1= 0.1, h2= 0.1.

    In this example, the equilibrium E6=(3.2026, 0.4114, 0.2784, 0.1608), and our main aim is to study the effect of time delay on the stability of the system.

    (1) From Figure 9, we can see that if τ=0, then the equilibrium E6 is stable with different sets of initial value: [X(0), Y1(0), Y2(0), Z(0)]=[0.5, 0.5, 0.5, 0.2], [2.8, 1.8, 3.3, 1.3].

    Figure 9.  (a)-(d) are the time series of the system (1.2), which show that the equilibrium E6 is stable for α= 0.89 and τ=0, the blue and line represents dynamics with initial value [0.5, 0.5, 0.5, 0.2];the red and .. line represents dynamics with initial value [0.8, 0.8, 0.3, 0.3].

    (2) From Figure 10, we will see that if τ=3<τ04.4671, then the equilibrium E6 is stable with different sets of initial value: [X(0), Y1(0), Y2(0), Z(0)]=[0.5, 0.5, 0.5, 0.2], [2.8, 1.8, 3.3, 1.3].

    Figure 10.  (a)-(d) are the time series of the system (1.2), which show that the equilibrium E6 is stable for α= 0.89 and τ=3<τ04.4671, the blue and line represents dynamics with initial value [0.5, 0.5, 0.5, 0.2];the red and .. line represents dynamics with initial value [0.8, 0.8, 0.3, 0.3].

    (3) From Figure 11, we can see that if τ=4.2<τ04.4671, then the equilibrium E6 is stable with different sets of initial value: [X(0), Y1(0), Y2(0), Z(0)]=[0.5, 0.5, 0.5, 0.2], [2.8, 1.8, 3.3, 1.3].

    Figure 11.  (a)-(d) are the time series of the system (1.2), which show that the equilibrium E6 is stable for α= 0.89 and τ=4.2<τ04.4671, the blue and line represents dynamics with initial value [0.5, 0.5, 0.5, 0.2];the red and .. line represents dynamics with initial value [0.8, 0.8, 0.3, 0.3].

    (4) From Figure 12, we will see that if τ=5>τ04.4671, then periodic oscillation occurs, and the equilibrium E6 will lose its stability and periodic solutions appear through Hopf bifurcation, with different sets of initial value: [X(0), Y1(0), Y2(0), Z(0)]=[0.5, 0.5, 0.5, 0.2], [0.8, 0.8, 0.3, 0.3].

    Figure 12.  (a)-(d) are time series of the system (1.2), which show that the equilibrium E6 is unstable, and periodic oscillation occurs, for α= 0.89 and τ=5>τ04.4671. The blue and line represents the dynamics with initial value [0.5, 0.5, 0.5, 0.2]; while the red and .. line represents the dynamics with initial value [0.8, 0.8, 0.3, 0.3].

    (5) From Figure 13, we can see that if τ=6>τ04.4671, then similar periodic oscillation occurs as that in Figure 12, but with larger amplitude.

    Figure 13.  (a)-(d) are the time series of the system (1.2), which show that the equilibrium E6 is unstable, and periodic oscillation occurs, for α= 0.89 and τ=6>τ04.4671. The blue and line represents the dynamics with initial value [0.5, 0.5, 0.5, 0.2]; while the red and .. line represents the dynamics with initial value [0.8, 0.8, 0.3, 0.3].

    Remark 4.4. From the above example, we can see that the time delay has an effect of destabilizing the equilibrium E6. In other words, the larger the value of time delay is, the more possible that the equilibrium E6 lose its stability.

    In this paper, a fractional-order mathematical model is constructed to describe the active of nutrient-phytoplankton-toxic phytoplankton-zooplankton.

    Through qualitative analysis, we get the following results.

    We figure out the sufficient conditions for the existence and local stability of E0, E1, E2, E3, E4, E5, E6 for τ=0.

    By using time delay as a bifurcation parameter, the existence of Hopf bifurcation is analyzed in detail. We find that if τ<τ0, then the equilibrium E6 is locally stable; while it is unstable if τ>τ0 and Hopf bifurcation may occur near τ0.

    Through numerical simulation we get the following results.

    Figure 1 shows the stability of equilibrium E0 for different values of α.

    Figure 2 shows the stability of equilibrium E1 for different values of α.

    Figure 3 shows the stability of equilibrium E2 for different values of α.

    Figure 4 shows the stability of equilibrium E4 for different values of α.

    Figure 5 shows the stability of equilibrium E5 for different values of α.

    Figure 6 shows the effect of parameter θ1 on the system (1.2).

    Figure 6.  (a)-(d) are the time series of the system (1.2), which show the influence of θ1. (The initial conditions: [0.2, 0.2, 0.3, 0.2], [0.8, 0.8, 0.6, 0.8]).

    Figure 7 and Figure 8 indicate that the value of α is closely relate to the stability of each equilibrium. The stability of each equilibrium becomes weaker as the value of α decreases.

    Figures 911 show that E6 is stable if τ[0,τ0); Figures 1213 show that E6 is unstable if τ>τ0 and periodic oscillation may occur; Figures 913 indicate that the impact of τ on the dynamics of the system is crucial.

    Table 2 and Figure 14 show that the value of τ0 arises as the value of α increases.

    Table 2.  The effect of α on the values ω0, τ0 in system(1.2).
    Fractional order α Critical frequency ω0 Bifurcation point τ0
    0.60 0.0221 49.8601
    0.65 0.0318 30.9419
    0.70 0.0436 20.0341
    0.75 0.0576 13.3431
    0.80 0.0740 8.9972
    0.85 0.0927 6.1119
    0.90 0.1140 4.1175
    0.95 0.1378 2.7095
    1.00 0.1642 1.6926

     | Show Table
    DownLoad: CSV
    Figure 14.  Illustration of bifurcation τ0 versus fractional order α for system (1.2). The bifurcation points are becoming smaller and smaller as the value of α increase.

    Remark 5.1. When τ=0, Y1=Y2 and h1=h2=0, then system (1.2) will degenerated to the model in [40].

    In the model of this paper, phytoplankton is divided into two class, non-toxic phytoplankton and toxic phytoplankton. Some experimental data suggested that some zooplankton are capable of selecting for nontoxic phytoplankton, a mechanism that allows toxic phytoplankton to coexist with nontoxic phytoplankton [47]. This is also consistent with the results in Figure 9. Although some zooplankton have the ability to distinguish between toxic and non-toxic plants, some other experiments have shown that some zooplankton might not be able to distinguish between toxic and nontoxic algae, even shows a slight preference for toxic strains [48,49]. We can find from Figure 6 that once non-toxic phytoplankton become scarce, the zooplankton start to eat the toxic phytoplankton, even if the toxicity is weak, the zooplankton may become extinct. This is dangerous for the ecosystem.

    Remark 5.2. Any ecosystem depends on the natural environment, and in the real natural environment there are more or less physical and chemical signals that interact with the ecosystem. This motivate us to consider stochastic effects to the ecosystem, that is, white noise should be included into the system. We leave this as our next work.

    Each of the authors, Ruiqing Shi, Jianing Ren, Cuihong Wang contributed to each part of this work equally and read and approved the final version of the manuscript.

    This work is partly supported by National Natural Science Foundation of China (No. 61907027). The authors would like to thank the anonymous reviewers for their helpful comments, which improved the quality of this paper greatly.

    The authors declare that they have no financial or non-financial competing interests.



    [1] P. Vaupel, A. Mayer, Hypoxia in cancer: significance and impact on clinical outcome, Cancer Metastasis Rev., 26 (2007), 225–239. https://doi.org/10.1007/s10555-007-9055-1 doi: 10.1007/s10555-007-9055-1
    [2] C. T. Lee, M. K. Boss, M. W. Dewhirst, Imaging tumor hypoxia to advance radiation oncology, Antioxid. Redox Signaling, 21 (2014), 313–337. https://doi.org/10.1089/ars.2013.5759 doi: 10.1089/ars.2013.5759
    [3] M. R. Horsman, J. Overgaard, The impact of hypoxia and its modification of the outcome of radiotherapy, J. Radiat. Res., 57 (2016), i90–i98. https://doi.org/10.1093/jrr/rrw007 doi: 10.1093/jrr/rrw007
    [4] J. A. Forsythe, B. H. Jiang, N. V. Iyer, F. Agani, S. W. Leung, R. D. Koos, et al., Activation of vascular endothelial growth factor gene transcription by hypoxia-inducible factor 1, Mol. Cell. Biol., 16 (1996), 4604–4613. https://doi.org/10.1128/MCB.16.9.4604 doi: 10.1128/MCB.16.9.4604
    [5] A. Ahluwalia, A. S Tarnawski, Critical role of hypoxia sensor-hif-1α in vegf gene activation. implications for angiogenesis and tissue injury healing, Curr. Med. Chem., 19 (2012), 90–97. https://doi.org/10.2174/092986712803413944 doi: 10.2174/092986712803413944
    [6] E. K. Rofstad, Microenvironment-induced cancer metastasis, Int. J. Radiat. Biol., 76 (2000), 589–605. https://doi.org/10.1080/095530000138259 doi: 10.1080/095530000138259
    [7] M. G. Binker, A. A. Binker-Cosen, D. Richards, H. Y. Gaisano, R. H. de Cosen, L. I. Cosen-Binker, Hypoxia–reoxygenation increase invasiveness of panc-1 cells through rac1/mmp-2, Biochem. Biophys. Res. Commun., 393 (2010), 371–376. https://doi.org/10.1016/j.bbrc.2010.01.125 doi: 10.1016/j.bbrc.2010.01.125
    [8] R. M. Phillips, Targeting the hypoxic fraction of tumours using hypoxia-activated prodrugs, Cancer Chemother. Pharmacol., 77 (2016), 441–457. https://doi.org/10.1007/s00280-015-2920-7 doi: 10.1007/s00280-015-2920-7
    [9] W. A. Denny, The role of hypoxia-activated prodrugs in cancer therapy, Lancet Oncol., 1 (2000), 25–29. https://doi.org/10.1016/S1470-2045(00)00006-1 doi: 10.1016/S1470-2045(00)00006-1
    [10] S. G. Peeters, C. M. Zegers, R. Biemans, N. G. Lieuwes, R. G. van Stiphout, A. Yaromina, et al., TH-302 in combination with radiotherapy enhances the therapeutic outcome and is associated with pretreatment [18f] hx4 hypoxia pet imaging, Clin. Cancer Res., 21 (2015), 2984–2992. https://doi.org/10.1158/1078-0432.CCR-15-0018 doi: 10.1158/1078-0432.CCR-15-0018
    [11] V. Liapis, A. Labrinidis, I. Zinonos, S. Hay, V. Ponomarev, V. Panagopoulos, et al., Hypoxia-activated pro-drug th-302 exhibits potent tumor suppressive activity and cooperates with chemotherapy against osteosarcoma, Cancer Lett., 357 (2015), 160–169. https://doi.org/10.1016/j.canlet.2014.11.020 doi: 10.1016/j.canlet.2014.11.020
    [12] Q. Liu, J. D. Sun, J. Wang, D. Ahluwalia, A. F. Baker, L. D. Cranmer, et al., TH-302, a hypoxia-activated prodrug with broad in vivo preclinical combination therapy efficacy: optimization of dosing regimens and schedules, Cancer Chemother. Pharmacol., 69 (2012), 1487–1498. https://doi.org/10.1007/s00280-012-1852-8 doi: 10.1007/s00280-012-1852-8
    [13] Hypoxia/Normoxia, in Encyclopedic Reference of Genomics and Proteomics in Molecular Medicine, Springer, Berlin, Heidelberg, (2006), 853–853. https://doi.org/10.1007/3-540-29623-9_7440
    [14] C. Wigerup, S. Påhlman, D. Bexell, Therapeutic targeting of hypoxia and hypoxia-inducible factors in cancer, Pharmacol. Ther., 164 (2016), 152–169. https://doi.org/10.1016/j.pharmthera.2016.04.009 doi: 10.1016/j.pharmthera.2016.04.009
    [15] J. D. Sun, Q. Liu, J. Wang, D. Ahluwalia, D. Ferraro, Y. Wang, et al., Selective tumor hypoxia targeting by hypoxia-activated prodrug TH-302 inhibits tumor growth in preclinical models of cancer, Clin. Cancer Res., 18 (2012), 758–770. https://doi.org/10.1158/1078-0432.CCR-11-1980 doi: 10.1158/1078-0432.CCR-11-1980
    [16] F. Meng, J. W. Evans, D. Bhupathi, M. Banica, L. Lan, G. Lorente, et al., Molecular and cellular pharmacology of the hypoxia-activated prodrug TH-302, Mol. Cancer Ther., 11 (2012), 740–751. https://doi.org/10.1158/1535-7163.MCT-11-0634 doi: 10.1158/1535-7163.MCT-11-0634
    [17] Y. Huang, Y. Tian, Y. Zhao, C. Xue, J. Zhan, L. Liu, et al., Efficacy of the hypoxia-activated prodrug evofosfamide (TH-302) in nasopharyngeal carcinoma in vitro and in vivo, Cancer Commun., 38 (2018), 1–9. https://doi.org/10.1186/s40880-018-0285-0 doi: 10.1186/s40880-018-0285-0
    [18] G. J. Weiss, J. R. Infante, E. G. Chiorean, M. J. Borad, J. C. Bendell, J. R. Molina, et al., Phase 1 study of the safety, tolerability, and pharmacokinetics of TH-302, a hypoxia-activated prodrug, in patients with advanced solid malignancies, Clin. Cancer Res., 17 (2011), 2997–3004. https://doi.org/10.1158/1078-0432.CCR-10-3425 doi: 10.1158/1078-0432.CCR-10-3425
    [19] I. Lohse, J. Rasowski, P. Cao, M. Pintilie, T. Do, M. S. Tsao, et al., Targeting hypoxic microenvironment of pancreatic xenografts with the hypoxia-activated prodrug th-302, Oncotarget, 7 (2016), 33571. https://doi.org/10.18632/oncotarget.9654 doi: 10.18632/oncotarget.9654
    [20] S. Matsumoto, S. Kishimoto, K. Saito, Y. Takakusagi, J. P. Munasinghe, N. Devasahayam, et al., Metabolic and physiologic imaging biomarkers of the tumor microenvironment predict treatment outcome with radiation or a hypoxia-activated prodrug in mice, Cancer Res., 78 (2018), 3783–3792. https://doi.org/10.1158/0008-5472.CAN-18-0491 doi: 10.1158/0008-5472.CAN-18-0491
    [21] K. J. Nytko, I. Grgic, S. Bender, J. Ott, M. Guckenberger, O. Riesterer et al., The hypoxia-activated prodrug evofosfamide in combination with multiple regimens of radiotherapy, Oncotarget, 8 (2017), 23702. https://doi.org/10.18632/oncotarget.15784 doi: 10.18632/oncotarget.15784
    [22] J. K. Saggar, I. F. Tannock, Activity of the hypoxia-activated pro-drug TH-302 in hypoxic and perivascular regions of solid tumors and its potential to enhance therapeutic effects of chemotherapy, Int. J. Cancer, 134 (2014), 2726–2734. https://doi.org/10.1002/ijc.28595 doi: 10.1002/ijc.28595
    [23] V. Liapis, I. Zinonos, A. Labrinidis, S. Hay, V. Ponomarev, V. Panagopoulos, et al., Anticancer efficacy of the hypoxia-activated prodrug evofosfamide (th-302) in osteolytic breast cancer murine models, Cancer Med., 5 (2016), 534–545. https://doi.org/10.1002/cam4.599 doi: 10.1002/cam4.599
    [24] J. K. Saggar, I. F. Tannock, Chemotherapy rescues hypoxic tumor cells and induces their reoxygenation and repopulation–-an effect that is inhibited by the hypoxia-activated prodrug TH-302, Clin. Cancer Res., 21 (2015), 2107–2114. https://doi.org/10.1158/1078-0432.CCR-14-2298 doi: 10.1158/1078-0432.CCR-14-2298
    [25] J. D. Sun, Q. Liu, D. Ahluwalia, W. Li, F. Meng, Y. Wang, et al., Efficacy and safety of the hypoxia-activated prodrug TH-302 in combination with gemcitabine and nab-paclitaxel in human tumor xenograft models of pancreatic cancer, Cancer Biol. Ther., 16 (2015), 438–449. https://doi.org/10.1080/15384047.2014.1003005 doi: 10.1080/15384047.2014.1003005
    [26] S. P. Chawla, L. D. Cranmer, B. A. Van Tine, D. R. Reed, S. H. Okuno, J. E. Butrynski, et al., Phase ii study of the safety and antitumor activity of the hypoxia-activated prodrug TH-302 in combination with doxorubicin in patients with advanced soft tissue sarcoma, J. Clin. Oncol., 32 (2014), 3299. https://doi.org/10.1200/JCO.2013.54.3660 doi: 10.1200/JCO.2013.54.3660
    [27] K. N. Ganjoo, L. D. Cranmer, J. E. Butrynski, D. Rushing, D. Adkins, S. H. Okuno, et al., A phase i study of the safety and pharmacokinetics of the hypoxia-activated prodrug th-302 in combination with doxorubicin in patients with advanced soft tissue sarcoma, Oncology, 80 (2011), 50–56. https://doi.org/10.1159/000327739 doi: 10.1159/000327739
    [28] J. Von Pawel, R. von Roemeling, U. Gatzemeier, M. Boyer, L. O. Elisson, P. Clark, et al., Tirapazamine plus cisplatin versus cisplatin in advanced non-small-cell lung cancer: A report of the international catapult I study group, J. Clin. Oncol., 18 (2000), 1351–1359. https://doi.org/10.1200/JCO.2000.18.6.1351 doi: 10.1200/JCO.2000.18.6.1351
    [29] F. W. Hunter, B. G. Wouters, W. R. Wilson, Hypoxia-activated prodrugs: paths forward in the era of personalised medicine, Br. J. Cancer, 114 (2016), 1071–1077. https://doi.org/10.1038/bjc.2016.79 doi: 10.1038/bjc.2016.79
    [30] J. Sun, Q. Liu, D. Ahluwalia, J. Curd, M. Matteucci, C. Hart, Complementary chemotherapies with th-302, a novel hypoxia activated prodrug: optimization of dosing regimens and schedules for study in phase 1/2 with docetaxel, gemcitabine, pemetrexed, and doxorubicin, Biosymposia: Hypoxia, Ischemia, Inflammation, 2008.
    [31] S. B. Reddy, S. K. Williamson, Tirapazamine: a novel agent targeting hypoxic tumor cells, Expert Opin. Invest. Drugs, 18 (2009), 77–87. https://doi.org/10.1517/13543780802567250 doi: 10.1517/13543780802567250
    [32] J. C. Forster, L. G. Marcu, E. Bezak, Approaches to combat hypoxia in cancer therapy and the potential for in silico models in their evaluation, Physica Med., 64 (2019), 145–156. https://doi.org/10.1016/j.ejmp.2019.07.006 doi: 10.1016/j.ejmp.2019.07.006
    [33] W. Tuckwell, E. Bezak, E. Yeoh, L. Marcu, Efficient Monte Carlo modelling of individual tumour cell propagation for hypoxic head and neck cancer, Phys. Med. Biol., 53 (2008), 4489. https://doi.org/10.1088/0031-9155/53/17/002 doi: 10.1088/0031-9155/53/17/002
    [34] W. M. Harriss-Phillips, E. Bezak, E. Yeoh, Monte Carlo radiotherapy simulations of accelerated repopulation and reoxygenation for hypoxic head and neck cancer, Br. J. Radiol., 84 (2011), 903–918. https://doi.org/10.1259/bjr/25012212 doi: 10.1259/bjr/25012212
    [35] W. M. Harriss-Phillips, E. Bezak, A. Potter, Stochastic predictions of cell kill during stereotactic ablative radiation therapy: Do hypoxia and reoxygenation really matter?, Int. J. Radiat. Oncol. Biol. Phys., 95 (2016), 1290–1297. https://doi.org/10.1016/j.ijrobp.2016.03.014 doi: 10.1016/j.ijrobp.2016.03.014
    [36] L. G. Marcu, D. Marcu, S. M. Filip, In silico study of the impact of cancer stem cell dynamics and radiobiological hypoxia on tumour response to hyperfractionated radiotherapy, Cell proliferation, 49 (2016), 304–314.
    [37] E. Lindblom, I. Toma-Dasu, A. Dasu, Accounting for two forms of hypoxia for predicting tumour control probability in radiotherapy: An in silico study, in Oxygen Transport to Tissue XL. Advances in Experimental Medicine and Biology, 1072 (2018), 183–187. https://doi.org/10.1007/978-3-319-91287-5_29
    [38] A. Foehrenbacher, K. Patel, M. Abbattista, C. P. Guise, T. W. Secomb, W. R. Wilson, et al., The role of bystander effects in the antitumor activity of the hypoxia-activated prodrug pr-104, Frontiers in oncology, 3 (2013), 263. https://doi.org/10.3389/fonc.2013.00263 doi: 10.3389/fonc.2013.00263
    [39] A. B. Foehrenbacher, T. W. Secomb, W. R. Wilson, K. O. Hicks, Design of optimized hypoxia-activated prodrugs using pharmacokinetic/pharmacodynamic modeling, Front. Oncol., 3 (2013), 314. https://doi.org/10.3389/fonc.2013.00314 doi: 10.3389/fonc.2013.00314
    [40] K. O. Hicks, F. B. Pruijn, T. W. Secomb, M. P. Hay, R. Hsu, J. M. Brown, et al., Use of three-dimensional tissue cultures to model extravascular transport and predict in vivo activity of hypoxia-targeted anticancer drugs, J. Nat. Cancer Inst., 98 (2006), 1118–1128. https://doi.org/10.1093/jnci/djj306 doi: 10.1093/jnci/djj306
    [41] C. Meaney, G. G. Powathil, A. Yaromina, L. J. Dubois, P. Lambin, M. Kohandel, Role of hypoxia-activated prodrugs in combination with radiation therapy: An in silico approach, Math. Biosci. Eng., 16 (2019), 6257. https://doi.org/10.3934/mbe.2019312 doi: 10.3934/mbe.2019312
    [42] C. Meaney, S. Rhebergen, M. Kohandel, In silico analysis of hypoxia activated prodrugs in combination with anti angiogenic therapy through nanocell delivery, PLoS Comput. Biol., 16 (2020), e1007926. https://doi.org/10.1371/journal.pcbi.1007926 doi: 10.1371/journal.pcbi.1007926
    [43] S. Hamis, M. Kohandel, L. J. Dubois, A. Yaromina, P. Lambin, G. G. Powathil, Combining hypoxia-activated prodrugs and radiotherapy in silico: Impact of treatment scheduling and the intra-tumoural oxygen landscape, PLoS Comput. Biol., 16 (2020), e1008041. https://doi.org/10.1371/journal.pcbi.1008041 doi: 10.1371/journal.pcbi.1008041
    [44] X. Mao, S. McManaway, J. K. Jaiswal, P. B. Patel, W. R. Wilson, K. O. Hicks, et al., An agent-based model for drug-radiation interactions in the tumour microenvironment: Hypoxia-activated prodrug sn30000 in multicellular tumour spheroids, PLoS Comput. Biol., 14 (2018), e1006469. https://doi.org/10.1371/journal.pcbi.1006469 doi: 10.1371/journal.pcbi.1006469
    [45] C. R. Hong, S. Y. Mehta, H. Liyanage, S. P. McManaway, H. H. Lee, J. K. Jaiswal, et al., Spatially-resolved pharmacokinetic/pharmacodynamic modelling of bystander effects of a nitrochloromethylbenzindoline hypoxia-activated prodrug, Cancer Chemother. Pharmacol., 88 (2021), 673–687. https://doi.org/10.1007/s00280-021-04320-3 doi: 10.1007/s00280-021-04320-3
    [46] K. O. Hicks, B. G. Siim, J. K. Jaiswal, F. B. Pruijn, A. M. Fraser, R. Patel, et al., Pharmacokinetic/pharmacodynamic modeling identifies sn30000 and sn29751 as tirapazamine analogues with improved tissue penetration and hypoxic cell killing in tumors, Clin. Cancer Res., 16 (2010), 4946–4957. https://doi.org/10.1158/1078-0432.CCR-10-1439 doi: 10.1158/1078-0432.CCR-10-1439
    [47] C. R. Hong, W. R. Wilson, K. O. Hicks, An intratumor pharmacokinetic/pharmacodynamic model for the hypoxia-activated prodrug evofosfamide (TH-302): Monotherapy activity is not dependent on a bystander effect, Neoplasia, 21 (2019), 159–171. https://doi.org/10.1016/j.neo.2018.11.009 doi: 10.1016/j.neo.2018.11.009
    [48] C. R. Hong, G. Bogle, J. Wang, K. Patel, F. B. Pruijn, W. R. Wilson, et al., Bystander effects of hypoxia-activated prodrugs: agent-based modeling using three dimensional cell cultures, Front. Pharmacol., 9 (2018), 1013. https://doi.org/10.3389/fphar.2018.01013 doi: 10.3389/fphar.2018.01013
    [49] D. Lindsay, C. M. Garvey, S. M. Mumenthaler, J. Foo, Leveraging hypoxia-activated prodrugs to prevent drug resistance in solid tumors, PLoS Comput. Biol., 12 (2016), e1005077. https://doi.org/10.1371/journal.pcbi.1005077 doi: 10.1371/journal.pcbi.1005077
    [50] S. Yonucu, D. Yilmaz, C. Phipps, M. B. Unlu, M. Kohandel, Quantifying the effects of antiangiogenic and chemotherapy drug combinations on drug delivery and treatment efficacy, PLoS Comput. Biol., 13 (2017), 1–17, https://doi.org/10.1371/journal.pcbi.1005724 doi: 10.1371/journal.pcbi.1005724
    [51] M. Kohandel, M. Kardar, M. Milosevic, S. Sivaloganathan, Dynamics of tumor growth and combination of anti-angiogenic and cytotoxic therapies, Phys. Med. Biol., 52 (2007), 3665. https://doi.org/10.1088/0031-9155/52/13/001 doi: 10.1088/0031-9155/52/13/001
    [52] G. Powathil, M. Kohandel, S. Sivaloganathan, A. Oza, M. Milosevic, Mathematical modeling of brain tumors: effects of radiotherapy and chemotherapy, Phys. Med. Biol., 52 (2007), 3291. https://doi.org/10.1088/0031-9155/52/11/023 doi: 10.1088/0031-9155/52/11/023
    [53] E. Dalah, D. Bradley, A. Nisbet, A mathematical approach towards simulating a realistic tissue activity curve of 64Cu-ATSM for the purpose of sub-target volume delineation in radiotherapy, Nucl. Instrum. Methods Phys. Res., Sect. A, 619 (2010), 283–286. https://doi.org/10.1016/j.nima.2009.10.160 doi: 10.1016/j.nima.2009.10.160
    [54] E. Dalah, D. Bradley, A. Nisbet, Simulation of tissue activity curves of 64Cu-ATSM for sub-target volume delineation in radiotherapy, Phys. Med. Biol., 55 (2010), 681.
    [55] F. Winkler, S. V. Kozin, R. T. Tong, S. S. Chae, M. F. Booth, I. Garkavtsev, et al., Kinetics of vascular normalization by vegfr2 blockade governs brain tumor response to radiation: role of oxygenation, angiopoietin-1, and matrix metalloproteinases, Cancer Cell, 6 (2004), 553–563. https://doi.org/10.1016/j.ccr.2004.10.011 doi: 10.1016/j.ccr.2004.10.011
    [56] W. R. Wilson, M. P. Hay, Targeting hypoxia in cancer therapy, Nat. Rev. Cancer, 11 (2011), 393–410. https://doi.org/10.1038/nrc3064 doi: 10.1038/nrc3064
    [57] L. J. Nugent, R. K. Jain, Extravascular diffusion in normal and neoplastic tissues, Cancer Res., 44 (1984), 238–244.
    [58] R. K. Jain, Transport of molecules in the tumor interstitium: a review, Cancer Res., 47 (1987), 3039–3051.
    [59] L. S. Goodman, Goodman and Gilman's the Pharmacological Basis of Therapeutics, McGraw-Hill New York, 1549 (1996).
    [60] K. M. Laginha, S. Verwoert, G. J. Charrois, T. M. Allen, Determination of doxorubicin levels in whole tumor and tumor nuclei in murine breast cancer tumors, Clin. Cancer Res., 11 (2005), 6944–6949. https://doi.org/10.1158/1078-0432.CCR-05-0343 doi: 10.1158/1078-0432.CCR-05-0343
    [61] A. W. El-Kareh, T. W. Secomb, A mathematical model for comparison of bolus injection, continuous infusion, and liposomal delivery of doxorubicin to tumor cells, Neoplasia, 2 (2000), 325.
    [62] A. J. Leu, D. A. Berk, A. Lymboussaki, K. Alitalo, R. K. Jain, Absence of functional lymphatics within a murine sarcoma: a molecular and functional evaluation, Cancer Res., 60 (2000), 4324–4327.
    [63] M. Wu, H. B. Frieboes, S. R. McDougall, M. A. Chaplain, V. Cristini, J. Lowengrub, The effect of interstitial pressure on tumor growth: coupling with the blood and lymphatic vascular systems, J. Theor. Biol., 320 (2013), 131–151. https://doi.org/10.1016/j.jtbi.2012.11.031 doi: 10.1016/j.jtbi.2012.11.031
    [64] G. Powathil, M. Kohandel, M. Milosevic, S. Sivaloganathan, Modeling the spatial distribution of chronic tumor hypoxia: implications for experimental and clinical studies, Comput. Math. Methods Med., 2012 (2012). https: //doi.org/10.1155/2012/410602
    [65] A. S. Fung, J. Jonkman, I. F. Tannock, Quantitative immunohistochemistry for evaluating the distribution of Ki67 and other biomarkers in tumor sections and use of the method to study repopulation in xenografts after treatment with paclitaxel, Neoplasia, 14 (2012), 324–IN6. https://doi.org/10.1593/neo.12346 doi: 10.1593/neo.12346
    [66] X. Mao, S. McManaway, J. K. Jaiswal, C. R. Hong, W. R. Wilson, K. O. Hicks, Schedule-dependent potentiation of chemotherapy drugs by the hypoxia-activated prodrug sn30000, Cancer Biol. Ther., 20 (2019), 1258–1269. https://doi.org/10.1080/15384047.2019.1617570 doi: 10.1080/15384047.2019.1617570
    [67] S. Harris, P. Mistry, C. Freathy, J. Brown, P. Charlton, Antitumour activity of xr5944 in vitro and in vivo in combination with 5-fluorouracil and irinotecan in colon cancer cell lines, Br. J. Cancer, 92 (2005), 722–728. https://doi.org/10.1038/sj.bjc.6602403 doi: 10.1038/sj.bjc.6602403
    [68] J. M. Saucier, J. Yu, A. Gaikwad, R. L. Coleman, J. K. Wolf, J. A. Smith, Determination of the optimal combination chemotherapy regimen for treatment of platinum-resistant ovarian cancer in nude mouse model, J. Oncol. Pharm. Pract., 13 (2007), 39–45. https://doi.org/10.1177/1078155207077948 doi: 10.1177/1078155207077948
    [69] M. R. Horsman, L. S. Mortensen, J. B. Petersen, M. Busk, J. Overgaard, Imaging hypoxia to improve radiotherapy outcome, Nat. Rev. Clin. Oncol., 9 (2012), 674–687. https://doi.org/10.1038/nrclinonc.2012.171 doi: 10.1038/nrclinonc.2012.171
    [70] J. Overgaard, Hypoxic modification of radiotherapy in squamous cell carcinoma of the head and neck–a systematic review and meta-analysis, Radiother. Oncol., 100 (2011), 22–32.
    [71] M. Kovacs, D. Hocking, J. Evans, B. Siim, B. Wouters, J. Brown, Cisplatin anti-tumour potentiation by tirapazamine results from a hypoxia-dependent cellular sensitization to cisplatin, Br. J. Cancer, 80 (1999), 1245–1251. https://doi.org/10.1038/sj.bjc.6690492 doi: 10.1038/sj.bjc.6690492
    [72] X. Zhang, J. W. Wojtkowiak, G. V. Martinez, H. H. Cornnell, C. P. Hart, A. F. Baker, et al., Mr imaging biomarkers to monitor early response to hypoxia-activated prodrug TH-302 in pancreatic cancer xenografts, PloS One, 11 (2016), e0155289. https://doi.org/10.1371/journal.pone.0155289 doi: 10.1371/journal.pone.0155289
    [73] J. Cárdenas-Rodríguez, Y. Li, J. P. Galons, H. Cornnell, R. J. Gillies, M. D. Pagel, et al., Imaging biomarkers to monitor response to the hypoxia-activated prodrug TH-302 in the miapaca2 flank xenograft model, Magn. Reson. Imaging, 30 (2012), 1002–1009. https://doi.org/10.1016/j.mri.2012.02.015 doi: 10.1016/j.mri.2012.02.015
    [74] A. Leimgruber, K. Hickson, S. T. Lee, H. K. Gan, L. M. Cher, J. I. Sachinidis, et al., Spatial and quantitative mapping of glycolysis and hypoxia in glioblastoma as a predictor of radiotherapy response and sites of relapse, Eur. J. Nucl. Med. Mol. Imaging, 47 (2020), 1476–1485. https://doi.org/10.1007/s00259-020-04706-0 doi: 10.1007/s00259-020-04706-0
    [75] B. V. Jardim-Perassi, W. Mu, S. Huang, M. R. Tomaszewski, J. Poleszczuk, M. A. Abdalah, et al., Deep-learning and mr images to target hypoxic habitats with evofosfamide in preclinical models of sarcoma, Theranostics, 11 (2021), 5313. https://doi.org/10.7150/thno.56595 doi: 10.7150/thno.56595
    [76] J. C. Forster, M. J. Douglass, W. M. Harriss-Phillips, E. Bezak, Development of an in silico stochastic 4D model of tumor growth with angiogenesis, Med. Phys., 44 (2017), 1563–1576.
    [77] W. M. Harriss-Phillips, E. Bezak, E. Yeoh, The HYP-RT hypoxic tumour radiotherapy algorithm and accelerated repopulation dose per fraction study, Comput. Math. Methods Med., 2012 (2012). https://doi.org/10.1155/2012/363564
  • This article has been cited by:

    1. Shuang Chen, Yuanjin Ren, Small amplitude periodic solution of Hopf Bifurcation Theorem for fractional differential equations of balance point in group competitive martial arts, 2022, 7, 2444-8656, 207, 10.2478/amns.2021.2.00152
    2. Yuanyuan Ma, Nan Dong, Na Liu, Leilei Xie, Hopf bifurcation of the model with terms of two time-delays and delay-dependent parameter based on the theory of crossing curves, 2022, 32, 1054-1500, 103114, 10.1063/5.0095794
    3. Songkran Pleumpreedaporn, Chanidaporn Pleumpreedaporn, Jutarat Kongson, Chatthai Thaiprayoon, Jehad Alzabut, Weerawat Sudsutad, Dynamical Analysis of Nutrient-Phytoplankton-Zooplankton Model with Viral Disease in Phytoplankton Species under Atangana-Baleanu-Caputo Derivative, 2022, 10, 2227-7390, 1578, 10.3390/math10091578
    4. Muhammad Farman, Aamir Shehzad, Kottakkaran Sooppy Nisar, Evren Hincal, Ali Akgul, Ahmed Muhammad Hassan, Generalized Ulam-Hyers-Rassias stability and novel sustainable techniques for dynamical analysis of global warming impact on ecosystem, 2023, 13, 2045-2322, 10.1038/s41598-023-49806-7
    5. Kaushik Dehingia, Salah Boulaaras, Suman Gogoi, On the dynamics of a nutrient–plankton system with Caputo and Caputo–Fabrizio fractional operators, 2024, 76, 18777503, 102232, 10.1016/j.jocs.2024.102232
    6. R.N. Premakumari, Chandrali Baishya, Mohammad Esmael Samei, Manisha Krishna Naik, A novel optimal control strategy for nutrient–phytoplankton–zooplankton model with viral infection in plankton, 2024, 137, 10075704, 108157, 10.1016/j.cnsns.2024.108157
    7. Yumei Lin, Yuan Ma, Yunxian Dai, Stability and Hopf Bifurcation Control for Fractional-Order Two-Gene Regulatory Network With Multiple Delays, 2023, 11, 2169-3536, 58389, 10.1109/ACCESS.2023.3283401
  • Reader Comments
  • © 2022 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(2289) PDF downloads(104) Cited by(0)

Figures and Tables

Figures(5)  /  Tables(3)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog