Dynamics of a stoichiometric phytoplankton-zooplankton model with season-driven light intensity

Zhenyao Sun; Da Song; Meng Fan; Zhenyao Sun; Da Song; Meng Fan

doi:10.3934/mbe.2024301

Mathematical Biosciences and Engineering

2024, Volume 21, Issue 8: 6870-6897. doi: 10.3934/mbe.2024301

Previous Article Next Article

Research article Special Issues

Dynamics of a stoichiometric phytoplankton-zooplankton model with season-driven light intensity

1.
School of Mathematics and Statistics, Northeast Normal University, 5268 Renmin Street, Changchun, Jilin, 130024, China
2.
School of Science, Shanghai Maritime University, 1550 Haigang Avenue, Shanghai, 201306, China

Received: 05 March 2024 Revised: 10 July 2024 Accepted: 06 August 2024 Published: 20 August 2024

Chemical heterogeneity significantly influences the dynamics of phytoplankton and zooplankton interactions through its effects on phytoplankton carrying capacity and zooplankton ingestion rates. Our central objective of this study was to develop and examine a nonautonomous model of phytoplankton-zooplankton growth, which incorporates season-driven variations in light intensity and chemical heterogeneity. The dynamics of the system is characterized by positive invariance, dissipativity, boundary dynamics, and internal dynamics. Subsequently, numerical simulations were conducted to validate the theoretical findings and to elucidate the effects of seasonal light intensity, nutrient availability, and zooplankton loss rates on phytoplankton dynamics. The outcomes of our model and analysis offer a potential explanation for seasonal phytoplankton blooms.

Keywords:

Citation: Zhenyao Sun, Da Song, Meng Fan. Dynamics of a stoichiometric phytoplankton-zooplankton model with season-driven light intensity[J]. Mathematical Biosciences and Engineering, 2024, 21(8): 6870-6897. doi: 10.3934/mbe.2024301

Related Papers:

[1]	Shuguang Li, Yong Sun, Muhammad Ijaz Khan . Single machine Pareto scheduling with positional deadlines and agreeable release and processing times. Electronic Research Archive, 2023, 31(5): 3050-3063. doi: 10.3934/era.2023154
[2]	Chuanyao Li, Yichao Lu, Yuqiang Wang, Gege Jiang . Congestion behavior and tolling strategies in a bottleneck model with exponential scheduling preference. Electronic Research Archive, 2023, 31(2): 1065-1088. doi: 10.3934/era.2023053
[3]	Xiaopeng Yi, Huey Tyng Cheong, Zhaohua Gong, Chongyang Liu, Kok Lay Teo . Dynamic optimization of a two-stage fractional system in microbial batch process. Electronic Research Archive, 2024, 32(12): 6680-6697. doi: 10.3934/era.2024312
[4]	E. A. Abdel-Rehim . The time evolution of the large exponential and power population growth and their relation to the discrete linear birth-death process. Electronic Research Archive, 2022, 30(7): 2487-2509. doi: 10.3934/era.2022127
[5]	Pan Zhang, Lan Huang . Stability for a 3D Ladyzhenskaya fluid model with unbounded variable delay. Electronic Research Archive, 2023, 31(12): 7602-7627. doi: 10.3934/era.2023384
[6]	Ming Wei, Congxin Yang, Bo Sun, Binbin Jing . A multi-objective optimization model for green demand responsive airport shuttle scheduling with a stop location problem. Electronic Research Archive, 2023, 31(10): 6363-6383. doi: 10.3934/era.2023322
[7]	Jun Liu, Yue Liu, Xiaoge Yu, Xiao Ye . An efficient numerical method based on QSC for multi-term variable-order time fractional mobile-immobile diffusion equation with Neumann boundary condition. Electronic Research Archive, 2025, 33(2): 642-666. doi: 10.3934/era.2025030
[8]	Yu Shen, Hecheng Li . A multi-strategy genetic algorithm for solving multi-point dynamic aggregation problems with priority relationships of tasks. Electronic Research Archive, 2024, 32(1): 445-472. doi: 10.3934/era.2024022
[9]	Li-Bin Liu, Ying Liang, Jian Zhang, Xiaobing Bao . A robust adaptive grid method for singularly perturbed Burger-Huxley equations. Electronic Research Archive, 2020, 28(4): 1439-1457. doi: 10.3934/era.2020076
[10]	Jun Pan, Yuelong Tang . Two-grid $H^1$ -Galerkin mixed finite elements combined with $L1$ scheme for nonlinear time fractional parabolic equations. Electronic Research Archive, 2023, 31(12): 7207-7223. doi: 10.3934/era.2023365

Abstract

1. Introduction

Tumourigenesis is a complex, multi-step process which leads to the formation of solid malignancies. A critical step in this process is angiogenesis in which tumours develop their own vasculature and blood supply. In normal tissues, angiogenesis is carefully self-regulated and tightly controlled, whereas in tumours, the vessels are structurally and functionally abnormal. Accordingly, they are often characterized by defective endothelia, basement membranes, and pericyte coverage leading to inefficient nutrient delivery to tumour cells despite a high global blood flow in the environment ^[1]. In particular, the poor delivery of oxygen results in regions of severe hypoxia—a trait observed in nearly all solid malignancies ^[2,3]. Tumour hypoxia has been linked to the increase of many cancerous behaviors such as genomic instability, malignant progression, and metastasis formation ^[4,5].

The clinical implications of hypoxia for cancer therapy are also of great importance, most evidently in the context of radiation and chemotherapy. It has been established that resistance to radiation is conferred by hypoxia, and the mechanism by which this is achieved is well understood ^[4,6]. The radioresistance arises in the hypoxic regions because of the relatively fewer oxygen-derived free radicals generated by ionizing radiation that compromise the cytotoxic effects of radiotherapy. In the case of chemotherapy, hypoxia has also been shown to elicit an overall decrease in efficacy in vivo. This is likely related to inefficient drug delivery due to the abnormal tumour vasculature, though the underlying mechanisms dictating this are less well understood ^[4,7].

