A Mathematical Model of Cysteine-Driven Metabolic Adaptation to Hypoxia in Ovarian Cancer

TL;DR

This paper presents a mathematical model showing how cysteine metabolism helps ovarian cancer cells adapt to low oxygen conditions.

Contribution

A novel reduced mechanistic model integrating cysteine allocation, redox control, and energy maintenance under hypoxia in ovarian cancer.

Findings

The model accurately reproduces extracellular fluxes with deviations below 7%.

Glutathione synthesis capacity is the main factor maintaining ATP under hypoxia.

Hydrogen sulfide production has a secondary stabilizing effect on energy state.

Abstract

Ovarian cancer progression is strongly influenced by tumour hypoxia and associated oxidative stress. Experimental evidence indicates that cysteine availability supports ovarian cancer cell fitness under hypoxic conditions, yet the quantitative integration of cysteine metabolism, redox control, and energetic maintenance remains incompletely understood. We present a reduced mechanistic mathematical model describing intracellular cysteine allocation between glutathione (GSH) synthesis and hydrogen sulfide production under experimentally imposed hypoxia. The model integrates extracellular cysteine uptake, GSH-dependent reactive oxygen species (ROS) detoxification, hypoxia-amplified ROS generation, and redox-modulated ATP maintenance. Parameter estimation was performed using experimentally derived extracellular metabolite fluxes measured over a 24 h interval. Uncertainty was assessed via bootstrap resampling, and variance-based sensitivity analysis was conducted within (patho)physiologically constrained parameter domains. The calibrated model reproduces extracellular fluxes with relative deviations below 7% and identifies GSH synthesis capacity as the dominant determinant of ATP maintenance within experimentally supported ranges. Hydrogen sulfide (H2S) production exerts a secondary stabilising influence, whereas hypoxia-driven ROS amplification negatively impacts energetic state. Numerical continuation across hypoxia levels reveals distinct qualitative response regions but does not imply a formal bifurcation structure. Importantly, intracellular metabolite dynamics are inferred as latent variables consistent with extracellular constraints and established biochemical knowledge; the model does not uniquely identify intracellular pool sizes or enzyme kinetics. The framework therefore provides flux-consistent mechanistic plausibility rather than direct intracellular validation. This systems-level analysis supports cysteine allocation as a quantitatively influential control point in hypoxic adaptation and establishes a constrained modelling framework for subsequent metabolic network expansions and experimental validation.

Figures

Equations

Full text

1. Introduction

Ovarian cancer remains one of the deadliest gynaecological malignancies, with survival rates that have improved only modestly over recent decades [1]. A hallmark of advanced ovarian tumours is the presence of hypoxic regions, arising from abnormal vasculature and rapid tumour growth. Hypoxia is a potent driver of tumour progression, metabolic reprogramming, and resistance to therapy [2,3].

One major consequence of hypoxia is increased oxidative stress, resulting from impaired mitochondrial electron transport and altered redox homeostasis. To survive under such conditions, cancer cells activate adaptive metabolic programmes that enhance antioxidant capacity and sustain energy production [4].

Cysteine plays a central role in these adaptations. It is the rate-limiting precursor for glutathione (GSH), the major intracellular antioxidant, and a substrate for hydrogen sulfide (H_2_S) production, which has been implicated in redox signalling and mitochondrial bioenergetics [5]. Experimental studies have demonstrated that ovarian cancer cells increase their dependence on cysteine under hypoxia, promoting survival and resistance to platinum-based chemotherapy [6,7].

Despite these insights, a quantitative understanding of how cysteine-dependent metabolic pathways interact to determine cellular fitness under hypoxia is still lacking. Mathematical modelling provides a powerful approach to integrate experimental data, formalise mechanistic hypotheses, and identify key metabolic control points.

Here, we develop a mechanistic mathematical model describing cysteine-driven metabolic adaptation to hypoxia in ovarian cancer cells. The model is calibrated using experimental metabolomic data and analysed through numerical simulation and global sensitivity analysis.

2. Experimental Data and Biological Basis

The model is informed by the levels of extracellular metabolites obtained from the exometabolome of ovarian cancer cell lines (ES2; CRL-1978, and OVCAR-3; HTB-161 from American Type Culture Collection-ATCC) cultured under normoxic and hypoxic conditions, with and without cysteine supplementation. The culture conditions and the experiments needed to obtain the exometabolome profile are already published [7]. Metabolite concentrations were measured in the culture medium at defined time points, including amino acids, such as cysteine, methionine, glutamate, glutamine, and glucose.

Changes in extracellular concentrations over a 24 h period were used to infer apparent uptake and secretion fluxes according to

$$[eqn]$$

where $[eqn]$ denotes the extracellular concentration of metabolite i.

The experimental dataset used in this study consists of extracellular metabolite concentrations measured at the beginning (0 h) and end (24 h) of the culture period under defined experimental conditions, including hypoxia and cysteine supplementation. These data provide quantitative information on the net metabolite uptake and secretion rates over the experimental time window.

