Clearance Asymmetry, Not Aquaporin Abundance, Sets Plasma-Activated-Water Selectivity for Glioblastoma: A Transfer-Function Dissection

1. Introduction

Glioblastoma (GBM), the most common and most aggressive primary malignant brain tumour in adults, remains among the least tractable of human cancers. Despite maximal safe surgical resection followed by concomitant and adjuvant temozolomide chemoradiotherapy—the standard of care established two decades ago [1] and refined chiefly through molecular stratification such as MGMT promoter methylation status [2]—median overall survival remains on the order of 15 months, and the five-year relative survival rate for malignant brain and other central nervous system tumors is approximately 35% [3]. Glioblastoma alone accounts for roughly half of all malignant primary CNS tumors [3]. The diffuse, infiltrative growth of the tumor, the protection afforded by the blood–brain barrier, and the profound intratumoral heterogeneity that underlies therapeutic resistance together place a hard ceiling on the efficacy of conventional cytotoxic and radiotherapeutic approaches. This persistent clinical impasse has motivated sustained interest in physically based and redox-based modalities that might exploit vulnerabilities orthogonal to those targeted by alkylating chemotherapy.

Among these emerging modalities, cold atmospheric plasma (CAP) and its liquid derivatives—plasma-activated water (PAW) and plasma-activated medium (PAM) have attracted particular attention. A partially ionized gas generated at near-ambient temperature, CAP delivers to biological targets a complex and tunable cocktail of reactive oxygen and nitrogen species (RONS), together with charged particles, photons, and electric fields [4,5,6].The defining and most therapeutically tantalizing observation in plasma oncology is one of selectivity: under appropriately chosen conditions, plasma exposure induces death in malignant cells while leaving their non-transformed counterparts comparatively unharmed [7]. This selectivity, first articulated as a possible paradigm shift in cancer therapy [7], has since been reported across a wide range of tumor types and experimental configurations, including ovarian [8] and metastatic breast [9] carcinoma models, and has been reviewed extensively as a defining feature of the field [10]. The same selective behavior persists when plasma is not applied directly but instead through a liquid intermediary, indicating that long-lived species carried in solution—rather than transient plasma-phase agents alone—are sufficient to reproduce the differential cytotoxicity [11].

For glioblastoma specifically, the selectivity phenomenon has been documented in directly comparative settings. Plasma-activated medium has been shown to kill GBM cells preferentially relative to normal astrocytes [11], and direct CAP exposure produces several-fold greater death in U87-MG glioblastoma cells than in matched non-malignant astrocytes under identical treatment [12]. Complementary work has demonstrated that plasma treatment can restore or enhance the sensitivity of otherwise resistant tumor cells [13]. The U87-MG line, despite well-recognized caveats regarding the fidelity of immortalized cell lines to primary tumor tissue, has accordingly become a standard in vitro model for mechanistic studies of plasma–glioblastoma interactions. The reproducibility of GBM-selective killing across laboratories and plasma sources establishes the phenomenon as robust; what remains unsettled, and what motivates the present work, is its mechanistic origin.

The chemistry underlying these effects is comparatively well characterized at the level of bulk composition. Aqueous-phase plasma chemistry generates, among other species, hydrogen peroxide (H₂O₂), nitrite and nitrate, and peroxynitrite, with the longer-lived members of this set surviving long enough to mediate biological effects far from the point of generation [14]. Quantitative measurements in plasma-treated buffers place H₂O₂ and nitrite concentrations in the micromolar-to-low-millimolar range under typical treatment protocols [15], and the generation, transport, and biological interactions of these species have been comprehensively reviewed [16]. Of the durable species, hydrogen peroxide occupies a privileged mechanistic position. It is not merely a cytotoxin but a bona fide signaling molecule whose biological action is governed by precise, compartmentalized chemistry rather than indiscriminate oxidation [17]. The recognition that the same oxidant can drive proliferative signaling at low fluxes and apoptosis at high fluxes—and that its effects depend critically on where and how rapidly it accumulates within the cell—reframes the selectivity question. If plasma selectivity is ultimately a question of intracellular H₂O₂ dose, then the determinants of that dose, namely the rates at which H₂O₂ enters and is cleared from each cell type, become the proper objects of mechanistic inquiry [4].

Three broad and non-exclusive hypotheses have been advanced to explain why malignant cells accumulate a more lethal RONS burden than normal cells under equivalent exposure. The first locates selectivity in the plasma membrane itself. Molecular dynamics studies have shown that the lipid peroxidation products generated by RONS reduce lipid order and, in the limit of extensive oxidation, can drive transient pore formation, thereby increasing membrane permeability to further species [18,19,20]. Because the structural and dynamic consequences of oxidation depend on the underlying lipid composition, differences in membrane makeup between cancer and normal cells have been proposed as a basis for differential permeabilization [21,22] Cholesterol content is a particularly important modulator: simulations indicate that cholesterol raises the free-energy barrier for RONS permeation and stabilizes the bilayer against oxidative injury, leading to the explicit proposal that the relative cholesterol content of cancer versus normal membranes could underlie plasma selectivity [23,24,25]. The distribution and mobility of reactive species within bilayers, which set the kinetic stage for permeation, have likewise been characterized in detail [26]. This membrane hypothesis is attractive because it requires no specialized machinery, but it carries a specific and testable corollary: the malignant membrane must be intrinsically more permeable, or more readily permeabilized, than the normal one.

The second hypothesis assigns the decisive role to facilitated transport through aquaporin water channels. A subset of aquaporins, termed peroxiporins, conduct hydrogen peroxide across the plasma membrane in addition to water. The foundational demonstration that specific aquaporins facilitate transmembrane H₂O₂ diffusion identified the unorthodox channel AQP8—but not the classical AQP1—as a conduit for peroxide [27], a pattern subsequently confirmed in mammalian cells, where AQP3 and AQP8 promote H₂O₂ uptake and downstream redox signaling [28]. Aquaporin-mediated peroxide transport has since been shown to be required for physiological signaling events in several systems [29], and additional plant and mammalian peroxiporins have been characterized [30]. Critically, peroxiporin conductance is not static: H₂O₂ permeation through AQP8 is gated by the redox state of conserved cysteine residues, providing a mechanism by which the channel can be reversibly closed under oxidizing conditions [31]. The biophysical permeability of individual channels to peroxide has been quantified for several isoforms [32]. This hypothesis acquires oncological relevance from two observations: aquaporins are frequently over-expressed in tumors, where they contribute to growth, migration, and angiogenesis [33], and AQP8-linked redox signaling has been implicated directly in malignant cells [34]. If selectivity is channel-mediated, then the abundance, isoform identity, and gating state of peroxiporins in each cell type should govern the differential influx.

The third hypothesis shifts the emphasis from influx to the cell’s capacity to detoxify and clear incoming oxidant and is most fully developed in the work of Bauer and colleagues. In this account, malignant cells are protected from RONS-driven apoptosis by a membrane-associated catalase that decomposes extracellular H₂O₂ and intercepts peroxynitrite-derived signaling [35]. Plasma-derived species transiently inactivate a fraction of this protective catalase; once the enzymatic shield is weakened, tumor-cell-derived H₂O₂ enters the cytoplasm—explicitly via aquaporins—and depletes intracellular glutathione, establishing a self-amplifying, apoptosis-inducing RONS signaling state that does not arise in non-malignant cells [36]. The cooperative interaction between catalase inactivation, aquaporin/H₂O₂-mediated glutathione depletion, and the countervailing glutathione system has been dissected in detail, leading to the proposal that tumor cells bearing more aquaporins are positively selected during oncogenesis precisely because they admit the H₂O₂ flux required for proliferative signaling [37]. This framework is notable because it unifies all three candidate elements—membrane composition, aquaporin transport, and clearance capacity—into a single causal chain, and because it makes the clearance asymmetry between cell types, rather than any single transport step, the proximate determinant of selective killing.

These three hypotheses are not mutually exclusive, and indeed the Bauer framework explicitly couples them; yet the field lacks an integrated quantitative treatment capable of weighing their relative contributions to the selectivity that is actually observed. The membrane hypothesis rests in part on assumed compositional differences between cancer and normal cells, which for glioblastoma are only partially characterized. Lipidomic and mass-spectrometry-imaging studies consistently report cholesterol enrichment in glioblastoma tissue relative to normal brain [38,39] and reduced sphingomyelin in tumor-cell regions [40], but the phospholipid and saturation profile of GBM is heterogeneous across molecular subtype, tumor region, and stem-cell state, complicating any attempt to specify a single representative malignant membrane against a normal reference such as the canonical mammalian plasma-membrane composition [41]. The aquaporin hypothesis, for its part, is hampered for AQP8 by the absence of an experimentally determined channel structure, which forces structural reasoning to rely on homology and computational models whose selectivity filters may not faithfully reproduce peroxide conductance [42]. And while the clearance hypothesis is mechanistically the most complete, it has not been cast in a form that predicts, quantitatively, how a given asymmetry in clearance and influx rates translates into a fold-difference in cytotoxic outcome.