On the other hand, hypoxia is a theoretically ideal target for guided drug delivery. To overcome hypoxia-induced treatment resistance, hypoxia-activated prodrugs (HAPs) have been developed that are nontoxic under physiological oxygen concentrations but are activated through a bioreduction specifically in hypoxic tumour regions. Once HAPs become activated at hypoxic sites, some of them are able to diffuse back out of hypoxic areas allowing them to attack nearby, non-hypoxic cells as well—a phenomenon called the bystander effect. These bystander effects play a vital role in the overall cytotoxic effectiveness of HAPs. Furthermore, the failure of several recent clinical studies involving HAPs suggests that one of the several factors contributing to this failure is insufficient classification of patients based on their hypoxia and nitroreductase expression statuses ^[8,9]. Hence, not only the presence of hypoxia, but also its spatial distribution, might significantly influence the efficacy of treatments involving HAPs. Preclinical studies have examined the effect of HAPs, administered singly or in combination with conventional treatments such as chemotherapy and radiotherapy, resulting in very promising outcomes in the control of tumour growth ^[10,11,12]^*. Despite the theoretical and preclinical results, several HAPs have recently failed to demonstrate efficacy in clinical trials. The disconnect between theory and practice is hypothesized to be mainly due to sub-optimal patient selection and heterogeneities in distributions of hypoxia, in addition to other issues common to chemotherapies such as poor delivery due to the insufficient blood supply and high interstitial fluid pressures ^[9,13]. A better understanding of these therapies and their interactions with tumour cells and the surrounding micro-environment for individual patients is required before they can be successfully translated into clinical practice.

^*https://www.youtube.com/watch?v=LjbfbQ_OON8&t=s

Mathematical modelling approaches have been used to study the interaction between hypoxia and external beam radiation therapy ^{[14,15,16,17]} as well as the action of HAPs ^{[18,19,20,21,22]}. Here, we developed an experimentally motivated mathematical model to predict the effect of HAPs by incorporating the amount of tumour hypoxia and the action of the prodrug as well as the sensitivity of activated metabolites. Once the model is calibrated and validated with appropriate experimental data, it can be used to determine the efficient administration of HAPs for a given tumour and to predict the outcome of a combination of HAPs and radiotherapy. Moreover, the developed mathematical model can be further used to analyse and generate testable hypotheses that can be further studied experimentally.

2. Mathematical model

Our mathematical model is developed based on the experimental observations of TH-302 administration by Peeters et al. ^[11]. In this study, Rhabdomyosarcoma R1 and H460 NSCLC xenografts were treated with TH-302 and external beam radiation therapy, and the effects on tumour size and hypoxic fraction were examined. Furthermore, the TH-302 efficacy based on initial hypoxia level was assessed as well as efficacy of its combination with radiation as well as order of treatment administration. In order to study the effect of HAPs and their bystander effects on tumour cells— particularly hypoxic tumour cells—we consider the spatio-temporal evolution of tumour cell concentration, oxygen distribution, inactive HAP distribution and active metabolite distribution. If we denote $n({\bf x}, t), \; k({\bf x}, t), \; c({\bf x}, t)$ and $c_a ({\bf x}, t)$ as the concentrations of tumour cells, oxygen, HAP (inactive) and active metabolite at the location ${\bf x}$ and time $t$ , we can write the model equations as given below.

Cancer cells

$\begin{eqnarray} &&\frac{\partial n({\bf x}, t)}{\partial t} ~ = ~ \underbrace{D_n\nabla^{2}n({\bf x}, t)}_{\text{Diffusion}} + \underbrace{\frac{r\, k({\bf x}, t)^p}{k_g^p+k({\bf x}, t)^p} \, n({\bf x}, t)}_{\begin{smallmatrix} \text{Oxygen-dependent} \\ \text{growth} \end{smallmatrix}}- \underbrace{\delta_{a} c_a({\bf x}, t) \, n({\bf x}, t)}_{\begin{smallmatrix} \text{Cell death by} \\ \text{activated HAPs}\end{smallmatrix}} - \underbrace{\alpha_k \, d_t \, n({\bf x}, t)}_{\begin{smallmatrix} \text{Cell death by} \\ \text{radiaion}\end{smallmatrix}}\end{eqnarray}$

(1)

$\begin{eqnarray} &&\alpha_k = \alpha_h \frac{(\alpha_n/\alpha_h) \, k({\bf x}, t) + k_s}{k({\bf x}, t)+k_s}~~, ~~d_t = \frac{d}{\tau}f(t)\end{eqnarray}$

(2)

Oxygen

$\begin{eqnarray} \frac{\partial k({\bf x}, t)}{\partial t} ~ = ~ 0 ~ = ~ \underbrace{D_k\nabla^{2} k({\bf x}, t)}_{\text{Diffusion}} + \underbrace{r_k v({\bf x}, t)}_{\text{Production}} - \underbrace{\frac{q_k k({\bf x}, t)^m}{k_c^m+ k({\bf x}, t)^m} \, n({\bf x}, t)}_{\begin{smallmatrix} \text{Oxygen-dependent} \\ \text{consumption} \end{smallmatrix}}\end{eqnarray}$

(3)

Hypoxia activated prodrugs (HAPs)

$\begin{eqnarray} \frac{\partial c({\bf x}, t)}{\partial t} ~ = ~ 0 ~ = ~ \underbrace{D_c\nabla^{2} c({\bf x}, t)}_{\text{Diffusion}} + \underbrace{r_c g(t) v({\bf x}, t)}_{\text{Production}} - \underbrace{\frac{q_a k_a^u}{k_a^u+k({\bf x}, t)^u} \, c({\bf x}, t)}_{\text{Activation}} - \underbrace{\lambda_c c({\bf x}, t)}_{\text{Decay}}\end{eqnarray}$

(4)

Activated HAP (AHAPs)

$\begin{eqnarray} \frac{\partial c_a({\bf x}, t)}{\partial t} \; = \; 0 \; = \; \underbrace{D_{a}\nabla^{2} c_a({\bf x}, t)}_{\text{Diffusion}} + \underbrace{\frac{q_a k_a^u}{k_a^u+k({\bf x}, t)^u} \, c({\bf x}, t)}_{\text{Activation}} - \underbrace{\lambda_a c_a({\bf x}, t)}_{\text{Decay}} = 0 \end{eqnarray}$

(5)

The description of the relevant parameters and corresponding values are given in . Here, we assume that cancer cells ( $n({\bf x}, t)$ ) follow an oxygen-dependent, exponential growth law and are treated with HAPs alone or in combination with external beam radiation therapy. The oxygen-dependent growth is incorporated to account for the well-known phenomenon of cancer cells altering their growth pathway in the absence of oxygen ^[23]. Mathematically, this is incorporated by the use of a Hill function, which is commonly used to model the dependence of a rate on a specific quantity. Notice that as the oxygen concentration increases, the Hill function in the oxygen-dependent growth term approaches $r$ (its maximum growth rate) and that as oxygen decreases, the growth rate approaches zero. The parameter $\delta_a$ , the rate of HAP induced cell death, is assumed to be cell-specific and reflects the death rate of the cancer cells to the activated HAPs (sensitivity to warhead). The radiation induced cell kill is incorporated into the model using a linear term with a radiosensitivity parameter $\alpha_k$ and dose rate $d_t$ . While $\alpha_k$ is a constant dependent on the tissue of interest, we employ an oxygen enhancement ratio to show the decreased effect of radiation in hypoxic areas of tumours. This is common in modelling the interaction of radiation with hypoxia and is shown above in Eq (2) where oxygen concentration is incorporated into the equation for $\alpha_k$ ^[19]. When modelling external beam radiation therapy, the linear-quadratic (LQ) model is often used to describe tumour survival. Here, we make the simplification that the linear term is dominant and exclusion of the quadratic term does not impact the qualitative results (as done in ^[14] for example). Radiation schedule is incorporated with the model using the function, $f(t)$ , where $f(t) = 1$ , when the radiation is given.