For each measured extracellular metabolite i, an apparent mean flux was computed as

$$[eqn]$$

where concentrations are expressed in mM and time in hours. These apparent fluxes constitute the primary experimental observables used for quantitative calibration of the mathematical model.

Intracellular metabolite concentrations are not directly measured and are therefore treated as latent variables, inferred indirectly through consistency between the model-predicted fluxes and experimentally observed extracellular fluxes.

These data indicate enhanced cysteine consumption under hypoxia, accompanied by coordinated changes in redox-related metabolites, supporting the central role of cysteine in hypoxic adaptation.

3. Mathematical Model Formulation

3.1. State Variables

The model describes the temporal evolution of five intracellular variables:

• $[eqn]$ : intracellular cysteine concentration,
• $[eqn]$ : reduced glutathione (GSH),
• $[eqn]$ : reactive oxygen species (ROS),
• $[eqn]$ : hydrogen sulfide (H_2_S),
• $[eqn]$ : ATP concentration (proxy for cellular energetic state).

Hypoxia is represented by a dimensionless parameter $[eqn]$ , corresponding to normoxia and hypoxia, respectively. In this first modelling framework, hypoxia is treated as a fixed external condition, reflecting experimentally imposed oxygen levels rather than dynamically evolving oxygen availability. While oxygen concentration may vary in space and time in vivo, explicitly modelling oxygen transport and consumption is beyond the scope of the present work and will be addressed in future extensions of the model.

3.2. Model Equations

The system is formulated as a set of ordinary differential equations:

• Cysteine dynamics

$$[eqn]$$

where cysteine is taken up from the extracellular medium and consumed by glutathione and H_2_S synthesis, as well as basal metabolic processes.

Glutathione dynamics

$$[eqn]$$

Hypoxia enhances the flux of cysteine toward GSH synthesis, reflecting increased antioxidant demand.

ROS dynamics

$$[eqn]$$

The term $[eqn]$ represents mass-action GSH-dependent ROS neutralisation, where the detoxification rate scales proportionally with both reduced glutathione concentration and ROS abundance. The explicitly multiplicative structure removes notational ambiguity and ensures direct correspondence between model equations and biochemical interaction terms.

Hydrogen sulfide dynamics

$$[eqn]$$

ATP dynamics

$$[eqn]$$

ATP production is not assumed to be directly generated by GSH. Instead, GSH is modelled as a permissive regulator of mitochondrial function through its role in redox buffering. The term $[eqn]$ captures the experimentally supported requirement of adequate redox control for efficient oxidative phosphorylation. Thus, G acts as a modulatory factor rather than as a direct energetic substrate.

This formulation separates oxygen-dependent ATP production from hypoxia-tolerant H_2_S-supported energetic pathways, reflecting experimental evidence that H_2_S can contribute to mitochondrial electron transport even under limited oxygen availability.

It should be noted that the structure of the glutathione equation favours stability of GSH levels provided cysteine is available, reflecting a strong homeostatic control. The robustness of this behaviour is further explored numerically in Section 6.

The system of Equations (3)–(7) defines a reduced intracellular metabolic network in which cysteine serves as a branching substrate between two functionally distinct pathways: redox buffering through GSH synthesis and alternative mitochondrial support through H_2_S production. Hypoxia is incorporated as an external parametric driver that increases ROS generation and redox demand while simultaneously impairing oxygen-dependent ATP production.

The structural organisation of this reduced network is summarised in Figure 1. The schematic provides a graphical correspondence between the biochemical processes and the mathematical terms appearing in the differential equations. In particular, it visualises cysteine allocation, GSH-dependent ROS neutralisation, H_2_S-supported mitochondrial flux, and hypoxia-mediated modulation. No additional pathways are introduced beyond those explicitly represented in Equations (3)–(7); the diagram therefore functions as a structural abstraction of the model rather than as an expanded biological map.

All state variables (C, G, R, S, A) are expressed in mM. Time is expressed in h. Uptake and first-order turnover parameters ( $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ ) have units of h^−1^ when multiplied by concentrations. Maximal flux parameters ( $[eqn]$ , $[eqn]$ ) have units of mM h^−1^, and Michaelis constants ( $[eqn]$ , $[eqn]$ ) have units of mM. Each term in Equations (3)–(7) therefore has units of mM h^−1^, ensuring dimensional consistency throughout the dynamical system.

The dimensional homogeneity of each equation was independently verified by expressing all flux terms in units of mM h^−1^. No composite term combines parameters and variables with inconsistent units.

4. Parameter Estimation and Numerical Methods

4.1. Selection of Parameters for Calibration

Parameter estimation was restricted to a subset of model parameters that are directly informed by the available experimental data. Specifically, the following parameters were estimated from the data: the cysteine uptake rate $[eqn]$ , the maximal flux of cysteine toward glutathione synthesis $[eqn]$ , the maximal flux of cysteine toward hydrogen sulfide production $[eqn]$ , the basal rate of reactive oxygen species generation $[eqn]$ , and the hypoxia-induced amplification factor $[eqn]$ .