The central question, then, is not whether each of these mechanisms can operate in principle—each is supported by direct evidence—but which of them, and in what combination, governs the magnitude of the selectivity observed in glioblastoma, and how robustly. Answering this requires a framework that connects the upstream plasma chemistry to a downstream cytotoxic phenotype through the intervening transport and clearance steps, with each step parameterized from measurable quantities, so that the sensitivity of the predicted selectivity to each parameter can be assessed directly. In the present work we construct such a framework. We treat the determination of intracellular H₂O₂ dose as a transport problem, deriving a compact steady-state expression in which the cytoplasmic concentration is set by the balance of channel- and lipid-mediated influx against first-order clearance, and we use this expression to define a dose-based selectivity factor between malignant and normal cells. We embed within this treatment the redox gating of AQP8 and a survival model that maps cumulative cytoplasmic dose onto cell death. We ground the membrane component of the model in all-atom molecular dynamics simulations of compositionally distinct glioblastoma-like and astrocyte-like bilayers, from which we extract the passive permeability-relevant structural and dynamic properties directly. Finally, we subject the resulting model to global variance-based sensitivity analysis and Monte Carlo robustness testing, in order to identify which biophysical parameters actually control the predicted selectivity and which do not, and to determine whether the selective regime is a fragile or a robust feature of the parameter space.

This integrated treatment yields several contributions. First, it provides a transparent transfer-function decomposition of plasma selectivity into a small number of physically interpretable quantities. Second, it localizes, through three independent and mutually consistent computational results, the origin of selectivity away from the passive lipid membrane. Third, it identifies clearance asymmetry and the operating-point sensitivity of channel permeability—rather than the fold-change in aquaporin abundance between cell types—as the dominant levers of the predicted selectivity, a finding with direct implications for how one would seek to enhance or defeat it. Fourth, it shows that the membrane’s documented cholesterol enrichment in glioblastoma, far from favoring tumor uptake, acts in the direction of reduced rather than increased permeability, reinforcing the displacement of selectivity onto regulated transport and clearance. Fifth, it casts the resulting picture as a set of falsifiable predictions and a concrete experimental programmed, so that the model can be tested and, where necessary, refuted. Throughout, we are careful to distinguish what is grounded in direct simulation or measurement from what is assumed, and to state the conditions under which each conclusion would fail.

2. Methods

2.1. Molecular Dynamics of Glioblastoma-Like and Astrocyte-Like Membranes

To probe whether passive lipid-phase transport could account for the differential RONS sensitivity of malignant and normal cells, we constructed two compositionally distinct atomistic bilayers intended to represent, respectively, a glioblastoma-like (U87-MG) plasma membrane and a non-malignant astrocyte plasma membrane. Both systems were assembled with the CHARMM-GUI Membrane Builder [43,44,45], each containing 512 lipids arranged with physiological transbilayer asymmetry—the choline phospholipids and sphingomyelin in the outer leaflet, and phosphatidylethanolamine, phosphatidylserine, and phosphatidylinositol in the inner leaflet, with cholesterol distributed between the two—together with explicit solvent and neutralizing ions. The lipid inventories comprised cholesterol together with the principal glycerophospholipid and sphingolipid classes of the mammalian plasma membrane: phosphatidylcholine, phosphatidylethanolamine, phosphatidylinositol, phosphatidylserine, and sphingomyelin, in both di-unsaturated and mono-unsaturated acyl-chain variants. The base composition was anchored to the experimentally informed average mammalian plasma-membrane model of Ingólfsson et al. (2014) [41], which we take as the non-malignant reference. The malignant composition differed from the reference in a small number of deliberate respects, chosen to reflect documented features of the glioblastoma lipidome rather than measured directly from the specific cell line. The U87-MG-like membrane was assigned a modestly higher cholesterol fraction (approximately 42% versus 40% of total lipid), a shift within the phosphatidylcholine and phosphatidylethanolamine pools from di-unsaturated toward mono-unsaturated species, and a small reduction in sphingomyelin content. These choices were grounded in the consistent reports of cholesterol enrichment in glioblastoma tissue relative to normal brain [38,39] and of reduced sphingomyelin in tumor-cell regions [40]. We emphasize that the glioblastoma lipidome is heterogeneous across molecular subtype, tumor region, and stem-cell state [46,47] and that immortalized cell lines may diverge from primary tissue; accordingly, these two bilayers are best understood as one plausible glioblastoma-like membrane set against a canonical normal reference, with the cholesterol difference—the most robustly documented feature, and the one that dominates the biophysical comparison—being the principal point of distinction. The full per-system lipid counts are reported in Table S1.

Interactions were described by the CHARMM36 all-atom additive force field for lipids [48], with water represented by the TIP3P model [49] and standard CHARMM ion parameters. All simulations were performed in GROMACS [50]. Each system was energy-minimized by steepest descent and then equilibrated through the standard staged CHARMM-GUI protocol, in which positional and dihedral restraints on the lipids are progressively released, before unrestrained production. Production dynamics were integrated with a 2 fs timestep; bonds to hydrogen were constrained with the LINCS algorithm [51]. Long-range electrostatics were treated by the particle-mesh Ewald method [52,53] with Lennard-Jones and real-space Coulomb interactions truncated at a short-range cutoff with a force-switching function. The temperature was maintained at 310.15 K with a stochastic velocity-rescaling thermostat [54] and the pressure at 1 bar with a semi-isotropic stochastic cell-rescaling barostat [55], so that the membrane normal and lateral dimensions could fluctuate independently. For each of the two compositions, three independent replicates were propagated for 200 ns, yielding an aggregate of 1.2 µs of trajectory across the six runs. Because the structural and dynamic observables analyzed here are mass-independent, the use of a 2 fs rather than a hydrogen-mass-repartitioned 4 fs timestep does not affect the reported quantities. Trajectories were analyzed with MDAnalysis [56,57] using the CHARMM-GUI-generated protein-structure files as topology. Four observables relevant to passive permeability were extracted over an 80–200 ns analysis window, discarding the first 80 ns as equilibration. Acyl-chain ordering was quantified through the carbon–deuterium order parameter, −S_CD, computed per methylene segment from the orientation of the reconstructed C–H bond vectors relative to the bilayer normal and averaged over chains and frames for the glycerophospholipid species (sphingomyelin and cholesterol excluded): $-S_{\mathrm{CD}}=-{\displaystyle \frac{1}{2}}⟨3{\cos }^2\theta -1⟩$ where θ is the angle between the instantaneous C–H bond vector and the membrane normal z. The area per lipid was obtained by partitioning the lateral box area among the phosphate-bearing lipids; because cholesterol occupies interstitial volume without a headgroup phosphate, the raw area-per-phospholipid is confounded by the differing cholesterol fractions of the two systems, and we therefore additionally report a cholesterol-corrected area normalized to the total lipid count, defined as: $⟨A_{\mathrm{total}}⟩={\displaystyle \frac{2⟨A_{\mathit{xy}}⟩}{N_{\mathrm{total}}}}$ where A_xy is the instantaneous lateral cross-sectional area of the simulation box and N_total = 512 is the total number of lipids (phospholipids plus cholesterol) across both leaflets. Bilayer thickness was measured as the mean inter-leaflet phosphate-to-phosphate distance (d_P–P). Lateral mobility was characterized by the in-plane lateral diffusion coefficient, D_lat, obtained from the slope of the centre-of-mass-removed lateral mean-squared displacement of the lipids after unwrapping periodic boundaries, fitted over an intermediate lag-time interval (20–80 ns): MSD_xy(Δt) = ⟨|r_xy(t+Δt) − r_xy(t)| ²⟩. Because the lateral mean-squared displacement remains sub-diffusive on the accessible 200 ns timescale (scaling exponent α≈0.6 in MSD∝Δt^α), D_lat is reported as an effective lateral mobility used for the malignant-versus-normal comparison rather than a fully converged Einstein diffusivity. Differences between the malignant and normal systems were assessed across replicates by Welch’s unequal-variances t-test, with the three replicates per composition treated as the units of replication.

2.2. Transfer-Function Model of Intracellular Dose

To connect the extracellular plasma chemistry to a cytotoxic outcome, we treat the determination of intracellular hydrogen-peroxide concentration as a membrane-transport problem and derive a compact steady-state relation. We consider a single cell of surface area A and volume V exposed to a constant extracellular concentration C_out of the principal long-lived oxidant, taken here to be H₂O₂, and we follow the cytoplasmic concentration C_in. The total membrane permeability to peroxide is decomposed into a channel-mediated and a lipid-mediated contribution,

P = P_AQP8 ⋅ g + P_lip, where P_AQP8 = ρ_AQP8 ⋅ p_f

(1)

where p_f is the single-channel unitary permeability (cm³ s⁻¹) and ρ_AQP8 is the functional channel surface density (cm⁻²), yielding the macroscopic channel permeability P_AQP8 (cm s⁻¹), g ∈ (0,1] is a dimensionless gating factor (Section 2.3), and P_lip is the background permeability of the protein-free lipid bilayer. The channel term embodies the facilitated-diffusion hypothesis, in which AQP8—but not classical aquaporins such as AQP1—conducts H₂O₂ across the membrane [27,28]; the single-channel permeability is taken from biophysical determinations of peroxiporin conductance [32]. The lipid term embodies the passive-membrane hypothesis; its baseline magnitude is taken from direct experimental measurement of H₂O₂ permeability across protein-free, cholesterol-containing phospholipid membranes [58], while its dependence on membrane composition—in particular the suppression of peroxide permeation by cholesterol—is informed by molecular-dynamics free-energy studies [24,26].

Mass balance across the membrane, with first-order intracellular clearance at rate k_clear, gives:

$${\displaystyle \frac{d{C}_{\mathrm{in}}}{\mathit{dt}}}=\left({\displaystyle \frac{\mathit{PA}}{V}}\right)\left(C_{\mathrm{out}}-C_{\mathrm{in}}\right)-k_{\mathrm{clear}}\cdot C_{\mathrm{in}}$$