Table 1. List of model parameters (

$^*$ for assumed and estimated parameters, refer to parameters section).

Parameter	Symbol	Value (Unit)	Reference
Diffusion coefficient (tumour cells)	$D_n$	0.32 (mm $^2$ /day)	^[14]
Growth rate	$r$	0.35 (1/day)	^[14]
Hill coefficient	$p$	8	Assumed $^*$
pO $_2$ at which the growth is half the maximum rate	$k_g$	10 (mmHg)	Assumed $^*$
Rate of HAP-induced cell death	$\delta_a$	10 (1/day[ $C_a$ ])	Assumed $^*$
Normoxic radiosensitivity	$\alpha_n$	0.3 (1/Gy)	^[16]
Hypoxic radiosensitivity	$\alpha_h$	0.1 (1/Gy)	^[16]
Constant	$k_s$	3 (mmHg)	^[16]
Radiation dose	$d$	2-8 (Gy)	^[11]
Irradiation duration	$\tau$	3 (minutes)	^[16]
Radiation schedule	$f(t)$	1 or 0
Diffusion coefficient (oxygen)	$D_k$	$2.5 \times 10^{-5}$ (cm $^2$ /s)	^[16]
Rate of oxygen supply	$r_k$	$0.1167$ (O $_2$ /s)	Estimate $^*$
Cellular oxygen consumption	$q_k$	$0.0414$ (O $_2$ /cells/s)	Estimate $^*$
Hill coefficient	$m$	2	Assumed $^*$
pO $_2$ at which the consumption is half the maximum rate	$k_c$	10 (mmHg)	Assumed $^*$
Diffusion coefficient (HAP)	$D_c$	$4 \times 10^{-5}$ (cm $^2$ /s)	Assumed $^*$
Production rate of HAP	$r_c$	$0.0704$ ([HAP]/s)	Assumed $^*$
HAP schedule	$g(t)$	1 or 0
Activation rate of HAP	$q_a$	$0.0084$ (1/s)	Assumed $^*$
Hill coefficient (activation)	$u$	2	Assumed $^*$
pO $_2$ at which the activation is half the maximum rate	$k_a$	5 (mmHg)	Assumed $^*$
HAP decay rate	$\lambda_C$	$0.023$ (1/s)	Assumed $^*$
Diffusion coefficient (AHAP)	$D_a$	$4 \times 10^{-5}$ (cm $^2$ /s)	Assumed $^*$
AHAP decay rate	$\lambda_a$	$0.023$ (1/s)	Assumed $^*$

| Show Table

DownLoad: CSV

Oxygen and drug distributions follow reaction diffusion equations and are assumed to be in quasi-steady state with respect to the computational time-step length used in numerical simulations. Specifically, their distributions are time-dependent over the course of the simulation, but achieve equilibrium within the span of a single computational time step of 0.001 days. Moreover, oxygen and HAP are assumed to be delivered from the blood vessels. Here, $v({\bf x}, t)$ stands for the density of the blood vessels and $v({\bf x}, t) = 1$ for the presence of the blood vessel and zero otherwise. Notice that a Hill function is again utilized for modelling the dependence of the oxygen consumption on the oxygen concentration. The function $g(t)$ , controls the delivery of HAP and when HAP is delivered, $g(t) = 1$ . The inactive HAPs ( $c({\bf x}, t)$ ) supplied by the vessels are activated at the presence of hypoxia and nitroreductases to produce activated HAPs ( $c_a({\bf x}, t)$ ). Once again a Hill function is used to mathematically show the dependence of the HAP activation rate on oxygen concentration. The parameter $q_a$ accounts for the rate by which HAPs are activated by nitroreductases. Here, for simplicity, we consider a uniform decay for both HAP and AHAP as we are not incorporating the details of drug dynamics. This could be improved with further experimental details and information. The active metabolites diffuse farther to produce a bystander effect affecting both hypoxic and oxygenated tumour cells. The formulation of the model is then completed by prescribing no-flux boundary conditions and an initial condition on the tumour cell density ^[16].

2.1. Computational domain and parameters

The simulations are performed using the images (Figure 1) obtained with the help of immunofluorescence staining techniques on a H460 xenograft to obtain vascularity, perfusion, and hypoxia using CD31 (endothelial marker), Hoechst 33342, and pimonidazole respectively, as described previously ^[11]. The stained tumour sections were then scanned and a threshold was applied to generate the corresponding binary masks.

Figure 1. Binary images of the xenographt cross-section, illustrating total tumour area, blood vessels, perfusion and hypoxia, respectively.

DownLoad: Full-Size Img PowerPoint

Following Powathil et al. ^[16], we assume that the perfused vascular network is the source of oxygen and drug supply, and thus gives the initial spatial distribution of concentrations. The perfused vasculature (at a fixed time) is obtained by combining the images of perfused areas (Hoechst 33342) and vascular structures (CD31) using the logical ``AND'' operation. The images of hypoxic regions (pimonidozole positive) are then used to compare the simulated hypoxic area and thus to calibrate/validate the mathematical model (Figure 1). The simulation is carried out over the image of the tumour cross-section (tumour mask). Note that for the in vitro case, we may assume a uniform distribution of oxygen and drug concentrations (when given). Further details can be found in ^[16].

The parameters for the model are listed in the Table 1. When available, we have used suitable parameters from previous studies ^[14,16]. The parameter values for oxygen dynamics have been calibrated to match the experimental results as shown in Figure 2. The experimental findings ^[24,25,26], suggest that HAPs have the ability to diffuse into hypoxic areas and exert their effect. Importantly, the diffusion parameters used in our model are estimates based on experimental observations and therefore are reflective of both the diffusive properties of the substances as well as their interactions with cells near the blood vessel. Production or supply rate of the drug is assumed to be (around 1.6 times) lower than that of oxygen and the decay rates of HAPs and AHAPs are assumed to be 0.023 (1/s). Moreover, we assume that HAPs are activated when the oxygen concentration is less than 5 mmHg ^[27]. We have further analyzed the effects of estimated parameter values using a sensitivity analysis, the results of which are reported in the relevant sections below. The initial spatial distribution of oxygen as determined from the vascular distribution is assumed to be constant with maximum value. Furthermore, we assume a uniform cell distribution as our initial condition, since we are considering a two dimensional cross section of in vivo xenografts.