All remaining parameters were fixed to biologically plausible values obtained from the literature. This parameter selection strategy reduces model dimensionality, improves numerical identifiability, and ensures that estimated parameters are constrained by experimental observations.

4.2. Definition of Experimental Observables

The experimental observables used for model calibration are the apparent extracellular uptake and secretion fluxes computed from concentration changes over a 24 h period. For each experimental condition, model simulations were performed over the same time interval, and the mean model-predicted fluxes were computed as

$$[eqn]$$

These quantities were directly compared with the experimentally inferred fluxes $[eqn]$ , ensuring consistency between the numerical model and the experimental dataset.

4.3. Objective Function for Parameter Estimation

Parameter estimation was formulated as a constrained optimisation problem, in which the model parameter vector $[eqn]$ was determined by minimising an objective function that quantified the mismatch between experimentally inferred extracellular metabolite fluxes and those predicted by the mathematical model.

For each extracellular metabolite $[eqn]$ and each experimental condition $[eqn]$ , the experimental observable is the apparent mean flux $[eqn]$ , computed from extracellular concentration changes over a 24 h interval. Given a parameter set $[eqn]$ , the model was numerically integrated over the same time window, and the corresponding mean model-predicted flux was computed as

$$[eqn]$$

The objective function was defined as

$$[eqn]$$

and parameter estimation was performed by minimising $[eqn]$ over a biologically admissible parameter domain.

The use of relative squared errors ensures that all metabolites contribute comparably to the objective function, independently of their absolute concentration levels or units, thereby preventing dominant contributions from high-abundance metabolites. This formulation is particularly appropriate for extracellular metabolomic data, where concentration ranges may differ by several orders of magnitude across metabolites.

By calibrating the model to experimentally inferred fluxes rather than intracellular concentrations, the optimisation problem remains fully consistent with the available data while allowing intracellular metabolite dynamics to be inferred implicitly. This flux-based calibration strategy reduces sensitivity to unknown intracellular pool sizes and measurement noise, and yields parameter estimates that are directly anchored to experimentally observable quantities.

Finally, the structure of the objective function allows for direct comparison across experimental conditions, including normoxia versus hypoxia and cysteine-supplemented versus control media, thereby ensuring that the estimated parameters capture the coordinated metabolic adaptations observed experimentally.

4.4. Numerical Optimisation and Uncertainty Analysis

Minimisation of the objective function defined in Section 4.3 was performed using nonlinear least-squares optimisation. To mitigate the risk of convergence to local minima, the optimisation procedure was repeated from multiple initial parameter guesses sampled uniformly within predefined biologically admissible ranges. These ranges were selected on the basis of biochemical plausibility and prior knowledge reported in the literature, particularly regarding hypoxia-induced metabolic regulation and cysteine-centred redox control [2,5].

Parameter bounds were chosen to ensure the positivity and physiological relevance of all estimated quantities, while remaining sufficiently wide to avoid artificial constraint of the optimisation landscape. This approach is particularly important in nonlinear metabolic models, where parameter interdependencies may give rise to flat or weakly identifiable directions in parameter space.

Uncertainty in the estimated parameter values was assessed using a bootstrap resampling strategy based on experimental flux variability. Specifically, experimental extracellular fluxes were resampled assuming independent perturbations consistent with the observed experimental variability, and the calibration procedure was repeated for each resampled dataset. The resulting ensemble of parameter estimates was used to compute confidence intervals and summary statistics for each calibrated parameter.

To further characterise the parameter dependencies, pairwise correlations between estimated parameters were computed from the bootstrap ensemble. This analysis was used to identify potential compensatory mechanisms within the model, in particular between the maximal flux of cysteine toward glutathione synthesis and the maximal flux toward hydrogen sulfide production. Such compensatory behaviour is biologically plausible, as both pathways draw from the same intracellular cysteine pool and have been implicated in hypoxia-associated redox and metabolic adaptation [5].

Together, this optimisation and uncertainty analysis framework ensures that the reported parameter estimates are not only consistent with experimental observations but also quantitatively robust and interpretable within the known biological constraints of hypoxia-driven metabolic reprogramming [2].

While the multi-start optimisation strategy ensures convergence towards a consistent minimising parameter vector, numerical agreement alone does not guarantee the practical recoverability of individual parameters. In nonlinear metabolic models constrained by limited observables, assessment of the parameter dispersion under experimental variability is therefore required to evaluate practical identifiability. This motivates the bootstrap-based analysis presented below.

4.5. Practical Identifiability Analysis

Because parameter estimation relies on extracellular flux measurements evaluated at two time points (0 h and 24 h), assessment of practical identifiability is essential. To quantify the estimation robustness, a nonparametric bootstrap procedure was implemented by resampling the experimentally inferred extracellular fluxes within their observed variability and recalibrating the model for each synthetic dataset.