(2)

Defining the effective membrane transport rate constant as $k_m={\displaystyle \frac{\mathit{PA}}{V}}={\displaystyle \frac{\left(P_{\mathrm{AQP}8}\cdot g+P_{\mathrm{lip}}\right)A}{V}}$ (with units of s⁻¹), Equation (2) simplifies to ${\displaystyle \frac{d{C}_{\mathrm{in}}}{\mathit{dt}}}=k_m\left(C_{\mathrm{out}}-C_{\mathrm{in}}\right)-k_{\mathrm{clear}}{C}_{\mathrm{in}}$, and the transport number can be compactly expressed as $\mathsf{\Phi }={\displaystyle \frac{k_m}{k_{\mathrm{clear}}}}$.

At steady state (dC_in/dt = 0), this yields the central result:

$$C_{\mathrm{in},\mathrm{ss}}=C_{\mathrm{out}}\cdot {\displaystyle \frac{\mathsf{\Phi }}{1+\mathsf{\Phi }}}$$

(3)

with the dimensionless transport number:

$$\mathsf{\Phi }={\displaystyle \frac{\left(P_{\mathrm{AQP}8}\cdot g+P_{\mathrm{lip}}\right)\left(A/V\right)}{k_{\mathrm{clear}}}}={\displaystyle \frac{\left({\rho }_{\mathrm{AQP}8}\cdot p_f\cdot g+P_{\mathrm{lip}}\right)\left(A/V\right)}{k_{\mathrm{clear}}}}$$

(4)

Equation (4) compresses the transport problem into a single governing quantity: when Φ≪1 the cell is clearance-dominated and admits only a small fraction of the extracellular concentration, whereas when Φ≫1 the cell is influx-dominated and the cytoplasmic concentration approaches the extracellular value. The selectivity between two cell types follows from the ratio of their steady-state internal concentrations, and a regime in which one cell type lies below the transition (Φ < 1) while the other lies above it (Φ > 1) constitutes a switch-like selective window. The clearance term k_clear carries the contribution of the protective antioxidant machinery, including the membrane-associated catalase whose transient inactivation by plasma-derived species is central to the redox account of selectivity [35,36,37]; a lower effective clearance in malignant cells therefore enters the model as a reduced k_clear. The extracellular driving concentration C_out is set to a value representative of plasma-activated buffers [14,15].

2.3. Redox Gating of AQP8

The hydrogen peroxide conductance of AQP8 is not constitutive; rather, it is dynamically regulated by the redox state of conserved cysteine residues (specifically Cys53 in human AQP8, UniProt O94778), which undergo reversible oxidation under severe oxidative stress to restrict pore conduction [31]. In the steady-state and transfer-function formulations of our transport model (Equations 1–4), this negative-feedback behavior is represented phenomenologically via a cooperative Hill gating factor that decreases as the cytoplasmic oxidant concentration rises:

$$g\left(C_{\mathit{in}}\right)={\displaystyle \frac{1}{1+{\left(\frac{C_{\mathit{in}}}{C_{\mathit{half}}}\right)}^n}}$$

(5)

where C_half = 40 μM denotes the intracellular concentration yielding 50% channel closure, and n = 2 is the Hill coefficient governing the cooperativity of the structural transition [31]. Because n > 1, the gate introduces non-linear negative feedback that caps channel-mediated influx under high exposure. Crucially, this brake operates predominantly in the influx-dominated malignant cell (where C_in > C_half) while leaving the clearance-dominated astrocyte virtually unattenuated (C_in≪C_half), thereby partially compressing rather than amplifying the raw transport selectivity. To evaluate the physical validity of treating this covalent modification as an instantaneous equilibrium function of C_in, we formulated the explicit dynamic reaction-kinetic equation governing the open channel fraction g(t):

$${\displaystyle \frac{\mathit{dg}}{\mathit{dt}}}=-k_{\mathit{ox}}{C}_{\mathit{in}}\left(t\right)\hspace{0.17em}g\left(t\right)+k_{\mathit{red}}\left(1-g\left(t\right)\right)$$

(6)

where k_ox is the second-order rate constant for H₂O₂-mediated sulfenylation of the critical thiol (−SH + H₂O₂ → −SOH + H₂O), and k_red represents the cellular reduction rate mediated by the thioredoxin/glutaredoxin systems. In mammalian peroxiporins, typical thiol oxidation rate constants range between k_ox∼10¹ − 10² M⁻¹s⁻¹ [31]. At the micromolar operating concentrations reached inside glioblastoma cells (C_in≈50 − 70 μM), the characteristic relaxation time for channel closure is given by: ${\tau }_{\mathit{gate}}\approx {\displaystyle \frac{1}{k_{\mathit{ox}}{C}_{\mathit{in}}+k_{\mathit{red}}}}\sim {10}^1-{10}^2\mathrm{s}$. By comparison, the characteristic timescale for membrane transport and cytoplasmic concentration equilibration is governed by the cellular dilution and clearance rate: ${\tau }_{\mathit{transport}}={\left(P{\displaystyle \frac{A}{V}}+k_{\mathit{clear}}\right)}^{-1}\approx {\left(10+5\right)}^{-1}\mathrm{s}\approx 0.07\mathrm{s}$. Because τ_transport≪τ_gate, a pronounced separation of timescales exists: upon sudden exposure to plasma-activated water, the membrane reaches its unattenuated transport capacity within fractions of a second, whereas covalent oxidative gating requires tens of seconds to reach its steady-state closed conformation. Consequently, during the initial 10 − 60 s of treatment, AQP8 conducts in an uninhibited state (g≈1), admitting a transient, high-amplitude oxidant pulse into the malignant cell before the Cys53 brake engages.

The quasi-steady-state formulation of Equation (5)—which assumes instantaneous closure to g≈0.39—therefore underestimates the early cumulative influx and represents a strictly conservative lower bound on the intracellular dose delivered to glioblastoma cells.

2.4. Dynamics, Cumulative Dose, and Survival

Although the steady-state relation (3) captures the operating point, cytotoxicity depends on the integrated exposure. Integrating Equation (2) for a step exposure gives the approach to steady state,

$$C_{\mathrm{in}}\left(t\right)=C_{\mathrm{in},\mathrm{ss}}\left[1-e^{-\left({\displaystyle \frac{\mathit{PA}}{V}}+k_{\mathrm{clear}}\right)t}\right]$$

(7)

from which the cumulative cytoplasmic dose is obtained as the time integral:

$$D\left(t\right)={\int }_0^t{C}_{\mathrm{in}}\left(t\mathrm{\text{'}}\right)\hspace{0.17em}\mathit{dt}\mathrm{\text{'}}$$

(8)

We map dose onto outcome through a cooperative (Hill-type) survival relation:

$$S={\displaystyle \frac{1}{1+{\left(\frac{D}{D_{50}}\right)}^h}}$$

(9)

in which D₅₀ is the dose giving 50% survival and h is the Hill coefficient setting the steepness of the threshold; both parameters are fixed by least-squares calibration to a published glioblastoma kill curve (Section 2.5) rather than assumed. A single-hit exponential form,

$$S=\exp \left(-{\displaystyle \frac{D}{D_{\mathrm{lethal}}}}\right)$$

(10)

was also fitted but is rejected against the cooperative law by information criteria (ΔAICc≈+12), consistent with the sigmoidal, threshold-like dose-response observed experimentally; the cooperative form is not intended to resolve the detailed apoptotic cascade. Within the cooperative family, a Weibull form fitted the same calibration data marginally better than the Hill form (ΔAICc≈2), too small a gap for to be decisively preferred; the Hill form is retained for the explicit cooperativity interpretation of its shape parameter. The qualitative dose–survival behavior—a saturating, threshold-like law that converts a modest dose difference into a large survival difference—is common to both cooperative forms, whereas the specific window boundaries reported below are conditioned on the Hill calibration. The therapeutic window is then defined operationally as the range of exposure times and extracellular concentrations over which the normal-cell survival remains high while the malignant-cell survival falls below a chosen lethality threshold.

2.5. Parameterization

The baseline parameter values, their plausible ranges, and the empirical or biophysical basis for each are collected in Table 1. The values are intended to be order-of-magnitude representative rather than cell-line-specific measurements, and the sensitivity and robustness analyses of Section 2.6 and Section 2.7 are designed precisely to test how sensitive the model’s conclusions are to this choice. Briefly, the extracellular driving concentration was set to C_out = 100 μM, representative of durable peroxide levels measured in plasma-activated liquids [14,15]; the surface-to-volume ratio was set to A/V = 4000 cm⁻¹, consistent with a spherical or ellipsoidal mammalian cell of radius r≈7.5 μm; and the background lipid permeability was set to P_lip = 5 × 10⁻⁴ cm s⁻¹, grounded in direct biophysical measurements of protein-free, cholesterol-containing phospholipid bilayers [58]. The aquaporin contribution was parameterized through an effective macroscopic channel permeability expressed as a multiple of P_lip, taken as 1.0 × P_lip for the astrocyte and 4.0 × P_lip for the U87-MG membrane, the latter reflecting the multi-fold peroxiporin over-expression documented across astrocytoma grades [59,60,61,62].

The operational clearance rate constants were set to k_clear = 50 s⁻¹ for the astrocyte and 5 s⁻¹ for the malignant cell. While basal endogenous clearance capacities in resting non-stressed cells exhibit modest fold-differences (∼1.8× on average; Doskey et al., 2016) [63], this 10-fold operational asymmetry encodes the tumor-specific collapse of protective, membrane-associated catalase induced by plasma-derived species under the Bauer paradigm [35,36,37], which deprives the malignant cell of its primary enzymatic shield while leaving non-malignant clearance intact. Clearance is treated as first order in intracellular H₂O₂: catalase, the principal scavenger under the sustained, high peroxide flux modelled here, has a Michaelis constant (K_m) for H₂O₂ in the tens-to-hundreds of millimolar range, several orders of magnitude above the intracellular concentrations reached at the operating point (of order 10 μM), so removal remains strictly in its linear regime; a saturable, Michaelis–Menten treatment leaves the selectivity and normal-cell sparing unchanged for any Michaelis constant above roughly 10 μM (Figure S1).

The gating parameters were C_half = 40 μM and n = 2, capturing reversible oxidative closure of AQP8 at conserved cysteine residues [31]. The survival parameters were obtained by calibrating the cooperative Hill relation of Section 2.4 to the empirical hydrogen-peroxide kill curve of Kurake et al. (2016) [64] (D₅₀ = 3.06 × 10⁵ μM ⋅ s, h = 6.82, and R² = 0.98). Although this calibration curve was determined in U251SP glioblastoma cells, it provides the most rigorously quantified, concentration-dependent survival benchmark for H₂O₂ in high-grade malignant glioma available in the literature; it is employed here as a representative cytotoxic response function for malignant glioma rather than an idiosyncratic U87-MG parameter.

Because matched concentration-dependent survival curves for primary human astrocytes under identical pure H₂O₂ bolus exposures were not established in that dataset, the non-malignant reference threshold from Kurake et al. (2016) [64] (MCF10A, ≥93% survival at 64 μM) was utilized to benchmark the high-survival boundary. We treat this empirical datum strictly as a phenotypic proof-of-concept anchor demonstrating normal non-malignant cell tolerance at glioma-lethal peroxide exposures under identical assay conditions, whereas the astrocyte survival predicted by our model emerges independently from its physiological clearance capacity (k_clear = 50 s⁻¹) operating below the lethal threshold.

Finally, we deliberately separate two conceptually distinct quantities that are easily conflated: the absolute channel permeability of a given cell, which dictates its operating point on the saturating Φ curve, and the tumor-to-normal ratio of aquaporin abundance (ρ_ratio), which operates as a non-linear permissive threshold switch transitioning the malignant cell across the Φ = 1 boundary rather than a simple linear sensitivity lever.

2.6. Global Sensitivity Analysis

To determine which parameters actually control the predicted selectivity, we performed a variance-based global sensitivity analysis rather than a local one-at-a-time perturbation, since the model is nonlinear and its parameters may interact. We computed first-order and total-effect Sobol indices [65] (Equations 11 and 12), which apportion the variance of the output among the inputs as:

$$S_i={\displaystyle \frac{{\mathrm{Var}}_{X_i}\left[{\mathrm{E}}_{{\mathbf{X}}_{\sim i}}\left(Y\mid X_i\right)\right]}{\mathrm{Var}\left(Y\right)}}$$

(11)

$$S_{T_i}={\displaystyle \frac{{\mathrm{E}}_{{\mathbf{X}}_{\sim i}}\left[{\mathrm{Var}}_{X_i}\left(Y\mid {\mathbf{X}}_{\sim i}\right)\right]}{\mathrm{Var}\left(Y\right)}}$$

(12)

where Y is the selectivity factor, X_i is the i-th parameter, and X_∼i denotes all parameters except X_i. The first-order index S_i measures the contribution of X_i acting alone, while the total-effect index $S_{T_i}$ additionally captures its contributions through interactions. Input samples were drawn using the Saltelli extension of Sobol sequences [66], and the indices were estimated with the SALib library [67]. To avoid the operating-point and ratio quantities of Section 2.5 being aliased onto one another, the sensitivity design treated the absolute channel permeability of each cell and the tumour-to-normal abundance ratio as separate inputs, so that their distinct effects on the selectivity could be resolved independently.

2.7. Robustness Analysis

To assess whether the selective regime is a fragile or a generic feature of the parameter space, we complemented the Sobol analysis with a Monte Carlo robustness study. Each parameter was assigned a distribution spanning its plausible range from Table 1, and the selectivity factor and the survival outcomes were evaluated over a large ensemble of independent parameter draws. From the resulting distribution we report the median selectivity and its interquartile range, the probability that the selectivity exceeds chosen thresholds, and the probability that the two cell types fall on opposite sides of the Φ = 1 transition, the last of these quantifying how often the switch-like selective window is realized. This procedure tests robustness against simultaneous, correlated uncertainty in all parameters rather than against perturbation of any single one.

2.8. Software and Data Availability

Molecular dynamics simulations were performed with GROMACS [50] on membranes built with CHARMM-GUI [43,44,45] Trajectory analysis used MDAnalysis [56,57] The transfer-function model, sensitivity analysis, and robustness study were implemented in Python using NumPy [68], SciPy [69], Matplotlib [70], and SALib [67].

3. Results

3.1. The Malignant and Normal Membranes are Near-Indistinguishable in the Properties that Govern Passive Permeability

We first asked whether the passive lipid bilayer could, in its pristine state, furnish an intrinsic permeability advantage to the malignant cell. The two simulated membranes differed only in the modest compositional respects described in Section 2.1: both contained 512 lipids, with the astrocyte-like bilayer comprising 205 cholesterol, 95 POPC, 51 POPE, 51 PSM, 38 DOPC, 31 POPS, 26 DOPE, and 15 POPI molecules, and the U87-MG-like bilayer comprising 215 cholesterol, 98 POPC, 61 POPE, 46 PSM, 31 DOPC, 31 POPS, 15 DOPE, and 15 POPI molecules. The malignant membrane thus carried slightly more cholesterol (approximately 42% versus 40%), a small redistribution within the phosphatidylcholine and phosphatidylethanolamine pools from di-unsaturated toward mono-unsaturated species, and marginally less sphingomyelin, reflecting documented features of the glioblastoma lipidome [38,39,40]. Each composition was sampled in three independent 200 ns replicates, and all observables were averaged over the 80–200 ns interval with the three replicates treated as the units of statistical replication.

Acyl-chain ordering, the most direct structural proxy for hydrophobic core packing and the free-energy barrier encountered by a neutral permeant, was very similar between the two systems. The mean glycerophospholipid carbon–deuterium order parameter, −S_CD, was 0.251 ± 0.001 in the astrocyte membrane and 0.257 ± 0.003 in the malignant membrane, a small (about 2%) difference that does not reach conventional significance (Welch’s t-test, p = 0.07) yet separates the two cell types without overlap across replicates (Mann–Whitney U, p = 0.10). The malignant membrane was consistently the more ordered—the opposite direction required for a passive permeability mechanism, since elevated chain order tightens packing and suppresses solute permeation. The per-segment order-parameter profiles, which overlay closely across both acyl chains, are shown in Figure S2. Bilayer thickness, measured as the mean phosphate-to-phosphate distance (d_P–P) across the leaflets, was likewise indistinguishable: 4.56 nm for the astrocyte and 4.58 nm for the malignant membrane (p = 0.14; Figure S3).

The raw area per lipid initially appeared to separate the two systems—73.8 Ų per phospholipid for the astrocyte versus 75.3 Ų per phospholipid for the malignant membrane (p = 0.005)—but this parameter is confounded by cholesterol content. Because cholesterol occupies lateral cross-sectional area without contributing a headgroup phosphate, normalizing strictly to phospholipid count artificially inflates the apparent area of the more cholesterol-rich malignant membrane. When the cross-sectional area was normalized to total lipid count (phospholipids plus sterol), the ordering reversed: the cholesterol-corrected area per total lipid, ⟨A_total⟩, was 44.28 ± 0.06 Ų for the astrocyte and 43.71 ± 0.15 Ų for the malignant membrane (Figure S3). On this confound-free measure, the malignant bilayer is marginally more condensed. Finally, lateral mobility, extracted from the center-of-mass-removed, periodic-boundary-unwrapped lateral mean-squared displacement, was 1.80 ± 0.42 μm² s⁻¹ for the astrocyte and 1.72 ± 0.18 μm² s⁻¹ for the malignant membrane, showing no statistically resolvable difference (Welch p = 0.79; Figure S4). On the accessible 200 ns timescale the displacement remains sub-diffusive (α≈0.59), so these values represent effective mobilities whose malignant-versus-normal comparison, rather than absolute Fickian diffusivity, constitutes the relevant metric.

Taken together, these four observables converge on a single baseline conclusion. In their unoxidized states, the two membranes are near-identical across every structural property governing passive permeation. In the two metrics where a difference is resolvable—cholesterol-corrected packing (⟨A_total⟩) and trend-level chain order (−S_CD)—the malignant membrane is marginally more ordered and condensed. This is precisely the wrong sign for a passive permeability account of selectivity, which requires the malignant membrane to be more permeable.

Crucially, this finding must be interpreted with methodological rigor. Because the elevated cholesterol fraction was assigned a priori to reflect lipidomic literature [38,39] rather than measured directly from this specific cell line, the marginal increase in order and condensation is in substantial part the direct, predictable physical consequence of that modeled sterol enrichment. We therefore designate this result the unoxidized bilayer baseline null (NULL#1): under intact, non-peroxidized conditions, glioblastoma-like compositional shifts provide no passive permeability advantage to the tumor. This null demonstrates that pristine bilayer biophysics cannot explain selective killing. Consequently, the mechanistic burden shifts away from pristine membrane order toward two alternative possibilities: facilitated transport through regulated protein conduits (AQP8), or potential differences in cellular susceptibility to RONS-induced lipid peroxidation and pore formation under sustained plasma challenge.

While formal two one-sided equivalence testing (TOST; Figure S5) robustly confirms structural equivalence within ±5% for bilayer thickness and total area per phospholipid, lateral diffusion (D_lat) is not tightly bounded within the ±10% equivalence margin at n = 3 (90% CI:[−29.4%, +35.0%]); hence, dynamic fluidity shows no statistically resolvable difference (p = 0.79) but represents the least statistically constrained membrane observable.

3.2. Transfer-Function Model

Having found no passive-membrane basis for selectivity, we turned to the transfer-function model, in which the cytoplasmic dose is governed by the single dimensionless transport number Φ of Equation (4). Evaluated at the baseline parameters of Table 1, the two cell types occupy qualitatively different regimes. The astrocyte, with its high clearance rate and modest channel permeability, has Φ_astro = 0.08, placing it deep in the clearance-dominated regime; its steady-state cytoplasmic concentration is correspondingly only C_in,astro = 7.4 μM, or 7.4% of the imposed extracellular concentration of 100 μM. The malignant cell, with its reduced clearance and elevated channel permeability, has Φ_u87 = 2.0, placing it in the influx-dominated regime; its steady-state concentration is C_in,u87 = 66.7 μM, two-thirds of the extracellular value. The two cells therefore sit on opposite sides of the Φ = 1 transition that separates the clearance-dominated from the influx-dominated regime (Figure 1).

This regime separation is the structural heart of the result. Selectivity in this model is not a gentle graded difference between two cells responding similarly to exposure; it is a consequence of the two cells lying on opposite sides of a saturating nonlinearity. The astrocyte, operating on the steep initial portion of the Φ/(1 + Φ) curve, converts even substantial extracellular exposure into only a small internal concentration, while the malignant cell, operating past the knee of the same curve, admits most of the extracellular concentration into its cytoplasm. The biological inputs that place each cell where it sits—the clearance asymmetry attributed to differential antioxidant and membrane-catalase activity [35,36] and the elevated peroxiporin permeability of the malignant cell [27,28]—thereby act jointly to straddle the transition rather than merely to shift both cells along a common response.

3.3. The Dose-Based Selectivity Factor

To make the origin of the selectivity transparent, we factorized the dose-based selectivity factor, defined as the ratio of steady-state cytoplasmic concentrations:

$$S_{\mathrm{dose}}={\displaystyle \frac{C_{\mathrm{in},\mathrm{u}87}}{C_{\mathrm{in},\mathrm{astro}}}}$$

At baseline and in the absence of channel gating, S_dose = 9.0. This value decomposes cleanly into two multiplicative contributions at the level of the transport number. The total membrane permeability of the malignant cell exceeds that of the astrocyte by a factor of 2.5—the influx amplification—while its clearance rate is lower by a factor of 10—the clearance asymmetry. The product of these is a 25-fold ratio in the transport number Φ. Crucially, however, the 25-fold ratio in Φ does not translate into a 25-fold ratio in concentration: because the malignant cell operates beyond the knee of the saturating transfer function, further increases in its Φ yield diminishing increases in its cytoplasmic concentration, whereas the astrocyte responds nearly linearly. The saturating nonlinearity therefore compresses the 25-fold Φ ratio into a 9.0-fold concentration ratio. This compression is not a weakness of the mechanism but a structural feature of it: the same nonlinearity that compresses the ratio is what produces the switch-like separation of regimes in the first place.

It is essential to be explicit about which numerical selectivity is being quoted, because the model yields three distinct values that are easily conflated. The first is the ungated steady-state dose ratio at baseline, S_dose = 9.0, just derived. The second, S_dose = 6.9, is the value obtained once the redox gating of AQP8 is included (Section 3.4); gating lowers the malignant-cell dose and therefore reduces the selectivity. The third, a median of 6.2 with an interquartile range of 4.6 to 8.6, is the central tendency of the selectivity over the full Monte Carlo ensemble of plausible parameter sets (Section 3.6). Throughout, we treat 6.9 as the representative point estimate, since it incorporates the gating that the channel biology requires, and 6.2 as the ensemble-robust estimate. The selectivity map of Figure 2, which shows S_dose across the plane of channel permeability and clearance rate, makes clear that the baseline operating point lies within a broad region of high selectivity rather than on a narrow ridge.

3.4. Redox Gating of the Channel Attenuates Rather Than Amplifies Selectivity

We next incorporated the cysteine-dependent gating of AQP8 described by Equation (5), with the half-closing concentration and Hill coefficient of Table 1. Solving the model self-consistently—the gating factor depends on the cytoplasmic concentration, which in turn depends on the gating factor—the malignant-cell steady-state concentration falls from 66.7 μM in the ungated case to 50.4 μM once gating is active, corresponding to a residual gating factor of approximately 0.39, that is, the channel is roughly 60% closed at the malignant cell’s operating concentration. The astrocyte, whose cytoplasmic concentration of 7.4 μM lies well below the half-closing concentration of 40 μM, is essentially ungated, with a gating factor near 0.97. The dose-based selectivity consequently falls from 9.0 to 6.9. The full time courses of cytoplasmic concentration are shown in Figure 3: the malignant cell with and without gating, and the astrocyte, which is essentially ungated and therefore drawn once.

The direction of this effect is mechanistically informative. Because the gating responds to the internal oxidant level, it acts most strongly on the cell that admits the most oxidant—the malignant cell—and barely at all on the cell that admits little. The redox gating of the peroxiporin [31] therefore functions as a direction-dependent negative-feedback brake that preferentially throttles influx into the very cell the mechanism is selecting against, attenuating rather than reinforcing the selectivity. This stands against any intuition that channel gating might sharpen the discrimination between cell types; in this model it does the opposite, and a substantial fraction of the raw selectivity is lost to it. That selectivity nonetheless survives gating at a factor of 6.9 indicates that the underlying clearance and operating-point asymmetries are large enough to absorb the attenuation.

3.5. Determinants of Selectivity: Clearance and Operating Point Dominate

The factorization of Section 3.3 identifies the individual contributors to selectivity but does not establish their relative importance across the multidimensional parameter space. To determine which parameters quantitatively control the cytotoxic divergence, we computed variance-based Sobol total-effect indices (S_T) and first-order indices (S₁) across the complete parameter ensemble using the decoupled sampling design of Section 2.6, which explicitly separates the baseline channel permeability of the normal cell from the tumor-to-normal channel abundance ratio (ρ_ratio).

The resulting sensitivity ranking is illustrated in Figure 4. The largest share of output variance is governed by the absolute channel permeability of the astrocyte (P_AQP8,astro, S_T = 0.34, S₁ = 0.19), followed by the astrocyte clearance rate (k_astro, S_T = 0.28, S₁ = 0.20), the clearance ratio between the two cell types (k_ratio = k_astro/k_U87, S_T = 0.23, S₁ = 0.17), and the background lipid permeability (P_lip, S_T = 0.18, S₁ = 0.08). In sharp contrast, the tumor-to-normal AQP8 abundance ratio (ρ_ratio) carries a total-effect index of only S_T = 0.078 (S₁ = 0.041), comparable to that of the extracellular exposure concentration (C_out, S_T = 0.076, S₁ = 0.031), while the gating parameters exhibit minor contributions (n_Hill, S_T = 0.048; C_half, S_T = 0.047).

This order requires careful physical and biological interpretation. The dominant parameter, P_AQP8,astro, represents the absolute conductance that fixes the operating point of the non-malignant cell near the steep, linear foot of the saturating transfer function (Φ_astro = 0.08). Because the astrocyte operates in this sub-saturation regime, small absolute shifts in its influx capacity produce large relative perturbations in the resulting concentration ratio. Conversely, the modest sensitivity index of ρ_ratio does not indicate that channel over-expression is biologically dispensable or secondary. Rather, because the malignant cell operates beyond the knee of the transfer function (Φ_U87 = 2.0), incremental additions to channel density yield diminishing marginal increases in C_in, as internal concentration asymptotically approaches C_out. AQP8 over-expression functions as an indispensable non-linear permissive threshold switch: it is the primary physical prerequisite that transitions the malignant cell from the clearance-dominated into the influx-dominated regime (Φ > 1). Once across this saturation threshold, the mathematical sensitivity to further channel additions naturally plateaus, allowing the clearance asymmetry to govern the ultimate divergence in intracellular dose.

To evaluate whether the modeled 10-fold clearance asymmetry (k_astro/k_U87 = 10) reflects physiologically achievable conditions, we grounded the survival model directly in the empirical single-cell clearance framework of Doskey et al. (2016) [63] (Figure S6). In native biological systems, resting H₂O₂ clearance rate constants (k_cell) exhibit modest baseline disparities between normal tissues and tumors, with an average ratio of approximately 1.8× and intra-tumoral extremes reaching ∼5× [63]. Simulating cellular survival across this empirical clearance plane demonstrates that a basal 1.8× clearance deficit is entirely insufficient on its own to achieve tumor eradication: at 100 μM for 3 h, a 1.8× asymmetry leaves >80% of U87-MG cells viable (Figure S6). Marked cytotoxic selectivity (S_dose > 5, tumor kill >90%) emerges only when the clearance ratio crosses a critical non-linear threshold of approximately 4× to 5×, becoming profound at the 10× operational baseline (1.9% tumor survival; Figure S6).

This threshold behavior demonstrates that resting, endogenous antioxidant deficits in glioblastoma cannot account for the selectivity of plasma-activated water. Instead, high selectivity strictly requires an active, exposure-dependent amplification—specifically, the tumor-specific inactivation of protective membrane-associated catalase formalized by the Bauer mechanism, which drives the operational clearance ratio from its modest basal state (1.8×) across the 5× threshold into the high-kill regime.

3.6. The Selective Regime is Robust to Joint Uncertainty in All Parameters

A selectivity that held only at a single hand-chosen parameter set would be of little interest. We therefore tested whether the selective regime persists under simultaneous uncertainty in every parameter, drawing a large Monte Carlo ensemble in which each parameter varied independently across the plausible range of Table 1, and evaluating the selectivity and the regime classification for each draw. The distribution of the selectivity factor over this ensemble (Figure 5) is centred on a median of 6.2 with an interquartile range from 4.6 to 8.6. The probability that the selectivity exceeds two-fold is 0.99, and the probability that it exceeds five-fold is 0.69. The probability that the two cell types fall on opposite sides of the Φ = 1 transition—that is, that the switch-like selective window is realized, with the astrocyte clearance-dominated and the malignant cell influx-dominated—is 0.92 (Figure 6).

These figures indicate that selectivity in this model is a generic rather than a fine-tuned property. Across the overwhelming majority of plausible parameter combinations, the malignant cell accumulates at least twice the cytoplasmic dose of the normal cell, and in roughly two cases out of three it accumulates at least five times as much; in more than nine cases out of ten the two cells occupy distinct transport regimes. The selectivity is therefore not contingent on a precise choice of parameters but emerges robustly from the qualitative structure of the model—a saturating transfer function straddled by two cells of differing clearance and permeability. The interquartile range, spanning a factor of roughly two around the median, quantifies the residual sensitivity: the existence of selectivity is robust, while its precise magnitude remains parameter-dependent, which is the appropriate level of confidence for a reduced model of this kind.

3.7. Cumulative Dose Integration Amplifies Cytotoxic Selectivity

The steady-state concentration difference translates into differential cytotoxicity through the cumulative accumulation of intracellular oxidant dose over exposure time. Integrating the cytoplasmic concentration over a standard 3-h treatment duration [10] and applying the cooperative Hill survival relation calibrated against the empirical glioblastoma kill curve of Kurake et al. (2016) [64] (Section 2.5) yields the cumulative-dose and viability trajectories shown in Figure 7. At an extracellular exposure of C_out = 100 μM, the astrocyte accumulates an integrated dose of D≈7.9 × 10⁴ μM ⋅ s, remaining well below the calibrated transition threshold (D₅₀ = 3.06 × 10⁵ μM ⋅ s) and retaining virtually complete viability (≈100%). In contrast, the malignant cell—even with negative-feedback channel gating active—accumulates D≈5.4 × 10⁵ μM ⋅ s, traverses the steep foot of the cooperative survival curve, and is reduced to approximately 1.9% survival, generating a ∼52-fold survival selectivity. Consequently, an exposure protocol that spares the non-malignant population eliminates roughly 98% of the malignant burden.

Because the calibrated cellular response is sigmoidal and cooperative rather than single-hit exponential—a formulation decisively favored by information criteria (ΔAICc≈12 relative to single-hit; Figure S7) and statistically comparable to a Weibull model (ΔAICc ≤ 2)—the underlying ∼7-fold steady-state concentration ratio is non-linearly amplified through the threshold regime into a near-complete viability divergence. To verify that the predicted sparing of non-malignant cells does not represent an unconstrained numerical extrapolation, we contextualize this response against the non-malignant experimental anchor from Kurake et al. (2016) (MCF10A, ≥93% survival at 64 μM; Figure S8).

While MCF10A serves as a general phenotypic benchmark for non-malignant mammalian cell tolerance rather than an astrocyte-specific lineage match, it empirically confirms that non-transformed cells withstand extracellular peroxide fluxes that are lethal to malignant glioma under identical assay configurations. In our transport model, astrocyte sparing is mechanistically rooted in its intrinsic, uncollapsed clearance capacity (k_clear = 50 s⁻¹), which restricts the accumulated cytoplasmic dose to the sub-lethal shoulder of the response curve.

Finally, we assembled the individual results into the unified conceptual framework summarized in Figure 9. Extracellular H₂O₂ delivered by plasma-activated water enters the cell through both background lipid permeation (P_lip) and facilitated peroxiporin channels (P_AQP8), with the latter subjected to negative-feedback attenuation via Cys53 oxidation (Figure 9A). The intracellular concentration is governed by mass-balance kinetics parameterized by the effective membrane transport rate constant k_m = (P_AQP8 ⋅ g + P_lip)A/V and first-order enzymatic clearance k_clear, defining the dimensionless transport number Φ = k_m/k_clear that yields a saturating steady-state dose relation (Figure 9B). This framework segregates cells into two distinct operational regimes: normal astrocytes (catalase intact) remain protected deep within the clearance-dominated regime (Φ = 0.08≪1), whereas glioblastoma cells undergo a plasma-induced collapse of protective membrane-associated catalase that transitions them across the threshold into the influx-dominated, lethal regime (Φ = 2.0 > 1) (Figure 9C). Global sensitivity analysis confirms that clearance asymmetry (k_astro/k_U87 = 10) is the dominant lever dictating this selectivity, whereas AQP8 over-expression acts as a permissive non-linear switch (S_T≈0.08) and background lipid permeability establishes a hard ceiling (Figure 9C). Integrating cytoplasmic concentration over time ($D={\int }_0^T{C}_{\mathrm{in}}\left(t\right)\hspace{0.17em}\mathit{dt}$) and mapping it through a Hill-calibrated survival threshold (D₅₀ = 3.06 × 10⁵ μM ⋅ s,h = 6.82) establishes an exposure contrast that defines a calculable, platform-testable therapeutic window (Figure 9D). This causal chain reframes PAW selectivity not as a static property of unoxidized membranes or isolated aquaporin abundance, but as an emergent property of a regulated transport-clearance switch.

3.8. Therapeutic Window Mapping Across Dose–Time Plane

Generalizing from single representative exposure, we mapped survival across the plane of exposure time and extracellular concentration for both cell types and identified the therapeutic window as the region in which malignant survival falls below a lethality threshold while normal survival remains high (Figure 8). The window is bound below by exposures too brief or too diluted to accumulate a lethal malignant dose and bounded above by exposures severe enough to begin compromising the normal cell. Between these bounds lies a substantial region in which the model predicts selective elimination of the malignant cell with preservation of the normal cell, consistent with the existence of a treatment regime exploiting the regime separation. The breadth of this window in the model reflects the robustness documented in Section 3.6; because the survival relation is calibrated to a measured kill curve (Section 2.5), its boundaries (here approximately 60 to 350 µM over a three-hour exposure) are quantitative predictions, though they remain specific to the calibration data and the model’s operating-point assumptions rather than prescriptive clinical dosing.

3.9. Unified Causal Chain of Plasma Selectivity

Finally, we assembled the individual results into the unified causal chain summarized in Figure 9, which traces the path from the extracellular plasma chemistry, through the membrane transport step composed of channel and lipid contributions, to the resulting cytoplasmic dose, and onward through cumulative dose to the survival phenotype, with the redox gating entering as feedback on the transport step. Read along this chain, the results localize the origin of selectivity unambiguously. The first link—the passive lipid membrane—contributes no selectivity and a slight counter-selective bias (Section 3.1). The decisive links are the regime separation produced jointly by the clearance asymmetry and the operating-point permeabilities (Section 3.2 and Section 3.3 and Section 3.5), partially attenuated but not abolished by channel gating (Section 3.4), and robust to joint parameter uncertainty (Section 3.6). The over-expression of aquaporins between cell types, often invoked as the proximate cause of plasma selectivity, enters the chain as a necessary but quantitatively secondary contributor (Section 3.5). The chain thereby reframes plasma selectivity for glioblastoma not as a property of any single membrane component but as an emergent consequence of regulated transport operating across a saturating nonlinearity, with the clearance asymmetry as its dominant driver.

4. Discussion

4.1. Three Convergent Nulls Localize the Selectivity Away from the Passive Membrane

We set out to determine whether the documented selectivity of plasma-activated water for glioblastoma over normal astrocytes [11,12] originates in membrane biophysics, and we tested three candidate formulations of that hypothesis. Each returned a null result, and the three nulls are convergent rather than redundant: they rule out distinct passive routes and collectively displace the origin of selectivity onto regulated transport and enzymatic clearance.

The first null concerns the baseline structural state of the unoxidized membranes. Across six independent 200 ns trajectories, the glioblastoma-like and astrocyte-like bilayers were statistically indistinguishable in bulk bilayer thickness (4.56 versus 4.58 nm) and effective lateral mobility (1.80 versus 1.72 µm² s⁻¹), while in acyl-chain order (−S_CD of 0.251 versus 0.257) and cholesterol-corrected area (44.28 versus 43.71 Ų per total lipid), the malignant membrane was marginally more ordered and condensed. A bilayer that is slightly more condensed is the opposite of what an intrinsic permeability hypothesis requires.

This result must be assigned its precise epistemic weight. Because the malignant bilayer was constructed to reflect the documented cholesterol enrichment of glioblastoma [38,39] rather than profiled from an isogenic cell pair, the slight increase in structural condensation is the direct thermodynamic outcome of incorporating additional sterol molecules. Rather than claiming an unexpected biological discovery, this simulation confirms an essential physical constraint: the documented lipidomic direction of glioblastoma (elevated cholesterol) works in the wrong direction to promote passive transport, as cholesterol systematically raises the free-energy barrier to H₂O₂ permeation across phospholipid bilayers [23,24,58].

Equally important is what this unoxidized baseline null (NULL#1) does and does not exclude. It firmly eliminates the hypothesis that glioblastoma cells are selectively killed because their pristine plasma membranes possess higher basal permeability. However, it does not rule out dynamic membrane mechanisms that occur only during plasma exposure, such as differential accumulation of oxidized phospholipid species (e.g., hydroxy- or hydroperoxy-lipids) and subsequent transient pore formation [18,19,20]. If malignant membranes possess lower lipophilic antioxidant defenses or altered unsaturated fatty acid pools that make them more vulnerable to peroxidative permeabilization, that mechanism operates downstream of initial stress. What our MD data demonstrate is that pristine membrane biophysics provides zero discriminatory power, establishing a rigorous negative control for resting bilayer states.

The second null concerns the intrinsic permeation barrier of the lipid core. Even setting structural equivalence aside, the lipid-phase permeability of small polar oxidants is modest and only weakly modulated by acyl chain variations, meaning the lipid bilayer does not chemically discriminate among physiological RONS fluxes. Direct functional measurements confirm that H₂O₂ permeability across cholesterol-containing phospholipid vesicles is moderate and independent of aquaporin expression in erythrocytes [58], a conclusion corroborated by free-energy simulations of neutral reactive species diffusing across hydrocarbon cores [19,26]. The lipid route thus operates in the model as a non-selective background conductance (P_lip) rather than an active discriminatory filter.

The third null concerns bulk antioxidant clearance treated as a simple linear scalar. A basic argument asserting that “tumor cells have lower antioxidant capacity and therefore accumulate more oxidant” fails to generate the switch-like survival selectivity observed experimentally: in the absence of a saturating transport nonlinearity, a linear difference in clearance scales internal concentrations proportionally across both cell types without bifurcating them into protected versus lethal states. Clearance asymmetry is indeed the decisive driver of selectivity, but only when coupled to the saturating transport number Φ, which translates a quantitative asymmetry into a qualitative regime switch.

Read in combination, the three nulls rule out pristine bilayer biophysics, non-specific lipid permeation, and simple scalar clearance. They localize the selectivity of plasma-activated water to the coordinated interaction of three specific elements formalized in the redox oncology framework [36,37]: channel-mediated H₂O₂ influx via peroxiporins (AQP8), the tumor-specific, plasma-induced inactivation of protective membrane-associated catalase, and the resulting transport regime separation across the cellular plasma membrane.

4.2. A Transfer Function for AQP8-Gated Influx, and the Regime Switch

The surviving mechanism is captured compactly by the steady-state transfer function of Equations (1)—(4), in which the killing-relevant cytoplasmic concentration is set by the balance of channel- and lipid-mediated influx against first-order clearance, and is governed by the single dimensionless transport number Φ. The value of this formulation is that it converts a qualitative, multi-factor narrative into a quantity with a sharp structural feature: a saturating nonlinearity with a transition at Φ = 1 separating a clearance-dominated regime, in which the cell admits only a small fraction of the extracellular oxidant and is protected, from an influx-dominated regime, in which the cell floods toward the extracellular concentration. At the baseline operating point, the astrocyte sits deep in the protected regime (Φ = 0.08, cytoplasmic concentration 7.4 µM) and the malignant cell sits in the flooded regime (Φ = 2.0, cytoplasmic concentration 66.7 µM). Selectivity is therefore not a graded difference between two similarly behaving cells but the consequence of the two cells lying on opposite sides of a threshold—the astrocyte with intact catalase remaining below it, the tumor cell, combining elevated AQP8 with plasma-collapsed catalase, driven across it. This regime-switch picture is the natural quantitative expression of the Bauer mechanism [35,36]in which the transient inactivation of protective catalase is precisely what collapses the malignant cell’s clearance and tips it from the protected to the flooded regime, and in which the subsequent H₂O₂ influx that depletes glutathione is explicitly aquaporin-mediated [37]. The dynamic form of the model (Figure 3) makes the same point as a live-cell readout: the malignant trajectory floods while the astrocyte trajectory does not.

4.3. Global Sensitivity Reassigns the Dominant Role from Channel Abundance to Clearance Asymmetry

The decomposition of steady-state dose into a 2.5-fold influx amplification and a 10-fold operational clearance asymmetry identifies the mathematical components of selectivity, but the global variance-based sensitivity analysis fundamentally revises their mechanistic hierarchy. Across the plausible parameter space, the primary determinants of selectivity are the baseline operating point on the Φ switch curve and the clearance asymmetry, while the tumor-to-normal AQP8 abundance ratio (ρ_ratio) accounts for less than 8% of the total output variance (S_T≈0.08).

This finding challenges the prevailing assumption that peroxiporin over-expression is the direct, linear driver of cancer selectivity in plasma oncology. The low sensitivity index of ρ_ratio is a direct mathematical consequence of the saturating transfer function: once the malignant cell crosses the Φ = 1 threshold, the membrane becomes transport-saturated, and further increases in channel density push the internal concentration toward C_out with rapidly diminishing returns. AQP8 over-expression is therefore a mandatory enabling condition—the physical conduit through which exogenous peroxide accesses the cytoplasm [27,28]—rather than the primary tuning knob that governs the magnitude of kill selectivity.

Crucially, our grounding of clearance rates in the empirical measurements of Doskey et al. (2016) resolves an important mechanistic question regarding the origin of this asymmetry. Unstressed tumor cells do not possess a static 10-fold clearance deficit; their intrinsic removal capacity is only ∼1.8-fold lower than normal counterparts on average. Our model reveals that if plasma-activated water relied solely on this native, basal difference, therapeutic selectivity would fail, sparing over four-fifths of the tumor population under standard treatment doses (Figure S6).

High selectivity is therefore an induced, dynamic state rather than a pre-existing physiological feature. It requires the extracellular chemical cocktail generated by plasma—specifically the synergistic interaction between H₂O₂ and ${\mathrm{NO}}_2^{-}$ generating transient peroxynitrite species—to selectively inactivate the membrane-associated catalase shield of tumor cells, as demonstrated experimentally by Bauer and colleagues [35,36,37]. This dynamic collapse widens the effective clearance gap from the basal 1.8× regime beyond the non-linear tipping point (≥5×), driving the malignant cell into the flooded, lethal state while leaving the non-malignant astrocyte below the critical threshold.

We emphasize two biological caveats. First, established glioma cell lines cultured in vitro may alter aquaporin expression profiles relative to primary tumor tissue in situ [60]. Second, AQP8 also localizes to the inner mitochondrial membrane, where it participates in organellar redox signaling [62]. These considerations do not compromise the model’s core architecture, but they underline the necessity of directly measuring functional cell-surface channel densities and real-time clearance decay kinetics under active plasma exposure (Section 4.8).

4.4. Lipid Permeability as a Ceiling, and an Unexpected Coupling Between the Membrane and Channel Work

The sensitivity analysis also revises the status of the lipid route, in a way that connects the two computational strands of this work. Although the membrane simulations establish that lipid permeability is non-selective, the background lipid permeability is not merely a null term: it sets a ceiling on the achievable channel-based selectivity (Figure 4). Because the lipid route admits H₂O₂ into both cell types indiscriminately, a sufficiently large background permeability drives both cells toward the influx-dominated regime and erases the contrast carried by the channel and the clearance asymmetry. The measured lipid permeability is therefore a quantitative input to the selectivity model rather than an isolated negative result. This is an unexpected and useful coupling between the membrane analysis and the transport model that should be stated explicitly: the selectivity window exists only because the non-selective background is small relative to the gated channel route—precisely the regime in which aquaporins, rather than the lipid phase, dominate the delivery of H₂O₂ to the cytoplasm [27,30]. The two strands therefore play complementary roles: the membrane simulations establish that this background is non-selective, while its magnitude—the quantity that bounds the channel-based selectivity—is set by direct experimental measurement (Section 2.2) rather than by the simulations themselves.

4.5. Redox Gating is a Direction-Dependent Brake, not a Driver

A central biophysical subtlety resolved by our model concerns the functional role of AQP8 redox gating. Under oxidative challenge, AQP8 permeability is reversibly restricted through the oxidation of conserved cysteine residues (Cys53) [31]. When incorporated into our transfer-function framework as a steady-state negative-feedback regulator, this gating depresses the malignant steady-state cytoplasmic concentration from 66.7 μM to 50.4 μM and attenuates the dose-based selectivity factor from 9.0 to 6.9, without abolishing the underlying regime separation (Figure 3).

The biological interpretation of this effect depends decisively on the vectorial direction of the transmembrane H₂O₂ gradient. Because aquaporins are passive, bidirectional channels, the protective versus cytotoxic valence of channel gating is inverted between endogenous and exogenous stress regimes:

Internal oxidative stress (Endogenous gradient): When H₂O₂ is generated intracellularly (e.g., via mitochondrial respiration or metabolic dysregulation), the concentration gradient drives peroxide efflux out of the cell. In this setting, an open channel facilitates toxic waste clearance; channel closure traps oxidants within the cytosol, whereas a non-oxidizable cysteine mutant maintains efflux and enhances stress survival [31].

Plasma-activated water exposure (Exogenous flood): Under therapeutic PAW challenge, the chemical gradient is radically reversed (C_out≫C_in). The open channel now operates as a conduit for lethal influx. Under these conditions, Cys53 oxidation acts as an emergency cellular brake that throttles further entry, partially shielding the malignant cell from external flooding.

This gradient-dependent inversion reconciles why non-oxidizable AQP8 mutants confer protection against endogenous metabolic stress [31] while simultaneously predicting that such mutants would exacerbate toxicity under plasma treatment by disabling the tumor’s protective brake.

Furthermore, the kinetic timescale separation established in Section 2.3 provides critical insight into the real-time dynamics of this mechanism. Because covalent thiol oxidation is kinetically bounded (τ_gate∼10¹ − 10² s), channel closure is non-instantaneous. During the initial phase of exposure, the glioblastoma cell is subjected to unrestricted channel transport (g≈1), driving rapid cytoplasmic accumulation at the maximal rate of 66.7 μM before the Cys53 brake throttles influx down to its steady-state operating level of 50.4 μM (Figure 3). The quasi-steady-state selectivity factor (S_dose = 6.9) and the resulting therapeutic window (Figure 8) are therefore highly robust: they represent a conservative baseline that survives channel attenuation, whereas real-time non-equilibrium kinetics deliver an even larger early lethal burden to the tumor.

We treat the quantitative parameterization of the gating transition as a priority for future atomistic refinement (Section 4.7). However, we note that previous molecular dynamics simulations of peroxiporin bilayers failed to resolve spontaneous oxidative gating on accessible nanosecond-to-microsecond trajectories [26], consistent with our finding that channel closure is governed by slow chemical kinetics rather than instantaneous conformational fluctuation.

4.6. From Modest Biophysical Selectivity to a Calculable Therapeutic Window

The biophysical selectivity the model produces is modest—a roughly seven-fold difference in cytoplasmic dose—but its translation into a survival difference is not, because the cytotoxic response to cumulative peroxide is cooperative rather than single-hit. Integrating the cytoplasmic concentration over the exposure and applying the Hill survival relation calibrated to a measured glioblastoma kill curve (Section 2.5) amplifies the seven-fold dose ratio through the threshold region into a near-complete viability gap: for a three-hour exposure at 100 µM, the model predicts essentially full astrocyte survival (about 100%) against roughly 1.9% malignant survival, a roughly 52-fold selectivity (Figure 7). Sweeping the extracellular concentration against the treatment time then defines a therapeutic window (Figure 8) bounded below, by a minimum dose-and-time below which the malignant kill is incomplete (here near 60 µM at three hours), and above, by a maximum beyond which the astrocyte too is killed (here near 350 µM)—the dose–time trade-off characteristic of a genuine therapeutic index. This connects the model directly to the experimental literature, in which plasma-activated media kill glioblastoma cells preferentially over astrocytes at intermediate doses [11,12] the window is the quantitative expression of that observation. Because the survival relation is now fixed by calibration rather than by an assumed lethal-dose scale, the boundaries of Figure 8 are quantitative predictions of the safe operating region—directly testable against kill curves on the existing plasma platform—rather than an illustration of its structure; a bootstrap over the calibration kill curve leaves the malignant-to-normal selectivity at least 17-fold across resampled fits.

4.7. A Structural Ceiling on Atomistic Refinement, which is Itself a Result

We originally sought to replace the literature-bracketed single-channel permeability and gating factor with first-principles values by computing the H₂O₂ potential of mean force through AQP8 directly. However, executing this free-energy calculation was precluded by a decisive structural limitation that is itself informative. No experimental atomic structure of AQP8 exists; the available model places the aromatic/arginine selectivity filter (Arg213) in a closed, pore-occluding conformation, and a side-chain repack of that residue followed by restrained equilibration relaxed back to the closed state. Decisively, a conformational ensemble generated by deep-learning structure prediction with input subsampling and inference dropout—240 models spanning three input depths—never opened the filter, with a maximum selectivity-filter radius of 0.63 Å against the roughly 1.0–1.1 Å required for single-file conduction. The closed filter is therefore robust across the model’s entire learned conformational distribution, not a single-seed artefact. This vindicates the prior structural caution that AQP8 models are built closed by pore-blocking residues [42]—we note that the occluding residue differs between that work and ours, an aromatic block versus the arginine of the aromatic/arginine region, but the conclusion that the resting filter is conformationally restricted is shared—and it is consistent with AQP8 being a tightly regulated peroxiporin whose permeability is gated rather than constitutive [31]. The practical consequence is that an absolute atomistic permeability for AQP8 cannot at present be obtained from structure, and the experiment-anchored transfer function used here is therefore the rigorous choice rather than a fallback. This negative result is worth reporting in its own right, because it bears on the growing use of predicted structures for channels that lack experimental coordinates. Should an experimental AQP8 structure become available—as has recently occurred for other aquaporins—the present framework specifies exactly which two quantities, the single-channel permeability and the cysteine gating factor, an atomistic calculation would refine.

4.8. Limitations, Falsifiability, and an Experimental Programmed

The model makes deliberately simple assumptions, and we state them plainly. It treats the cytosol as well-mixed, clearance as first-order, the gate as quasi-equilibrium, and the cell as a single plasma-membrane compartment, neglecting the mitochondrial pool of AQP8 and the spatial organization of redox relays that more detailed treatments of peroxide signaling emphasize [17]. The analysis is framed around hydrogen peroxide as the effector species; peroxynitrite and other reactive oxygen and nitrogen species, and any peroxide-independent cytotoxic routes, lie outside its scope and are treated only qualitatively. The compositions of the two membranes are modelled rather than measured, so the membrane null is conditional on the chosen lipidome, as discussed in Section 4.1, even though the direction of the cholesterol effect makes it robust. The variance-based attribution, though not its two robust conclusions, depends on the chosen parameter ranges. None of these simplifications are hidden, and each defines a direction for refinement.

Crucially, the framework is falsifiable, and its predictions are specific. It predicts that knockdown of AQP8 should attenuate plasma-driven intracellular accumulation and killing in a dose-dependent manner; that if the tumor-to-normal AQP8 ratio is near unity on the actual cell lines, the density term should vanish, leaving clearance asymmetry as the sole driver; that a non-oxidizable gating-cysteine mutant should alter the flood kinetics under plasma exposure only to the extent the brake is operative; and that a sufficiently permeable lipid background would abolish the selectivity. These predictions map onto a five-experiment validation programmed executable on the existing platform within twelve months: quantification of AQP8 surface density on the malignant versus normal cells, distinguishing the surface from the mitochondrial pool [62]; measurement of catalase and glutathione-peroxidase activity and H₂O₂-decay kinetics with and without plasma exposure, to fix the clearance asymmetry [36]; live-cell ratiometric imaging of intracellular peroxide across a dose series, to locate the Φ switch; AQP8 knockdown, to establish causality; and the gating-cysteine mutant, to isolate the brake. The two dominant model parameters—the clearance asymmetry and the absolute channel permeability—are precisely those the programmed measures first, so the experiments are ordered by the sensitivity analysis itself.

A standard consideration in translational plasma oncology is whether the assumption of constant extracellular oxidant concentration (C_out = 100 μM) holds during multi-hour in vitro incubations, as cancer cells progressively consume extracellular H₂O₂ (Yan et al., 2017) [10]. To address this, we developed a self-consistent depletion model coupling finite extracellular bath consumption to cytoplasmic uptake (Figure S10, left). Recalibrating the cooperative survival threshold under this decaying regime demonstrates that selective kill efficacy is strictly preserved or marginally enhanced (62× to 64× selectivity vs 52× in the unperturbed system; Figure S10, right), because malignant cells absorb their lethal cumulative dose during the initial high-flux transient before bath exhaustion occurs. Consequently, the constant-C_out framework adopted in our primary derivations represents a conservative lower bound on therapeutic selectivity.

5. Conclusions

The picture that emerges is that plasma selectivity for glioblastoma is not a property of any single molecule but an emergent property of a regulated transport system. A non-selective lipid background, held small, admits peroxide into both cell types without discrimination and indeed slightly favors the normal cell on account of the tumor’s cholesterol enrichment.

A gated AQP8 conduit provides the indispensable physical pathway for cellular influx; while its differential expression behaves as a saturating variable in variance-based sensitivity analysis, channel over-expression serves as a mandatory permissive threshold switch that transitions the malignant cell into the influx-dominated regime (Φ > 1), while its redox gating acts as a negative-feedback brake that partially attenuates the flood. The decisive element governing the magnitude of the therapeutic gap is the tumor-specific, exposure-dependent collapse of catalase-based clearance, which traps the malignant cell in a lethal state while the normal cell, with its clearance intact, remains protected below the threshold.

The biophysical selectivity this produces is modest at steady state, but its cooperative amplification into survival yields a calculable therapeutic window whose boundaries the present model turns into testable, platform-ready predictions. By contextualizing AQP8 over-expression as an enabling threshold switch and establishing dynamic clearance collapse as the primary divergence lever, and by quantifying a structural ceiling that justifies an experiment-anchored rather than a first-principles treatment, the analysis offers both a specific mechanistic claim about glioblastoma and a transferable template for dissecting selectivity in plasma oncology more broadly.