Figure 2. Comparison of hypoxic images from (a) pimonidazole binding and (b) in silico model (white voxels in the in silico image correspond to hypoxic areas where

$pO_2 < 10\; mmHg$ ). The amount of hypoxic area predicted by the model differs by 3.06% from hypoxic area calculated by pimonidazole binding. Furthermore, the images have an overall 90.63% agreement in their pixel colours.

DownLoad: Full-Size Img PowerPoint

To our best of knowledge, most of the parameter values related to the evolution of HAPs and AHAPs are generally unknown or not clear and hence, these parameters are assumed phenomenologically based on the experimental observations and information. The clinical and experimental investigation into the effects and efficacy of HAPs alone or in combination with other treatment modalities is still in its infancy and consequently, modelling studies can play a vital role in providing a better overall qualitative or quantitative understanding. However, since these experimental and clinical studies are ongoing, the limited information available may not always help to accurately estimate all the necessary parameters needed for any modelling approaches. Nevertheless, they may help to give some qualitative, phenomenological estimates of the parameters that eventually help to provide qualitative inferences. Moreover, these modelling approaches can further inform experimental studies on accurately estimating appropriate parameters to provide quantitative insights. There are existing modelling works in this direction but each are focused on different compounds and settings and similarly use assumed values when accurate estimates are unavailable ^[15]. However, these assumed values can be modified or refined when accurate values of these parameters are estimated. For the values that are assumed in the current study, a sensitivity analysis is given to illustrate the sensitivity of the model results on the parameter values. We believe that unavailability of such data should not hinder modelling studies that might inform and motivate related experimental research; although, one should consider these results with some reservations until accurate parameters are available.

The partial differential equations governing each of the concentrations are numerically solved by using Finite Difference Method (FDM) over the tumour cross-section and the model is implemented using Matlab ^[28]. In addition, an iterative scheme is used to solve the equations for Oxygen, HAPs and AHAPs. Here, we assume that the growth is restricted within the cell-map obtained from images and the aim is to understand the effects of various treatment within this tumour cross section. Hence, for computational implementation, the tissue cross sections are embedded on a square domain with boundary conditions prescribed on the square domain. The dynamics of various distributions are only active where the cellmap () is non-zero. A uniform distribution is used for the initial condition of tumour cells across the mask and the other equations are solved in the steady state, so an initial condition is not necessary. The tumour cross-section is considered to be a 200 $\times$ 200 grid and is treated over 5 days with the simulation consisting of 5000 time steps each of size 0.001 days. The time step length remains constant throughout the treatment and particularly during the treatment times. The perfused vasculature network is read into the MATLAB script through the imread function.

3. Results and discussions

The identification of hypoxia in many types of tumours and its role in treatment resistance, especially in radiation therapy, prompted the exploration of other treatment options—in particular, HAPs ^[9]. However, as discussed earlier, recent clinical studies have failed to demonstrate the success of various HAP-based treatments. One possible reason for this failure could be lack of detailed information about the hypoxia levels of each patient and how that affects the success of the therapy ^[9]. In this section, following Powathil et al. ^[16], we discuss a theoretical, non-invasive approach to predict the spatial and temporal distribution of hypoxia (oxygen concentration) and analyze the effects of HAPs—alone and in combination with radiation.

In earlier work, Powathil et al. ^[16] discussed a computational approach to simulate and explore the spectrum of spatial hypoxia distributions at a snapshot in time that will result from a given vascular distribution. We used two-dimensional binary images of tumour cross-sections together with a vascular distribution as a computational domain on which a model for the oxygen distribution and tumour cell density was solved. The resulting hypoxic area was quantified and compared against the hypoxic proportions determined from images ^[16]. Here, since we are only considering viable tumour mass (tumour mask), necrotic regions are not included in the analysis. The methodology and results were discussed in detail in Powathil et al. ^[16]. Here, we use a similar approach to estimate the hypoxia on a xenograft sample and use it to study the effects of HAPs under varying hypoxic conditions. Figure 2 compares the binary images of hypoxia estimated using the mathematical model simulations and experiments (immunohistochemical analysis) and re-illustrates the usefulness of mathematical and theoretical approaches to predict hypoxic distributions. While the experimental and in silico tumours are not identical in their images, the model predicts hypoxic area to within 3.06% accuracy. As mathematical models are idealized by their nature, small discrepancies are expected. Furthermore, the model is still able to predict phenomena shown in experiments like efficacy based on hypoxic level and differences in sequencing.

3.1. HAPs and hypoxic levels

To study the role of spatial hypoxia distributions on HAP effectiveness, we used the above-described mathematical model to simulate the HAP response under two different hypoxic levels—low and high. In the low hypoxic case, we used all the perfused vessels obtained via image analysis of xenograft images ^[16] and in the high hypoxic case, we used a hypothetical scenario with only 50% of the perfused vessels. Figure 3 shows the comparison of total cell number and corresponding hypoxic area before and after 5 days of HAP treatment. In both cases, HAPs are given as daily doses over 5 days. We assume that in each dose, the drug is active in the system (delivered via vessels) for up to 3 hours. This can be changed or adapted appropriately depending on the given drug or mode of delivery. Since the primary aim of this study is to qualitatively study the effects of HAPs, we believe this schedule to be sufficient.

Figure 3. Plots showing the (a) total (normalized) tumour cell population and (b) (normalized) hypoxic area on day 0 and day 5, when the tumour is treated with HAP only (4 daily doses with active metabolites for up to 3 hours). Cell number is normalized using the initial cell number and hypoxic area is scaled as a percentage of total tumour area.

DownLoad: Full-Size Img PowerPoint

The simulation results shown in Figure 3 are qualitatively in agreement with experimental results ^[11]. The plots show that HAPs are more effective in reducing the number of cells under high hypoxic conditions: HAP treatment reduces cell number by 60.4% in the high hypoxic case and only by 31.1% in the low hypoxic case. Note that one possible reason for a relatively higher hypoxic area after treatment in simulation results might be due to the fact that experimental results account only for hypoxic area with viable cells (as binding is required), while simulation hypoxic area accounts for the hypoxic area within the entire domain (including regions with low cell density due to the cytotoxic effect of the released drug). This clearly indicates that information about hypoxic area can be a vital factor in determining the successfulness of treatments involving HAPs.

3.2. Combining HAPs with radiation

Radiation therapy is widely known to be very effective at targeting well-oxygenated tumour cells while hypoxic cells often show a degree of radiation resistance. Therefore, the biological rationale in combining radiation with HAPs is very promising as HAPs target hypoxic cells, leaving radio-sensitive tumour cells to be killed by radiation. Furthermore, recent preclinical experimental studies ^[8,29,11] suggest that the antitumour effectiveness of HAPs is enhanced by combining HAPs with radiotherapy. Here, we use the formulated mathematical model to study the effects of HAPs alone and in combination with radiation in low- and high-hypoxic tumours. Figure 4 shows the number of cells (Figure 4a) and hypoxic area (Figure 4b) before and after treatment under low- and high-hypoxic conditions. These results show that depending on the tumour type (hypoxia level), a combination treatment of radiation combined with HAPs can be very effective (compared to radiation alone).