Let $[eqn]$ denote the parameter vector estimated from bootstrap replicate b, with $[eqn]$ . In the present study, $[eqn]$ bootstrap replicates were generated. For each calibrated parameter $[eqn]$ , empirical 95% confidence intervals were computed using the percentile method, yielding

$$[eqn]$$

The resulting bootstrap-based summary statistics are reported in Table 1. All calibrated parameters exhibit bounded and non-degenerate confidence intervals.

The coefficients of variation remain below 8% for all calibrated parameters, indicating practical identifiability within the flux-constrained reduced model. The moderate negative correlation between $[eqn]$ and $[eqn]$ ( $[eqn]$ ) reflects the partial compensatory allocation of intracellular cysteine flux. However, no near-singular parameter dependencies (defined as $[eqn]$ ) were detected, indicating the absence of severe compensation or practical non-identifiability within the biologically admissible parameter domain.

The hypoxia amplification parameter $[eqn]$ exhibits particularly low dispersion across bootstrap replicates. This reflects the strong constraint imposed by the hypoxia-dependent ROS amplification term within the flux-calibrated objective function. The narrow interval does not arise from boundary saturation, as the optimisation bounds allow substantially wider admissible variation.

From a structural perspective, identifiability of the hypoxia amplification parameter $[eqn]$ follows directly from the steady-state form of the ROS equation. At equilibrium,

$$[eqn]$$

so that, for $[eqn]$ ,

$$[eqn]$$

Hence, $[eqn]$ is algebraically recoverable from the steady-state redox balance once $[eqn]$ is independently estimated, and simulations include at least two distinct hypoxia levels (e.g., normoxia $[eqn]$ and hypoxia $[eqn]$ ). The narrow bootstrap interval therefore reflects the strong structural constraint imposed by the redox steady-state relation rather than optimisation artefacts or boundary saturation effects.

To further characterise the local curvature properties of the objective function near the optimum, the Hessian matrix of $[eqn]$ was numerically approximated at the minimising parameter vector. The resulting condition number ( $[eqn]$ ) indicates anisotropic but well-conditioned local curvature, with no evidence of flat directions compatible with practical ill-posedness in the calibrated parameter neighbourhood.

Taken together, these results demonstrate controlled parameter compensation without degeneracy of the inverse problem under the imposed cysteine-centred reduction architecture.

5. Sensitivity Analysis

Sensitivity analysis was conducted to quantify how parameter uncertainty propagates to model outputs associated with energetic and redox state. Variance-based Sobol indices were computed within physiologically admissible parameter domains. Sampling was performed using quasi-random Sobol sequences across bounded intervals reflecting both biochemical plausibility and bootstrap-derived uncertainty ranges. Unlike purely structural global sensitivity analysis over arbitrarily wide domains, the present approach quantifies variance propagation within experimentally constrained regions of parameter space. The resulting indices therefore represent biologically relevant control coefficients rather than unrestricted structural sensitivities.

5.1. Sensitivity Method

Global sensitivity analysis was performed using first-order Sobol indices to quantify the relative contribution of individual model parameters to the variance of selected model outputs. This variance-based approach is well suited for nonlinear dynamical systems, as it captures parameter influence over the full parameter range while accounting for non-additive effects.

Sensitivity analysis was carried out with respect to the subset of parameters calibrated against the experimental data, while all remaining parameters were fixed at their nominal values. Parameter sampling was performed around the calibrated parameter vector, with each parameter independently varied within a restricted interval of $[eqn]$ – $[eqn]$ of its estimated value. These bounds were chosen to reflect the experimentally plausible variability, as inferred from the bootstrap-based uncertainty analysis.

For each sampled parameter set, the model was numerically integrated over a 24 h period under hypoxic conditions, and selected outputs were computed from the resulting trajectories. First-order Sobol indices were then computed for each parameter–output pair, measuring the fraction of the total output variance attributable to variation in that parameter alone.

By conditioning the sensitivity analysis on the calibrated parameter values and experimentally constrained parameter ranges, the resulting Sobol indices reflect physiologically relevant control mechanisms rather than purely structural properties of the model. This approach ensures that parameters identified as influential exert a quantitatively meaningful impact on model behaviour within the biological context considered.

To assess robustness beyond the physiologically constrained interval, an auxiliary wide-domain sensitivity scan was performed by varying the calibrated parameters within $[eqn]$ of their estimated values while preserving positivity. Although absolute variance magnitudes increased under this broader exploration, the ranking of first-order Sobol indices remained qualitatively preserved, with $[eqn]$ retaining dominant influence on ATP-related outputs. This indicates that the identified control hierarchy is not exclusively an artefact of narrow parameter bounds.

5.2. Sensitivity Outputs

Sensitivity was evaluated with respect to four complementary outputs:

• (i)Final ATP level at 24 h:

$$[eqn]$$

• (ii)Cumulative ATP availability:

$$[eqn]$$

• (iii)Cumulative ROS burden:

$$[eqn]$$

• (iv)Energetic–redox balance ratio:

$$[eqn]$$

Explicitly,

$$[eqn]$$

This multidimensional output set distinguishes the parameters predominantly affecting energy maintenance from those modulating oxidative burden or energetic–redox trade-offs.

Across all outputs, $[eqn]$ consistently exhibits the largest first-order Sobol index within the explored admissible domain, followed by $[eqn]$ . The hypoxia amplification parameter $[eqn]$ contributes negatively to ATP-related outputs, reflecting its role in ROS escalation.

The dominance of glutathione synthesis therefore persists within experimentally constrained domains, although it should be interpreted as conditional upon the reduced cysteine-centred architecture and bounded parameter space.

6. Numerical Results

6.1. Model Calibration Against Experimental Extracellular Fluxes

The calibrated model reproduces the experimentally inferred extracellular uptake and secretion fluxes across all experimental conditions, with relative deviations below 7% for all measured metabolites (Table 2). Experimental fluxes ( $[eqn]$ ) were computed from extracellular metabolite concentration changes over a 24 h period, while the mean model-predicted fluxes ( $[eqn]$ ) correspond to the time-averaged reaction rates obtained from numerical simulations of the model equations. The close agreement between $[eqn]$ and $[eqn]$ , including cysteine uptake rates under hypoxia and cysteine supplementation, provides quantitative validation of the parameter estimation procedure and supports the use of the calibrated model for subsequent dynamical and sensitivity analyses.

6.2. Dynamic Behaviour of the Cysteine-Centred Metabolic Network Under Hypoxia

All dynamic simulations reported below were performed using parameters calibrated against experimental extracellular flux data.

We first analysed the baseline dynamic behaviour of the calibrated cysteine-centred metabolic network under hypoxic conditions, using cysteine availability consistent with the experimental setting. Numerical integration of the model over a 24 h period reveals a consistent sequence of metabolic responses (Figure 2). Intracellular cysteine concentrations decline rapidly during the early phase of the simulation and remain low thereafter, reflecting sustained net consumption within the model structure.

As cysteine availability decreases, GSH levels decline continuously over time, indicating that the antioxidant capacity is constrained by limited precursor supply and that no compensatory mechanism is sufficient to stabilise the GSH concentrations under hypoxia. As a consequence, reactive oxygen species (ROS) accumulate monotonically throughout the simulation window. The concurrent decrease in GSH and increase in ROS indicate a progressive weakening of the redox buffering capacity rather than maintenance of redox homeostasis.

The ROS dynamics do not reach saturation or equilibrium within the simulated time window, suggesting that detoxification processes represented in the model are insufficient to offset hypoxia-associated ROS production. The system therefore evolves away from a redox steady state, displaying sustained oxidative pressure rather than dynamic balance.

Hydrogen sulfide (H_2_S) exhibits a transient increase at early and intermediate times, followed by stabilisation at lower concentrations. Within the model, this behaviour arises from the temporary allocation of cysteine flux toward H_2_S-producing reactions, combined with degradation and declining substrate availability at later times. The cellular energetic state, represented by ATP concentration, shows a gradual and continuous decline over time, reflecting the reduced efficiency of energy-generating processes under hypoxic conditions and the absence of energetic recovery within the simulated timeframe.

While these simulations characterise the intrinsic hypoxic dynamics of the calibrated network, they do not isolate the specific contribution of cysteine availability from other structural aspects of the model. To explicitly disentangle the causal role of cysteine, we next simulated the system under hypoxic conditions with and without extracellular cysteine supplementation, while keeping all intracellular parameters fixed (Figure 3).

Under cysteine-depleted conditions, the GSH levels decline progressively, resulting in sustained ROS accumulation and a marked reduction in ATP availability. In contrast, cysteine supplementation preserves the intracellular redox buffering capacity by sustaining the GSH levels, limits ROS accumulation, and maintains higher ATP levels throughout the simulation window (Figure 3). Notably, ATP stabilisation under hypoxia emerges from coordinated effects on both antioxidant defence and alternative energy-supporting pathways, including enhanced H_2_S availability. Together, these results demonstrate that cysteine acts as a protective metabolic input under hypoxic stress, rather than merely reflecting condition-dependent parameter shifts.

Although direct intracellular measurements are not used for calibration, the qualitative ordering of simulated intracellular trends is consistent with the experimentally reported directionality in ovarian cancer hypoxia models. In particular, studies have observed the depletion of intracellular GSH and accumulation of oxidative stress markers under cysteine limitation, as well as the partial restoration of redox balance upon cysteine supplementation. The present simulations reproduce these directional patterns (GSH ↓, ROS ↑ under depletion; attenuation under supplementation), thereby providing indirect qualitative validation of the inferred intracellular dynamics. No claim of quantitative intracellular matching is made.

6.3. Global Sensitivity Analysis Identifies Key Metabolic Control Points