Figure 4. Plots showing the (a) total (normalized) cell number and (b) (normalized) hypoxic area at the end of day 5 when the tumour is treated with radiation only or combination of radiation and HAP (8 Gy of radiation given over a 3 minute fraction on the first day of treatment with 4 daily doses of active metabolites for up to 3 hours). Cell number is normalized using the initial cell number and hypoxic area is scaled as a percentage of total tumour area.

DownLoad: Full-Size Img PowerPoint

3.3. Effects of treatment scheduling

Although combination therapies with HAPs have shown promising results in controlling a growing tumour, the positive response of this combination treatment often depends on the treatment sequencing and scheduling parameters. For example, Liu et al. showed in xenograft models that TH-302 in combination with a conventional chemotherapeutic achieved the greatest tumour inhibition when TH-302 was administered 2–8 hours prior to its companion (either cisplatin, docetaxel, gemcitabine, or doxorubicin) ^[10]. This work hypothesized that administration of the conventional therapeutic prior to TH-302 caused a reoxygenation of the tumour microenvironment, therefore decreasing the activation rate of the HAP and leading to decreased efficacy. Conversely, HAP-first treatment was thought to lead to a reoxygenation allowing further penetration of the accompanying drug. Peeters et al. examined the combination of TH-302 with radiotherapy applying TH-302 in different schedules with a single 8 Gy dose of radiation (as done in our simulations) ^[11]. They showed that TH-302 combined with radiation therapy increased tumour control from the radiation monotherapy case. Furthermore, it was shown that pretreatment with TH-302 prior to radiation lead to optimal efficacy. They speculated that HAP treatment reduced the hypoxic fraction which enhanced the radiation effectiveness. Nytko et al. showed the same qualitative results experimentally by combining radiation and HAP therapy, determining that a neoadjuvant (HAP-first) was the optimal treatment schedule ^[29].

Determining the optimal plan for administration of HAPs and radiation in a combination therapy is clearly an important goal. To help answer this question, the model was used to see the resulting tumour cell population and hypoxic area at the end of treatment for 5 different radiation schedules. In each schedule, HAPs are given on days 1 through 5 and the results are calculated at the end of the 5th day. An 8 Gy fraction of radiation is also given over a 3 minute interval on one of the first 4 days. Additionally, the case of splitting the dose into 4 fractions with 2 Gy/day was also included.

The analysis shows that the best option is to apply radiation on the 4th day of treatment, both in considering tumour cell population and hypoxic area. Though we only carried out this analysis for the 5-day treatment plan, we expect that this result would generalize to the best option simply being to apply radiation as late as possible. This is because the radiation has its highest effectiveness when hypoxic area is lowest: So naturally, the way to have the lowest hypoxic area is to apply the highest amount of HAPs before radiation. Similarly, the further fractionated 2 Gy/day schedule lies at a medium effectiveness. However this may not be the case for rapidly growing tumours where the turnover rate of hypoxic cells is very high. The results of this analysis are shown in Figure 5.

Figure 5. Plots showing the (a) total (normalized) cell number and (b) (normalized) hypoxic area on day 5 when the tumour is treated with Radiation and HAP. HAPs are given on days 1 through 4 and remain active for up to 3 hours. 8 Gy of radiation is given over a 3 minute fraction on one of the first 4 days except for the final case in which 2 Gy radiation is given on each of the first 4 days. The cell number and hypoxic area are calculated at the end of the 5th day. Cell number is normalized using the initial cell number and hypoxic area is scaled as a percentage of total tumour area.

DownLoad: Full-Size Img PowerPoint

4. Parameter sensitivity analysis

Here, we use a simple sensitivity analysis to study the influence of assumed/chosen parameters by comparing the results (cell number and hypoxia) obtained by halving and doubling the chosen values (default) of the parameters for two treatment scenarios.

The effectiveness of HAPs is also influenced by other relevant measurable parameters. For example, the diffusivity of the inactive HAPs, $D_c$ , impacts the delivery the drug to the hypoxic zones. To see this, we simulated the mathematical model with the drug diffusivity scaled by factors of 0.5 and 2 ( $0.5D_c$ and $2D_c$ ) and study the effect on resulting tumour cell population and hypoxic area. This was done in combination with radiation therapy for radiation given either on day 1 or day 4. The results show that the larger diffusivity results in a decrease in tumour cells and hypoxic area while the smaller diffusivity results in an increase in tumour cells and hypoxic area. However, the overall effect due to either change in diffusivity is low.

Another parameter which can be examined is the activation rate of the HAPs, $q_a$ . While this activation rate can't easily be pharmacologically altered, it is potentially alterable through the expression of nitroreductases in the tumour cells. With measurements of the levels of these enzymes, this activation rate could be theoretically estimated. Like for the diffusivity, we simulated the model for the activation rate scaled by factors of 0.5 and 2 for radiation given on either day 1 or day 4. The results show that the higher activation rate leads to a decrease in cell number and hypoxic area and a decrease in activation rate leads to the opposite. The overall effect of the change in activation rate is more significant than the effect of changing the diffusivity.

The effectiveness of the HAP also depends on the sensitivity of the activated HAPs in inducing cytotoxic effects to hypoxic and surrounding cells. For example, if it is an alkylating agent, the DNA repair status will be an important factor to consider while studying the effectiveness. DNA crosslinker such as the drug evofosfamide TH-302, pose strong replication blocks and cellular survival requires the concerted action of nucleotide excision, Fanconi anemia (FA) and homologous recombination (HR) driven DNA repair pathways ^[30,31]. The sensitivity to the drug is included into the model using the parameter $\delta_a$ and we have further explored its effects by varying this parameter. The results are reported in the Figure 5. The results show that increased sensitivity improves final outcome (higher cell-kill and lower hypoxic area), as expected.

A similar analysis was performed to study the sensitivity of each of the parameters related to the HAPs including: $r_c, \lambda_c, D_a$ , and $\lambda_a$ in addition to $D_c$ and $q_a$ as mentioned above. This analysis consisted of simulating the model with each of the parameters sequentially scaled by a factor of either 0.5 or 2 for radiation given on either the first or the fourth day. The parameters are scaled one at a time to assess to which parameters the model is most dependent. Looking at , we can see that $q_a, r_c, \lambda_c$ , and $\lambda_a$ are the parameters which have the largest impact on tumour cell number and hypoxic area. Though this method is not a rigorous sensitivity analysis, it nonetheless is a good indicator of the stability and robustness of our model since scales of parameters do not lead to drastic changes in the overall results. As precise measurements of the parameter values are not reliably known, especially since they can very patient-to-patient, we seek only to make qualitative predictions rather than quantitative ones. Accordingly, the goal of testing simple parameter scalings is to establish that the model is robust to changes of parameter values within a reasonable range. When precise measurements are able to be obtained, they can be used with the model to make more accurate predictions. The next step is thus to carefully design experimental studies to obtain relevant data such as HAP (drug specific) diffusivity and half-life to better understand and parameterize the mathematical model.