To quantitatively identify the parameters governing cellular fitness under hypoxia, a variance-based global sensitivity analysis was performed using Monte Carlo sampling over the biologically admissible parameter ranges. The final ATP concentration at 24 h was used as a proxy for cellular fitness.

The resulting first-order Sobol indices are summarised in Table 3.

Small negative values obtained for $[eqn]$ and $[eqn]$ ( $[eqn]$ ) arise from numerical estimation noise and are therefore reported as zero.

First-order effects indicate comparable primary influence of $[eqn]$ (maximal cysteine-to-glutathione flux) and $[eqn]$ (cysteine uptake rate) on ATP-related outputs, with secondary contribution from $[eqn]$ (maximal cysteine-to-H_2_S flux). In contrast, parameters governing basal and hypoxia-amplified ROS production ( $[eqn]$ and $[eqn]$ ) exhibited negligible first-order variance contribution within the explored domain.

Total-order indices confirm that $[eqn]$ and $[eqn]$ remain the dominant contributors to output variance when parameter interactions are taken into account, while $[eqn]$ exerts a moderate influence, and the ROS-production parameters remain weakly contributory. These sensitivity magnitudes therefore reflect quantitatively relevant control within physiologically constrained parameter ranges, rather than unrestricted structural dominance.

Pairwise bootstrap correlations (Section 4.5) indicate moderate compensatory behaviour between $[eqn]$ and $[eqn]$ , consistent with the shared dependence of intracellular allocation on cysteine availability. Nevertheless, no near-singular dependencies were detected, and the global sensitivity ranking remains stable within the calibrated neighbourhood.

Structural Robustness to ATP Formulation

To evaluate the structural robustness of the identified control hierarchy, the simulations were repeated under an alternative ATP formulation in which oxidative phosphorylation efficiency was independent of the glutathione concentration (ATP production proportional only to $[eqn]$ ). Although the absolute ATP levels differed under this alternative structure, the relative ranking of Sobol indices remained qualitatively preserved, with $[eqn]$ and $[eqn]$ retaining primary influence. This indicates that the observed control hierarchy is not an artefact of the specific ATP representation adopted but rather a robust property within the imposed reduced cysteine-centred network architecture and bounded parameter domain.

6.4. Robustness of the Metabolic Adaptation Regime

Across the explored parameter space, the model displays robust qualitative behaviour, with no abrupt transitions or numerical instabilities.

The sensitivity hierarchy identified in Section 6.3 remains consistent across the admissible parameter domain, with the comparable primary influence of $[eqn]$ and $[eqn]$ and secondary contribution from $[eqn]$ . No parameter exhibited abrupt dominance shifts under moderate parameter variation.

Collectively, these results indicate that hypoxic adaptation in ovarian cancer cells is governed by coordinated cysteine uptake and allocation, rather than by a single isolated flux parameter. The robustness of this hierarchy supports the interpretation that metabolic control emerges from balanced cysteine partitioning within the reduced network architecture.

6.5. Regime Exploration Under Hypoxia

To investigate how oxygen availability shapes system-level metabolic behaviour, the hypoxia parameter H was systematically varied from normoxia ( $[eqn]$ ) to severe hypoxia ( $[eqn]$ ). For each value of H, steady-state levels of ATP and reactive oxygen species (ROS) were computed from the model.

As shown in Figure 4, this parameter scan reveals the emergence of three operationally distinct metabolic response regions. At low hypoxia, ATP levels remain largely stable while ROS levels are low, indicating a regime of sustained energetic capacity. At intermediate values of H, ATP exhibits a monotonic decline accompanied by a gradual accumulation of ROS, consistent with progressive energetic depletion. At high hypoxia, ATP levels collapse, and ROS increase sharply, defining a high oxidative burden regime.

These regions do not correspond to bifurcation points but rather to continuous yet steep variations in steady-state observables as H increases.

6.6. Linear Stability Analysis and Structural Neutrality

To formally assess whether the hypoxia-dependent transitions observed in Figure 4 correspond to true dynamical bifurcations, a linear stability analysis of the steady state was performed across the full hypoxia domain $[eqn]$ .

Let $[eqn]$ denote the state vector. The system (3)–(7) can be compactly written as $[eqn]$ .

For each value of H, the system was numerically integrated to convergence, and the Jacobian matrix

$$[eqn]$$

was evaluated at the computed steady state $[eqn]$ . Let $[eqn]$ denote the eigenvalues of $[eqn]$ . The dominant real part satisfies

$$[eqn]$$

indicating the absence of eigenvalues with a strictly positive real part and therefore the exclusion of local bifurcations within the explored physiological range. Representative numerical values across the hypoxia domain remained within floating-point precision of zero.

Interestingly, the eigenvalue spectrum consistently contains a structurally zero eigenvalue (Figure 5). This neutral direction arises from the bilinear redox interaction term $[eqn]$ , which imposes at steady state the constraint

$$[eqn]$$

thereby defining a one-dimensional invariant manifold within the $[eqn]$ subsystem. While this induces structural degeneracy in the redox variables, all transverse eigenvalues possess strictly negative real parts.

Consequently, the steady state is normally attractive: perturbations orthogonal to the neutral manifold decay exponentially, and no instability emerges across the hypoxia domain. The hypoxia-dependent regimes described in Section 6.5 therefore correspond to continuous quantitative transitions in steady-state levels rather than bifurcation-driven qualitative state changes.

This neutral direction reflects the minimal redox representation adopted in the reduced network and does not imply multistability or alternative stable states.

This result reinforces that the observed regime structure reflects graded metabolic adaptation within a globally stable reduced network architecture.

7. Discussion

This study develops a reduced flux-constrained mathematical framework to analyse cysteine-centred metabolic adaptation under hypoxic conditions in ovarian cancer cells. The model provides a quantitative systems-level explanation for the experimentally observed dependence of ovarian cancer cells on cysteine under hypoxia, a pro-oxidative stress condition [8] that characterises regions of the ovarian cancer microenvironment [9,10].

The model does not claim quantitative prediction of intracellular metabolite concentrations. Rather, intracellular dynamics are inferred as latent variables constrained by extracellular flux consistency and biologically supported qualitative behaviour reported in independent studies [6,7,7]. The framework therefore provides mechanistic plausibility under flux constraints rather than direct intracellular validation.

In ovarian cancer, hypoxia and oxidative stress coexist and are associated with the upregulation of regulatory factors such as hypoxia inducible factor–1 $[eqn]$ (HIF–1 $[eqn]$ ) and nuclear factor erythroid 2-related factor 2 (NRF2) [11]. These transcriptional programmes contribute to redox adaptation and may promote resistance to conventional chemotherapy [12,13]. Within this biological context, the present model isolates cysteine allocation as a determinant of redox–energetic coupling under experimentally imposed hypoxia.

It is essential to state explicitly the epistemic scope of the framework. Calibration relies exclusively on extracellular metabolite fluxes measured over a 24 h interval. Consequently, intracellular state variables (GSH, ROS, H_2_S, ATP) are not directly validated experimentally within this study. Rather, their trajectories represent dynamically inferred latent variables constrained by (i) extracellular flux consistency, (ii) biochemical plausibility, and (iii) qualitative agreement with the independent literature. The model therefore provides constrained mechanistic plausibility rather than quantitative intracellular validation.

Within experimentally supported parameter ranges, GSH synthesis capacity emerges as the dominant determinant of ATP stability. This observation is consistent with previous experimental studies demonstrating that cysteine-sustained GSH turnover is crucial for ovarian cancer cell adaptation to hypoxia and resistance to carboplatin cytotoxicity [6]. The model does not imply that GSH directly produces ATP; rather, it captures the requirement of adequate redox buffering for sustained oxidative phosphorylation efficiency. Thus, the apparent predominance of GSH synthesis reflects the central role of redox control under pro-oxidative hypoxic stress.

In the clinical setting, cysteine is the predominant thiol in peripheral blood serum and ascitic fluid of cancer patients [6], and both cysteine and homocysteine levels have been proposed as biomarkers distinguishing healthy women, women with benign ovarian tumours, and women with ovarian cancer [6]. The present modelling results therefore provide a mathematical formalisation of experimentally observed cysteine dependency. However, no direct quantitative comparison with patient-level metabolomic ranges is performed here; extrapolation to clinical stratification remains a hypothesis requiring independent validation.

The extracellular cysteine concentrations explored in the simulations lie within the same order of magnitude as values reported in serum and ascitic fluid metabolomic studies, supporting the physiological plausibility of the imposed input ranges. However, no direct quantitative patient-level calibration is performed, and translational inference remains hypothetical pending dedicated clinical validation.

Beyond GSH synthesis, the model quantifies the contribution of cysteine-derived H_2_S. The results underscore the metabolic redundancy and flexibility in cysteine allocation. In the reduced network considered, H_2_S exerts a secondary stabilising influence on ATP levels. This is consistent with the reported roles of H_2_S as both an antioxidant molecule contributing to redox regulation [14] and an electron donor to the mitochondrial electron transport chain [15]. Importantly, in the present framework, H_2_S contribution is represented phenomenologically rather than via explicit modelling of electron transport kinetics.

ATP levels in the simulations remain above critical thresholds in parameter regimes where cysteine allocation toward both GSH and H_2_S is preserved. Within the confines of the reduced network, this supports the hypothesis that cysteine metabolism contributes to ovarian cancer cell fitness under oxidative stress conditions such as hypoxia and drug exposure [16]. However, this conclusion is conditional upon the imposed network architecture and bounded parameter ranges.

A critical methodological point concerns metabolic simplification. The model intentionally excludes glycolysis, fatty acid oxidation, NADPH-producing reactions, and one-carbon metabolism. These pathways undeniably influence ATP production and redox control. Their omission should not be interpreted as biological irrelevance but as a deliberate reduction strategy to avoid underdetermined parameter spaces and uncontrolled compensatory flux redistribution. Consequently, the identified control hierarchy is conditional on this reduced cysteine-centred subnetwork.