Figure 6. Plots showing the effect of scaling HAP diffusivity and activation rate on tumour cell population and hypoxic area for different radiation schedules. The top plot corresponds to HAPs given on each of days 1 through 4 which remain active up to 3 hours with radiation given on the first day. The bottom plot corresponds to the same HAP schedule with radiation given on day 4. The default parameters are given in Table 1. Cell number is normalized using the initial cell number and hypoxic area is scaled as a percentage of total tumour area.

DownLoad: Full-Size Img PowerPoint

5. Conclusion

We have developed here a spatio-temporal mathematical model to describe the effects of HAPs, alone and in combination with radiation, on a growing tumour. The model incorporates cell number, oxygen concentration, inactive HAPs, and activated metabolites. The system was solved on a domain obtained through imaging of perfused areas and vascular structures.

The model was solved for a low-hypoxia and high-hypoxia tumour to see which was more susceptible to treatment using HAPs. The model results showed a significant improvement in overall cell kill for tumours with high levels of hypoxia which was in agreement with experimental results ^[11]. This was extended to include a combination with radiation therapy which again predicted that the combination therapy is more effective for high-hypoxic tumours. Treatment sequencing was also altered to attempt to determine the optimal schedule of treatment methods. For our model, the best case scenario was to administer HAPs before radiation, which resulted in the highest cell kill and lowest hypoxic area. Additionally, scaling key model parameters can have an effect on the treatment efficacy. Simulating the results for different parameter scalings shows that increasing the diffusivity and activation rate of HAPs also improves the treatment efficacy, although the degree of change due to diffusivity scaling is relatively small.

Presently, the model is developed only using one cross-section of xenograft experimental image, showing its potential strength in studying the effects of HAPs in combination with radiation within a realistic microenvironmental setting. Further calibration and validation of the computational model using multiple images can help the model for further generalization. Also note that currently the model assumes a static vasculature over the treatment. Since the focus of this work is to study the effects of HAPs alone/in combination with radiation at a snapshot of the time with given spatial distribution of hypoxia resulted from a particular distribution of vasculature, the use of static vessels may be justified accordingly in a two dimensional settings. However, it would be ideal to simulate this in a three dimensional scenario, provided such spatial maps are available to implement the model. Furthermore, as the simulations take place over a maximum length of 5 days, the tumour vasculature is unlikely to change significantly over the course of treatment. For an in vivo tumour growing over a longer period of time, we would expect the vasculature to change throughout the treatment, thereby altering the transport of oxygen and drugs and effecting the treatment efficacy. This is another direction we hope to examine in the future. Another future direction is to incorporate a detailed pharmacodynamics of considered drug to improve the predictive ability of the model, moving a step closer to personalised therapeutic delivery. For example, for a drug such as PR-104, the intermediate molecules along its activation path may play a significant role as they too can have cytotoxic effects ^[13]. However, for a drug such as TH-302, as we considered here, bromo-isophosphoramide mustard (Br-IPM) is the lone dominant cytotoxic form (though the Cl-IPM metabolite also appears) ^[32]. Similarly, the simplified model that we used here to study radiation effects can be modified further to include multiple effects of radiation such as delayed cell death, indirect effects, cell-cycle sensitivity and so on, depending on the availability of the data and required complexity of the model. Spatial effects are also not examined to their full extent. Analysis of treatment efficacy as it relates to distance from vessels is an interesting question which we hope to examine in future works. As we use only one experimental image to generate the hypoxia distribution, the conclusions we can draw here are limited.

Although, current predictions from the model are qualitative in nature, we have shown that the model can be used to understand the combined effects of radiation and HAPs in general. Moreover, the results from the model are in good agreement with experimental observations and the inferences from the model results such as optimal combinations and sequencing, can be further analyzed and validated using appropriate experimental studies. Once the model is carefully calibrated with HAP specific experimental data, it can be used for quantitative predictions and treatment optimizations. Furthermore, this model can be adapted for clinical applications by using relevant hypoxic data obtained from biopsies and novel imaging techniques such as perfusion CT imaging to spatially reconstruct tumour hypoxia map for the model.

Acknowledgements

Financial supports by the ERC advanced grant (ERC-ADG-2015, n 694812 Hypoximmuno) and Kankeronderzoekfonds Limburg from the Health Foundation Limburg (AY, LJD, and PL), and the Natural Sciences and Engineering Research Council of Canada (NSERC) (MK) are gratefully acknowledged.

Conflicts of interest

All authors declare no conflicts of interest in this paper.

References