Although no bifurcation structure was detected, linear stability analysis confirmed that hypoxia-induced regime transitions correspond to continuous steady-state deformations within a globally stable reduced network.

Finally, the model predicts that the simultaneous disruption of cysteine uptake and downstream allocation pathways may destabilise the redox–energetic balance under oxidative stress. Within the restricted scope of the current framework, this suggests a potential systems-level rationale for targeting cysteine metabolism in combination therapeutic strategies aimed at overcoming hypoxia-associated chemoresistance. Nevertheless, predictive use in therapeutic planning would require integration with intracellular validation, expanded metabolic representation, and patient-derived quantitative data.

8. Limitations

While the proposed model captures key features of cysteine-centred metabolic adaptation to hypoxia, its scope and inferential limits must be clearly delineated.

First, the model calibration relies exclusively on extracellular metabolite fluxes derived from exometabolomic measurements over a finite experimental interval. Consequently, the model cannot identify intracellular metabolite pool sizes, transport kinetics, or enzyme-level regulatory parameters from exometabolome data alone. Intracellular states are inferred as dynamically consistent latent variables constrained by extracellular flux balance and biochemical plausibility rather than as independently validated quantitative measurements. Accordingly, the conclusions of the model are restricted to flux-consistent mechanistic inference under controlled in vitro conditions.

Second, exchange processes between extracellular and intracellular compartments are represented implicitly rather than through explicit transport equations. As a result, microenvironmental gradients, transporter saturation effects, and dynamic exchange kinetics are not resolved.

Third, oxygen availability is represented by a single abstract hypoxia parameter H, rather than by explicit modelling of oxygen transport, diffusion, or consumption. The framework therefore does not capture spatial or temporal heterogeneity in oxygen tension, nor feedback coupling between mitochondrial respiration and oxygen depletion.

Fourth, although the experimental data originate from large populations of ovarian cancer cells, including two distinct cell lines, the model describes an average intracellular metabolic state. It does not account for cell-to-cell heterogeneity, stochastic variation, proliferation dynamics, or cell death processes.

The metabolic structure itself is intentionally restricted to a minimal cysteine-centred subnetwork. Core pathways known to contribute to hypoxic adaptation, including glycolysis, fatty acid oxidation, NADPH-generating reactions, and the methionine one-carbon cycle, are not explicitly represented. This reductionist architecture was adopted to avoid underdetermined parameter spaces and uncontrolled compensatory flux redistribution. However, it necessarily limits the metabolic completeness of the framework. Incorporation of these pathways in future model extensions will be required to evaluate whether the identified control hierarchy remains stable under broader metabolic coupling conditions.

From an identifiability perspective, explicit incorporation of glycolysis, NADPH-generating reactions, fatty acid oxidation, and one-carbon metabolism would introduce a substantial increase in dimensionality, adding more than 15–20 additional kinetic parameters in minimal representations. Given calibration based solely on extracellular fluxes over a finite time interval, such expansion would lead to severely underdetermined parameter spaces and amplified compensatory flux redistribution. The present reduction strategy therefore reflects a controlled trade-off between the biological scope and inverse-problem well-posedness.

9. Conclusions

We present a reduced mechanistic mathematical framework linking cysteine metabolism to hypoxic adaptation in ovarian cancer cells. By calibrating the model exclusively against experimentally derived extracellular metabolite fluxes and constraining the parameter uncertainty within biologically admissible ranges, we establish a flux-consistent systems-level description of cysteine allocation between GSH synthesis and H_2_S production under hypoxic stress.

Within the imposed network structure and calibrated parameter domain, GSH synthesis emerges as the principal determinant of ATP maintenance, while H_2_S production provides a secondary stabilising contribution. These findings should be interpreted as conditional on the reduced cysteine-centred subnetwork considered and do not imply exhaustive metabolic control dominance.

The framework does not uniquely identify intracellular metabolite pool sizes, enzyme-level kinetics, or transport dynamics; rather, it defines dynamically feasible intracellular behaviours consistent with extracellular constraints and established biochemical mechanisms. Its primary contribution lies in formalising experimentally observed cysteine dependency into a mathematically disciplined inference structure.

Future developments should incorporate additional metabolic pathways, including glycolysis, NADPH-generating reactions, fatty acid oxidation, and one-carbon metabolism, as well as explicit oxygen dynamics and population-level effects, in order to evaluate compensatory flux redistribution and broader metabolic coupling. Integration with intracellular measurements and quantitative patient-derived metabolomic data, including peripheral blood serum and ascitic fluid profiling, will be required to assess the predictive validity and translational relevance.

In its present form, the model provides a constrained systems-level platform for hypothesis generation and controlled expansion toward more comprehensive representations of hypoxic cancer metabolism.