[1]	C. M. Lalli, T. R. Parsons, Biological oceanography: An introduction, Butterworth-Heinemann, Oxford, 1997.
[2]	S. J. Brentnall, K. J. Richards, J. Brindley, E. Murphy, Plankton patchiness and its effect on larger-scale productivity, J. Plankton Res., 25 (2003), 121–140. https://doi.org/10.1093/plankt/25.2.121 doi: 10.1093/plankt/25.2.121
[3]	H. W. Paerl, N. S. Hall, E. S. Calandrino, Controlling harmful cyanobacterial blooms in a world experiencing anthropogenic and climatic-induced change, Sci. Total Environ., 409 (2011), 1739–45. https://doi.org/10.1016/j.scitotenv.2011.02.001 doi: 10.1016/j.scitotenv.2011.02.001
[4]	H. W. Paerl, T. G. Otten, Harmful cyanobacterial blooms: causes, consequences, and controls, Microb. Ecol., 65 (2013), 995–1010. https://doi.org/10.1007/s00248-012-0159-y doi: 10.1007/s00248-012-0159-y
[5]	D. L. Roelke, J. P. Grover, B. W. Brooks, J. Glass, D. Buzan, G. M. Southard, et al., A decade of fish-killing prymnesium parvum blooms in texas: roles of inflow and salinity, J. Plankton Res., 33 (2010), 243–253. https://doi.org/10.1093/plankt/fbq079 doi: 10.1093/plankt/fbq079
[6]	M. J. Ulloa, P. Alvarez-Torres, K. P. Horak-Romo, R. Ortega-Izaguirre, Harmful algal blooms and eutrophication along the mexican coast of the gulf of mexico large marine ecosystem, Environ. Dev., 22 (2017), 120–128. https://doi.org/10.1016/j.envdev.2016.10.007 doi: 10.1016/j.envdev.2016.10.007
[7]	J. H. Landsberg, The effects of harmful algal blooms on aquatic organisms, Rev. Fish Sci., 10 (2002), 113–390. https://doi.org/10.1080/20026491051695 doi: 10.1080/20026491051695
[8]	K. W. Crane, J. P. Grover, Coexistence of mixotrophs, autotrophs, and heterotrophs in planktonic microbial communities, J. Theor. Biol., 262 (2010), 517–527. https://doi.org/10.1016/j.jtbi.2009.10.027 doi: 10.1016/j.jtbi.2009.10.027
[9]	H. Stickney, R. Hood, D. Stoecker, The impact of mixotrophy on planktonic marine ecosystems, Ecol. Modell., 125 (2000), 203–230. https://doi.org/10.1016/S0304-3800(99)00181-7 doi: 10.1016/S0304-3800(99)00181-7
[10]	J. Zhang, J. D. Kong, J. Shi, H. Wang, Phytoplankton competition for nutrients and light in a stratified lake: A mathematical model connecting epilimnion and hypolimnion, J. Nonlinear Sci., 31 (2021), 35. https://doi.org/10.1007/s00332-021-09693-6 doi: 10.1007/s00332-021-09693-6
[11]	M. Chen, M. Fan, R. Liu, X. Wang, X. Yuan, H. Zhu, The dynamics of temperature and light on the growth of phytoplankton, J. Theor. Biol., 385 (2015), 8–19. https://doi.org/10.1016/j.jtbi.2015.07.039 doi: 10.1016/j.jtbi.2015.07.039
[12]	K. Yoshiyama, H. Nakajima, Catastrophic transition in vertical distributions of phytoplankton: Alternative equilibria in a water column, J. Theor. Biol., 216 (2002), 397–408., https://doi.org/10.1006/jtbi.2002.3007 doi: 10.1006/jtbi.2002.3007
[13]	X. Zhao, L. Liu, H. Wang, M. Fan, Ecological effects of predator harvesting and environmental noises on oceanic coral reefs, Bull. Math. Biol., 85 (2023), 59. https://doi.org/10.1007/s11538-023-01166-z doi: 10.1007/s11538-023-01166-z
[14]	X. Zhao, L. Liu, M. Liu, M. Fan, Stochastic dynamics of coral reef system with stage-structure for crown-of-thorns starfish, Chaos Solitons Fractals, 181 (2024), 1–23. https://doi.org/10.1016/j.chaos.2024.114629 doi: 10.1016/j.chaos.2024.114629
[15]	C. C. Carey, K. L. Cottingham, K. C. Weathers, J. A. Brentrup, N. M. Ruppertsberger, H. A. Ewing, Experimental blooms of the cyanobacterium gloeotrichia echinulata increase phytoplankton biomass, richness and diversity in an oligotrophic lake, J. Plankton Res., 36 (2014), 364–377. https://doi.org/10.1093/plankt/fbt105 doi: 10.1093/plankt/fbt105
[16]	Y. Lehahn, I. Koren, S. Sharoni, F. d'Ovidio, A. Vardi, E. Boss, Dispersion/dilution enhances phytoplankton blooms in low-nutrient waters, Nat. Commun., 8 (2017), 14868. https://doi.org/10.1038/ncomms14868 doi: 10.1038/ncomms14868
[17]	R. W. Sterner, K. L. Reinl, B. M. Lafrancois, S. Brovold, T. R. Miller, A first assessment of cyanobacterial blooms in oligotrophic lake superior, Limnol. Oceanogr., 65 (2020), 2984–2998. https://doi.org/10.1002/lno.11569 doi: 10.1002/lno.11569
[18]	J. Zhang, J. Shi, X. Chang, A mathematical model of algae growth in a pelagic-benthic coupled shallow aquatic ecosystem, J. Math. Biol., 76 (2018), 1159–1193. https://doi.org/10.1007/s00285-017-1168-8 doi: 10.1007/s00285-017-1168-8
[19]	J. N. Boyer, S. K. Dailey, P. J. Gibson, M. T. Rogers, D. Mir-Gonzalez, The role of dissolved organic matter bioavailability in promoting phytoplankton blooms in florida bay, Hydrobiologia, 569 (2006), 71–85. https://doi.org/10.1007/s10750-006-0123-2 doi: 10.1007/s10750-006-0123-2
[20]	C. H. Chow, W. Cheah, J. H. Tai, A rare and extensive summer bloom enhanced by ocean eddies in the oligotrophic western north pacific subtropical gyre, Sci. Rep., 7 (2017), 6199. https://doi.org/10.1038/s41598-017-06584-3 doi: 10.1038/s41598-017-06584-3
[21]	V. Ligorini, N. Malet, M. Garrido, V. Derolez, M. Amand, B. Bec, et al., Phytoplankton dynamics and bloom events in oligotrophic mediterranean lagoons: seasonal patterns but hazardous trends, Hydrobiologia, 849 (2022), 2353–2375. https://doi.org/10.1007/s10750-022-04874-0 doi: 10.1007/s10750-022-04874-0
[22]	S. J. Thackeray, I. Maberly, Long-term change in the phenology of spring phytoplankton: Species-specific responses to nutrient enrichment and climatic change, J. Ecol., 96 (2008), 523–535. https://doi.org/10.1111/j.1365-2745.2008.01355.x doi: 10.1111/j.1365-2745.2008.01355.x
[23]	F. B. Wang, S. B. Hsu, W. Wang, Dynamics of harmful algae with seasonal temperature variations in the cove-main lake, Discr. Contin. Dyn. Syst. B, 21 (2016), 313–335. https://doi.org/10.3934/dcdsb.2016.21.313 doi: 10.3934/dcdsb.2016.21.313
[24]	A. Huppert, B. Blasius, R. Olinky, L. Stone, A model for seasonal phytoplankton blooms, J. Theor. Biol., 236 (2005), 276–290. https://doi.org/10.1016/j.jtbi.2005.03.012 doi: 10.1016/j.jtbi.2005.03.012
[25]	M. Chen, M. Fan, X. Yuan, R. Liu, X. Wang, H. Zhu, Effect of seasonal changing temperature on the growth of phytoplankton, Math. Biosci. Eng., 14 (2017), 1091. https://doi.org/10.3934/mbe.2017057 doi: 10.3934/mbe.2017057
[26]	L. Ratnarajah, R. Abu-Alhaija, A. Atkinson, S. Batten, N. J. Bax, K. S. Bernard, et al., Monitoring and modelling marine zooplankton in a changing climate, Nat. Commun., 14 (2023), 564. https://doi.org/10.1038/s41467-023-36241-5 doi: 10.1038/s41467-023-36241-5
[27]	S. G. Pitois, C. P. Lynam, T. Jansen, N. Halliday, M. Edwards, Bottom-up effects of climate on fish populations: data from the continuous plankton recorder, Mar. Ecol. Prog. Ser., 456 (2012), 169–186. https://doi.org/10.3354/meps09710 doi: 10.3354/meps09710
[28]	J. J. Ruzicka, R. D. Brodeur, R. L. Emmett, J. H. Steele, J. E. Zamon, C. A. Morgan, et al., Interannual variability in the northern california current food web structure: Changes in energy flow pathways and the role of forage fish, euphausiids, and jellyfish, Prog. Oceanogr., 102 (2012), 19–41. https://doi.org/10.1016/j.pocean.2012.02.002 doi: 10.1016/j.pocean.2012.02.002
[29]	V. Lauria, M. J. Attrill, A. Brown, M. Edwards, S. C. Votier, Regional variation in the impact of climate change: evidence that bottom-up regulation from plankton to seabirds is weak in parts of the Northeast Atlantic, Mar. Ecol. Prog. Ser., 488 (2013), 11–22. https://doi.org/10.3354/meps10401 doi: 10.3354/meps10401
[30]	R. F. Heneghan, J. D. Everett, J. L. Blanchard, A. J. Richardson, Zooplankton are not fish: Improving zooplankton realism in size-spectrum models mediates energy transfer in food webs, Front. Mar. Sci., 3 (2016). https://doi.org/10.3389/fmars.2016.00201
[31]	J. Luo, Phytoplankton-zooplankton dynamics in periodic environments taking into account eutrophication, Math. Biosci., 245 (2013), 126–136. https://doi.org/10.1016/j.mbs.2013.06.002 doi: 10.1016/j.mbs.2013.06.002
[32]	P. Lehette, A. Tovar-Sanchez, C. M. Duarte, S. Hernandez-Leon, Krill excretion and its effect on primary production, Mar. Ecol. Prog. Ser., 459 (2012), 29–38. https://doi.org/10.3354/meps09746 doi: 10.3354/meps09746
[33]	E. L. Cavan, A. Belcher, A. Atkinson, S. L. Hill, S. Kawaguchi, S. McCormack, et al., The importance of antarctic krill in biogeochemical cycles, Nat. Commun., 10 (2019), 4742. https://doi.org/10.1038/s41467-019-12668-7 doi: 10.1038/s41467-019-12668-7
[34]	K. Schmidt, C. Schlosser, A. Atkinson, S. Fielding, H. J. Venables, C. M. Waluda, et al., Zooplankton gut passage mobilizes lithogenic iron for ocean productivity, Curr. Biol., 26 (2016), 2667–2673. https://doi.org/10.1016/j.cub.2016.07.058 doi: 10.1016/j.cub.2016.07.058
[35]	J. J. Elser, J. Urabe, The stoichiometry of consumer-driven nutrient recycling: theory, observations, and consequences, Ecology, 80 (1999), 735–751. https://doi.org/10.1890/0012-9658(1999)080 doi: 10.1890/0012-9658(1999)080
[36]	X. Irigoien, K. J. Flynn, R. P. Harris, Phytoplankton blooms: a 'loophole' in microzooplankton grazing impact?, J. Plankton Res., 27 (2005), 313–321. https://doi.org/10.1093/plankt/fbi011 doi: 10.1093/plankt/fbi011
[37]	J. Ji, R. Milne, H. Wang, Stoichiometry and environmental change drive dynamical complexity and unpredictable switches in an intraguild predation model, J. Math. Biol., 86 (2023), 31. https://doi.org/10.1007/s00285-023-01866-z doi: 10.1007/s00285-023-01866-z
[38]	P. D. Jeyasingh, J. M. Goos, S. K. Thompson, C. M. Godwin, J. B. Cotner, Ecological stoichiometry beyond redfield : An ionomic perspective on elemental homeostasis, Front. Microbiol., 8. https://doi.org/10.3389/fmicb.2017.00722
[39]	A. Peace, Effects of light, nutrients, and food chain length on trophic efficiencies in simple stoichiometric aquatic food chain models, Ecol. Modell., 312 (2015), 125–135. https://doi.org/10.1016/j.ecolmodel.2015.05.019 doi: 10.1016/j.ecolmodel.2015.05.019
[40]	I. Loladze, Y. Kuang, J. Elser, Stoichiometry in producer-grazer systems: Linking energy flow with element cycling, Bull. Math. Biol., 62 (2000), 1137–1162. https://doi.org/10.1006/bulm.2000.0201 doi: 10.1006/bulm.2000.0201
[41]	A. Edwards, J. Brindley, Zooplankton mortality and the dynamical behaviour of plankton population models, Bull. Math. Biol., 61 (1999), 303–339. https://doi.org/10.1006/bulm.1998.0082 doi: 10.1006/bulm.1998.0082
[42]	H. Wang, R. W. Sterner, J. J. Elser, On the 'strict homeostasis' assumption in ecological stoichiometry, Ecol. Model., 243 (2012), 81–88. https://doi.org/10.1016/j.ecolmodel.2012.06.003 doi: 10.1016/j.ecolmodel.2012.06.003
[43]	M. Chen, M. Fan, Y. Kuang, Global dynamics in a stoichiometric food chain model with two limiting nutrients, Math. Biosci., 289 (2017), 9–19. https://doi.org/10.1016/j.mbs.2017.04.004 doi: 10.1016/j.mbs.2017.04.004
[44]	D. Song, M. Fan, M. Chen, H. Wang, Dynamics of a periodic stoichiometric model with application in predicting and controlling algal bloom in bohai sea off china, Math. Biosci. Eng., 16 (2019), 119–138. https://doi.org/10.3934/mbe.2019006 doi: 10.3934/mbe.2019006
[45]	K. Mischaikow, H. Smith, H. R. Thieme, Asymptotically autonomous semiflows: chain recurrence and lyapunov functions, Trans. Amer. Math. Soc., 347 (1995), 1669–1685. https://doi.org/10.1090/S0002-9947-1995-1290727-7 doi: 10.1090/S0002-9947-1995-1290727-7
[46]	R. E. Gaines, J. L. Mawhin, Coincidence degree and nonlinear differential equations, Springer, 1977.
[47]	I. Barhalat, Systems d'equations differential d'oscillations nonlinearies, Rev. Roum. Math. Pures. Appl., 4 (1959), 267–270.

This article has been cited by:

Peiqun Lin, Chenxing He, Lingshu Zhong, Mingyang Pei, Chuhao Zhou, Yang Liu, Bus timetable optimization model in response to the diverse and uncertain requirements of passengers for travel comfort, 2023, 31, 2688-1594, 2315, 10.3934/era.2023118

Reader Comments

Your name:*

Email:*
© 2024 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)