Background: Electrochemical therapy (ECT) offers a minimally invasive approach to tumor ablation by using direct current to generate cytotoxic pH extremes. However, clinical adoption has been hindered by a lack of predictive tools for real-time, patient-specific dose planning, as high-fidelity multi-physics models are computationally prohibitive for intraoperative use.

Objective: This work presents a comprehensive computational framework to elucidate the fundamental physics of ECT and to overcome the barrier to clinical translation. We develop a rigorous one-dimensional axisymmetric model of a needle electrode in hepatic tissue, solved with COMSOL Multiphysics, which couples Nernst-Planck ion transport with Butler-Volmer electrode kinetics. Furthermore, we introduce a novel Proper Orthogonal Decomposition/Discrete Empirical Interpolation Method (POD-DEIM) reduced-order model (ROM) designed to deliver real-time simulation capability.

Results: The full-order model reveals critical and non-intuitive dynamics. A therapeutically lethal pH below 2 is achieved after approximately 2000 seconds, consistent with experimental data. We uncover a surprising phenomenon: the peak proton concentration occurs not at the anode surface, but at an interior point, a result of the complex interplay between time-varying potential, proton migration, and chloride depletion. This chloride depletion drives a characteristic three-phase evolution in current density, marking a mechanistic shift from chlorine-dominated to oxygen-dominated electrolysis. The POD-DEIM ROM, utilizing only 20 modes, reproduces the full-order pH and concentration fields with errors below 1% while achieving a transformative speedup factor of 10²–10³.

Conclusion: This work provides a complete visual and mathematical atlas of ECT physics, from single-electrode dynamics to multi-needle arrays. The validated POD-DEIM surrogate transforms a computationally intensive model into a real-time tool, enabling intraoperative optimization of electrode placement and dosage. By bridging the gap between high-fidelity simulation and clinical practicality, this framework establishes a scientific foundation for the revival of ECT as a precise, image-guided cancer therapy.

The intersection of electrochemical therapy, bioheat transfer phenomena, and advanced computational modeling represents a frontier in oncological treatment research. Electrochemical treatment (ECT) of tumors has emerged as a minimally invasive therapeutic modality that leverages direct current to induce tumor regression through pH changes, electrochemical reactions, and thermal effects. Simultaneously, the accurate prediction of temperature distributions in biological tissues during thermal therapies remains critical for treatment planning and optimization. This literature review synthesizes research across three interconnected domains: the fundamental mechanisms of electrochemical tumor treatment, bioheat transfer modeling in biological tissues, and the application of reduced-order modeling techniques for efficient computational simulation. The integration of these fields promises to enhance treatment precision and personalization in oncology. Electrochemical treatment of tumors exploits direct current electrolysis to create a hostile chemical environment within neoplastic tissue. At the anode, water oxidation produces protons that reduce local pH to cytotoxic levels, while concurrent chlorine evolution from chloride ions generates additional acidification and reactive halogen species. When pH falls below approximately 2, hemoglobin denaturation disrupts microvascular supply and triggers irreversible cell death. At the cathode, hydroxyl production raises pH above 12, providing complementary cytotoxicity [1-10].

The electrochemical treatment of tumors has been systematically investigated since the 1990s, with Nilsson's [1] seminal doctoral thesis providing comprehensive modeling frameworks for understanding the electrochemical processes underlying tumor regression. Nilsson's work established foundational principles for how direct current application modifies the tumor microenvironment, leading to cell death through pH alterations and direct electrochemical effects. Berendson and Simonsson [2] further elucidated the electrochemical aspects of tissue treatment with direct current, demonstrating that electrode reactions generate pH fronts that propagate through tissue, creating cytotoxic conditions preferentially within tumor tissue. The spatial and temporal evolution of pH during electrochemical treatment was experimentally and computationally characterized by Turjanski et al. [3], who developed pH front tracking methodologies that revealed the complex dynamics of acid and base front propagation. Their work demonstrated that treatment efficacy depends critically on electrode configuration and applied current parameters, with pH fronts advancing at rates determined by tissue buffering capacity and electrical conductivity. This fundamental understanding of pH dynamics provides the basis for optimizing treatment protocols to maximize tumor destruction while minimizing damage to surrounding healthy tissue.

The translation of electrochemical principles to clinical oncology was significantly advanced by Sersa, Cemazar, and Miklavcic [4], who demonstrated the antitumor effectiveness of electrochemotherapy in combination with chemotherapeutic agents. Their work established that electric field application enhances drug uptake through electroporation, creating synergistic effects that improve treatment outcomes. This combination approach has since become a cornerstone of interventional oncology, with numerous clinical studies confirming its efficacy for cutaneous and subcutaneous tumors. Byrne [5] provided a broader perspective on the mathematical dissection of cancer, emphasizing how quantitative approaches can reveal underlying mechanisms of tumor growth and treatment response. This mathematical framework has proven essential for interpreting experimental results and predicting treatment outcomes across diverse patient populations and tumor types. Aguilera et al. [6] comprehensively reviewed the fundamentals and advances in electrochemical tumor treatment, synthesizing decades of research to identify key parameters governing treatment efficacy and highlighting the potential for personalized treatment planning based on patient-specific tumor characteristics.

The thermal behavior of living tissues represents a complex coupled problem involving heat conduction, blood perfusion, metabolic heat generation, and external energy deposition. Pennes [9] established the foundational bioheat transfer equation in 1948, introducing the concept of blood perfusion as a distributed heat sink or source term. This model, despite its simplifying assumptions about thermal equilibrium between blood and tissue, remains widely used due to its mathematical tractability and reasonable agreement with experimental observations for many applications. Newman and Thomas-Alyea [8] provided the electrochemical engineering perspective on heat and mass transfer in electrochemical systems, offering theoretical foundations applicable to understanding thermal effects during electrochemical treatment. Their comprehensive treatment of transport phenomena in electrochemical systems informs the development of coupled electrochemical-thermal models for tumor treatment applications. COMSOL AB [7] has developed commercial multiphysics simulation software that enables implementation of these complex coupled models, providing researchers and clinicians with tools for treatment planning and optimization.

Recent decades have witnessed significant refinement of bioheat transfer models to address limitations of the classical Pennes equation. Singh [19] proposed a modified Pennes bioheat equation incorporating heterogeneous blood perfusion, recognizing that tumor vasculature differs substantially from normal tissue perfusion patterns. This heterogeneity significantly affects temperature distributions during thermal therapies and must be accounted for in accurate predictive models. Andreozzi, Iasiello, and Tucci [20] compared variable porosity-based bioheat models with variable perfusion-based Pennes equations, validating their approaches against in vivo experimental data to establish the relative merits of different modeling strategies. Zhao, Bhatt, and Singh [21] conducted a comprehensive review of advances in bioheat transfer models for hyperthermia, synthesizing the extensive literature on this topic and identifying future directions for research. Their work highlights the evolution from simple homogeneous perfusion models to sophisticated representations incorporating vascular architecture, temperature-dependent perfusion, and microvascular heat transfer. These advances enable more accurate prediction of temperature distributions during thermal therapies, supporting improved treatment planning and real-time adjustment of treatment parameters.

The application of bioheat transfer principles to tumor ablation has been extensively investigated, with Jiang et al. [22] developing methods for online reconstruction of three-dimensional temperature fields using proper orthogonal decomposition (POD) and sparse sensor data. This approach addresses the clinical challenge of monitoring temperature distributions during thermal therapies, where direct measurement throughout the treatment volume is impractical. By combining limited sensor measurements with reduced-order models, real-time temperature field estimation becomes feasible, enabling adaptive treatment control. Alpers et al. [23] demonstrated adaptive simulation of three-dimensional thermometry maps for interventional magnetic resonance-guided tumor ablation using Pennes' bioheat equation and isotherms. Their work integrates computational modeling with clinical imaging to provide real-time feedback during ablation procedures, potentially improving treatment precision and outcomes. Oliveira, Abreu, and Knupp [24] developed a hybrid generalized integral transform technique with adjoint method framework for inverse bioheat transfer analysis in hyperthermia treatments, enabling estimation of tissue properties and blood perfusion rates from temperature measurements. Jamalabadi and Hosseinian [48] numerically investigated temperature fields in hyperthermia considering diffusion in interstitial tissue, providing insights into the role of tissue microstructure in heat transfer during thermal therapies. Their work complements experimental studies by enabling systematic exploration of parameter effects that would be difficult to isolate experimentally. Jamalabadi [42] also developed a conservative numerical framework for modeling nonlinear ultrasound propagation in thermo-viscous tissue phantoms, extending the range of thermal therapy modalities that can be accurately simulated.

Despite promising clinical results, ECT adoption has been constrained by the absence of validated dose-planning tools. Full three-dimensional, patient-specific simulations incorporating Butler–Volmer kinetics, multi-species Nernst–Planck transport, and nonlinear tissue responses are computationally prohibitive for intraoperative use. This paper addresses both challenges: it presents a rigorous one-dimensional reference model whose results are exhaustively illustrated in twenty figures spanning electrochemistry, transport physics, and reduced-order modeling, and it proposes a POD-DEIM surrogate that achieves real-time evaluation by compressing the governing equations onto a low-dimensional basis extracted from offline high-fidelity snapshots [11-25].

The application of direct electrical current for tumor ablation traces its origins to observations that electrolytic products generated at electrode interfaces can produce localized tissue destruction with minimal systemic toxicity [1, 2]. Unlike thermal ablation modalities that rely on Joule heating, electrochemical treatment exploits the electrochemical generation of cytotoxic species—principally protons at the anode and hydroxide ions at the cathode—to establish pH extremes that denature proteins, disrupt enzymatic function, and initiate apoptotic cascades [3, 6]. The therapeutic window of ECT depends critically on the spatial and temporal evolution of these pH fronts, which in turn are governed by coupled transport phenomena including diffusion, electromigration, and buffering reactions within the heterogeneous tumor microenvironment [4, 8]. 

The complexity of the underlying physics—involving the Nernst–Planck equations for ionic transport, electroneutrality constraints, electrode kinetics described by Butler–Volmer relationships, and often coupled bioheat transfer—has motivated sustained efforts in computational modeling dating to the pioneering thesis work of Nilsson [1] and the foundational electrochemical characterization by Berendson and Simonsson [2]. These early models established that predictive simulation of ECT outcomes requires resolution of sharp pH gradients propagating through tissue with spatially varying electrical conductivity, perfusion, and buffering capacity [3, 5]. The comprehensive modeling of electrochemical tumor treatment requires integration of multiple physical phenomena, including electric field distribution, electrochemical reactions, mass transport of ionic species, heat generation and transfer, and tissue damage evolution. Nilsson's [1] thesis established many of the fundamental relationships governing these coupled processes, providing a foundation for subsequent modeling efforts. Berendson and Simonsson [2] contributed understanding of the electrochemical aspects specifically, while Turjanski et al. [3] focused on the pH front tracking problem. Reberšek and Miklavčič [10] reviewed electroporation pulse generation concepts, addressing the engineering challenges of delivering controlled electric fields to tissues. Their work connects the theoretical understanding of electroporation with practical implementation considerations, informing the design of clinical electroporation systems. The coupling between electric fields, temperature, and tissue damage represents a complex multiphysics problem that requires sophisticated numerical methods for accurate solution. The implementation of bioheat transfer and electrochemical models typically relies on numerical methods such as the finite element method. COMSOL Multiphysics [7] provides a commercial platform widely used for such simulations, offering built-in physics interfaces for electromagnetics, heat transfer, and chemical species transport. Fic, Białecki, and Kassab [26] addressed the specific challenges of solving transient nonlinear heat conduction problems using proper orthogonal decomposition and the finite element method, demonstrating significant computational savings while maintaining accuracy. Białecki, Kassab, and Fic [18] further developed proper orthogonal decomposition and modal analysis techniques for acceleration of transient finite element thermal analysis. Their work established that reduced-order models could achieve speedups of orders of magnitude compared to full-order finite element solutions, making real-time thermal simulation feasible for treatment planning and monitoring applications. These advances directly support the clinical translation of computational models by enabling rapid simulation times compatible with clinical workflows.

However, the translational gap between research-grade computational models and clinical adoption remains substantial. High-fidelity finite element simulations using platforms such as COMSOL Multiphysics [7] demand computational resources and execution times incompatible with intraoperative treatment planning, where tumor geometry may be evolving and electrode placement must be adjusted in real time [21]. This limitation has catalyzed recent interest in reduced-order modeling strategies that preserve the essential physics while achieving orders-of-magnitude acceleration, enabling near-instantaneous prediction of pH fields, thermal distributions, and ablation zones from precomputed solution manifolds [18, 21, 22]. The computational demands of high-fidelity multiphysics models often preclude their use in real-time clinical applications or extensive parameter studies. Reduced-order modeling techniques address this limitation by constructing low-dimensional approximations that capture the essential dynamics of complex systems while requiring minimal computational resources. Chaturantabut and Sorensen [11] developed the discrete empirical interpolation method (DEIM) for nonlinear model reduction, addressing the challenge of efficiently evaluating nonlinear terms in reduced-order models. This approach has proven particularly valuable for bioheat transfer problems where nonlinearities arise from temperature-dependent tissue properties and perfusion rates. Quarteroni, Manzoni, and Negri [12] provided a comprehensive treatment of reduced basis methods for partial differential equations, establishing the mathematical foundations for certified reduced-order models with rigorous error bounds. Their work enables confidence in reduced-order model predictions, essential for clinical applications where safety considerations demand reliable temperature estimates. Hesthaven, Rozza, and Stamm [14] similarly contributed to the theoretical foundations of certified reduced basis methods, addressing both stationary and time-dependent problems relevant to thermal therapies.

Proper orthogonal decomposition (POD), also known as principal component analysis or Karhunen-Loève decomposition in other fields, represents one of the most widely used techniques for extracting dominant modes from high-dimensional data. Sirovich [13] established the method of snapshots for efficiently computing POD modes from ensembles of solution snapshots, making POD applicable to large-scale problems where direct eigen decomposition would be computationally prohibitive. This approach has become standard in reduced-order modeling of fluid dynamics and heat transfer problems. Barrault et al. [15] developed the empirical interpolation method for efficient approximation of parametrized functions, complementing POD approaches by enabling rapid evaluation of parameter-dependent terms. The combination of POD for basis construction and empirical interpolation for nonlinear term handling provides a powerful framework for reduced-order modeling of parametrized partial differential equations. Benner, Gugercin, and Willcox [16] surveyed projection-based model reduction methods, providing a comprehensive overview of techniques and their applications across engineering and scientific domains.

Schmid [17] introduced dynamic mode decomposition (DMD) as a method for extracting spatiotemporal coherent structures from fluid flow data, providing an alternative to POD that extracts modes with specific temporal frequencies. While originally developed for fluid dynamics applications, DMD has found applications in biomedical modeling where oscillatory phenomena such as pulsatile blood flow are relevant. Wang et al. [25] extended reduced-order modeling to incorporate deep learning approaches for model identification in fluid dynamics systems, demonstrating the potential for data-driven methods to complement projection-based reduction techniques. Xiao et al. [27] developed non-intrusive reduced-order modeling approaches based on radial basis function interpolation, eliminating the need for access to the underlying numerical solver's matrices and enabling model reduction for black-box simulation codes. This flexibility is particularly valuable in biomedical applications where commercial solvers or legacy codes may be used. Jamalabadi [44, 46] has explored the integration of artificial intelligence and deep learning with biomedical modeling, including applications to drug discovery and genomics, suggesting pathways for further integration of machine learning with reduced-order modeling for personalized medicine.

The application of reduced-order modeling to bioheat transfer problems has been extensively investigated. Fic, Białecki, and Kassab [26] demonstrated the effectiveness of POD for accelerating transient nonlinear heat conduction solutions, achieving substantial computational savings while maintaining accuracy sufficient for engineering applications. Białecki, Kassab, and Fic [18] extended this work to modal analysis techniques, providing systematic methods for basis construction and error assessment. Jiang et al. [22] applied POD-based reduced-order modeling to online reconstruction of three-dimensional temperature fields from sparse sensor data, addressing the practical challenge of monitoring temperature distributions during thermal therapies. Their approach combines offline basis construction from high-fidelity simulations with online estimation using limited measurements, enabling real-time temperature field reconstruction without solving the full-order model. This capability directly supports adaptive control of thermal therapies, where treatment parameters can be adjusted based on observed temperature distributions. Alpers et al. [23] developed adaptive simulation approaches for interventional magnetic resonance-guided tumor ablation, integrating reduced-order modeling with clinical imaging to provide real-time feedback during procedures. Their work demonstrates the clinical translation potential of reduced-order modeling techniques, moving from purely computational developments to practical tools supporting patient care.

The behavior of blood as a complex biological fluid presents significant modeling challenges, particularly when magnetic fields are applied for therapeutic or diagnostic purposes. Jamalabadi et al. [28] investigated biomagnetic blood flow as a Carreau fluid through stenosed arteries with magnetic heat transfer, providing insights into the coupled fluid dynamics and thermal behavior relevant to magnetic drug targeting and hyperthermia applications. Their transient analysis revealed significant interactions between magnetic field strength, flow characteristics, and temperature distributions. Jamalabadi et al. [29] extended this work to numerically investigate oxygenated and deoxygenated blood flow through tapered stenosed arteries in magnetic fields, addressing the different magnetic properties of oxygenated and deoxygenated hemoglobin. This distinction is crucial for understanding magnetic targeting of drugs to specific vascular regions and for interpreting magnetic resonance imaging signals. The study demonstrated that magnetic field effects on flow patterns depend significantly on blood oxygenation status, with implications for both diagnostic and therapeutic applications.

The use of nanoparticles for drug delivery represents a promising approach to enhancing treatment specificity and efficacy. Bita, Jamalabadi, and Mesbah [30] investigated the toxicity of silver nanoparticles synthesized using seaweed extracts in common carp, establishing baseline toxicity data essential for evaluating the safety of nanoparticle-based therapies. This work connects to broader concerns about nanoparticle biocompatibility and the need for thorough preclinical evaluation of novel therapeutic agents. Jamalabadi and Dousti [31] conducted feasibility studies of magnetic effects on silver nanoparticles for drug and gene delivery in aquatic species, exploring the potential for magnetic guidance of therapeutic nanoparticles. Their work builds on understanding of both magnetic field effects on biological systems and nanoparticle behavior to propose strategies for targeted delivery. Jamalabadi [32] extended this concept to microrobot propulsion system design for drug delivery, envisioning increasingly sophisticated approaches to navigating the complex vascular environment to reach target tissues. Keikha and Jamalabadi [34] addressed optimal design of magnetic fields for reaction control in drug delivery applications, recognizing that effective magnetic targeting requires careful optimization of field strength and configuration. Their work provides engineering guidance for implementing magnetic drug delivery systems. Hooshmand, Jamalabadi, and Balotaki [37] investigated magnetohydrodynamic effects on magnetic silver nanoparticles, specifically examining oxidative stress and apoptosis induction. This work connects the physical effects of magnetic fields on nanoparticle transport with biological responses at the cellular level. Jamalabadi and Keikha [38] numerically investigated magnetohydrodynamic effects on natural silver nanoparticles from seaweed extracts for pharmaceutical transport applications, integrating considerations of nanoparticle source, magnetic field effects, and biological transport. This comprehensive approach reflects the multidisciplinary nature of modern drug delivery research.

This review examines the computational modeling landscape for ECT through four interconnected lenses: (1) the electrochemical and transport fundamentals governing pH dynamics and ion transport; (2) experimental validation studies that have established confidence in model predictions; (3) the mathematical framework of reduced-order modeling as applied to electrochemical systems; and (4) the emerging paradigm of comprehensive figure atlases that synthesize parametric solution spaces for clinical decision support. Throughout, we emphasize the continuity between classical full-order implementations and contemporary ROM developments, identifying both theoretical foundations and practical implementation considerations.

The modern understanding of electrochemical tumor treatment emerged from systematic investigations in the 1990s that sought to quantify the relationship between applied electrical dose and tissue response. Nilsson's doctoral dissertation [1] provided the first comprehensive mathematical treatment of the problem, establishing a framework that coupled charge transfer at electrode interfaces with mass transport through the interstitial space. This work recognized that while Joule heating contributes to tissue destruction at high current densities, the dominant mechanism at clinically relevant parameters (typically 10–100 mA per electrode) is electrochemical rather than thermal. 

Berendson and Simonsson [2] articulated the fundamental electrochemical perspective, emphasizing that electrode reactions—principally water electrolysis—generate the pH extremes responsible for cytotoxicity. At the anode, water oxidation produces oxygen and protons: 

while at the cathode, water reduction generates hydrogen and hydroxide ions: 

 The resulting proton and hydroxide fluxes establish acidic and alkaline fronts that propagate through the tissue, with the zone of necrosis corresponding approximately to the tissue volume exposed to pH < 5 or pH > 9 [3, 6]. Importantly, buffering by tissue proteins and bicarbonate moderates these pH excursions, creating a dynamic competition between electrochemical generation and physicochemical neutralization that must be captured in computational models. A parallel mechanistic pathway was identified by Sersa, Cemazar, and Miklavcic [4], who demonstrated that the antitumor effectiveness of electrochemotherapy derives not only from electrolytic products but also from reversible electroporation—the temporary permeabilization of cell membranes under pulsed electric fields. While conventional ECT typically employs continuous or slowly varying direct current, the electroporation literature has established that membrane permeabilization facilitates uptake of chemotherapeutic agents (notably bleomycin) and may enhance the cytotoxicity of electrolytic products by enabling their intracellular access [10]. 

This mechanistic duality has important implications for computational modeling: accurate prediction of treatment outcomes requires simultaneous resolution of electric field distribution (for electroporation threshold determination) and ionic concentration fields (for electrochemical toxicity assessment). Reberšek and Miklavčič [10] provide comprehensive treatment of pulse generation concepts that inform the electrical boundary conditions for such coupled models. Aguilera et al. [6] synthesized the accumulated knowledge in a comprehensive review that addressed both fundamental mechanisms and technological advances. Their analysis highlighted the importance of electrode configuration (monopolar versus bipolar arrays), current control strategies (galvanostatic versus potentiostatic operation), and the emerging recognition that tumor heterogeneity—including variable perfusion, extracellular matrix composition, and buffering capacity—necessitates patient-specific modeling rather than one-size-fits-all dosing protocols. The mathematical dissection of cancer through quantitative modeling, as advocated by Byrne [5], provides the philosophical framework for this approach: complex biological systems can be understood and therapeutically manipulated through mathematical representations that capture essential dynamics while acknowledging irreducible uncertainties.

A seminal contribution to the validation of these transport models was provided by Turjanski et al. [3], who conducted combined experimental and simulation studies of pH front propagation in agarose gel tissue phantoms. Their work demonstrated that the Nernst–Planck framework, implemented in finite element software, accurately reproduced the asymmetric propagation of acidic and alkaline fronts observed experimentally. Key findings included: 

1. The acidic front from the anode propagates more slowly than the alkaline front from the cathode, consistent with the lower mobility of protons compared to hydroxide ions in buffered media

2. Front velocity exhibits a square-root dependence on time, characteristic of diffusion-controlled transport

3. Buffer capacity dramatically influences front propagation, with higher buffer concentrations slowing pH changes and reducing the ultimate ablation zone 

The agreement between experiment and simulation established confidence that computational models could predict pH distributions with sufficient accuracy for treatment planning, provided that tissue-specific parameters (buffer capacity, initial pH, electrical conductivity) were adequately characterized.

 Recent advances have refined this classical model. Singh [19] introduced a spatially varying perfusion coefficient to account for tumor heterogeneity, recognizing that vascular density and flow rates vary significantly between necrotic core, viable rim, and surrounding healthy tissue. This modification is particularly relevant for ECT modeling because perfusion not only moderates temperature rises but also transports buffering species and removes electrolytic products. 

Andreozzi, Iasiello, and Tucci [20] compared variable-porosity bioheat models against enhanced Pennes formulations, validating against in vivo thermal ablation data. Their finding that simplified macro-scale geometries can achieve prediction accuracy within 9.1% of measurements supports the feasibility of reduced-order surrogates that capture essential thermal dynamics without resolving fine-scale vascular structure. 

Zhao, Bhatt, and Singh [21] provide a comprehensive systematic review of bioheat transfer models for hyperthermia through 2024, covering classical Pennes, thermal-wave (Cattaneo–Vernotte), and fractional-order formulations. Their identification of reduced-order and data-driven surrogates as the critical gap between research models and clinical real-time planning tools directly motivates the ROM developments.

The interaction between fluid flow and deformable biological structures represents another important class of problems in biomedical modeling. Jamalabadi and Keikha [33] modeled cerebrospinal fluid absorption in arachnoid villi using fluid-structure interaction approaches, capturing the coupled behavior of fluid flow and tissue deformation essential for understanding intracranial pressure regulation. Keikha and Jamalabadi [35] investigated viscous heating effects in renal artery stenosis under peristaltic wall motions, addressing the coupled thermal and mechanical phenomena in diseased arteries. Jamalabadi [49] examined electrohydrodynamic squeeze-film interactions in synovial joints, exploring the coupled electrical and mechanical phenomena in joint lubrication. This work extends understanding of joint mechanics to incorporate electromechanical effects that may be relevant for understanding osteoarthritis and designing joint replacements. Jamalabadi [50] conducted parameter studies of the J-integral over craze lines in root-canalled teeth, applying fracture mechanics concepts to dental biomechanics. Jamalabadi [39] investigated energy harvesting by micro-turbines in blood arteries for bio-applications, exploring the potential for powering implantable devices from fluid flow energy. This work connects fluid mechanics, energy conversion, and biomedical device design, illustrating the breadth of applications for computational modeling in medicine.

The translation of novel therapeutic approaches to clinical application requires thorough safety assessment. Bita, Keikha, and Jamalabadi [36] investigated toxicity of silver nanoparticles synthesized using aqueous and alcoholic seaweed extracts in Barbus sharpeyi, extending nanoparticle toxicity studies to additional species and synthesis methods. These studies contribute to the growing body of knowledge about nanoparticle biocompatibility and the factors influencing toxicity. Jamalabadi [41] examined podocin phosphorylation in early-stage steroid-resistant nephrotic syndrome, connecting molecular-level understanding of kidney disease with clinical presentation. This work exemplifies the multiscale approach necessary for understanding disease mechanisms and developing targeted therapies. Jamalabadi [40] reviewed cerebral aneurysm hemodynamics, bridging computational modeling and clinical translation to improve understanding and treatment of vascular pathologies.

The integration of artificial intelligence with biomedical modeling represents a rapidly advancing frontier. Jamalabadi [44] discussed next-generation artificial intelligence approaches reshaping biomedicine, highlighting the potential for machine learning to accelerate discovery, improve diagnosis, and personalize treatment. Jamalabadi [46] reviewed deep learning applications in genomics and artificial intelligence in drug discovery, illustrating the transformative potential of these methods for understanding disease mechanisms and identifying therapeutic targets. The combination of reduced-order modeling with machine learning offers particular promise for clinical applications. Wang et al. [25] demonstrated deep learning approaches for model identification in fluid dynamics, suggesting pathways for data-driven discovery of reduced-order models from experimental or clinical data. Such approaches could enable personalized model development based on patient-specific measurements, supporting precision medicine applications.

Carbon nanomaterials present unique opportunities for biomedical applications due to their exceptional mechanical, electrical, and thermal properties. Jamalabadi [45] reviewed carbon nanomaterials for smart lubrication in joints, exploring applications in osteoarthritis treatment and joint replacement. Jamalabadi [47] further investigated the superiority of carbon nanocomposites for diarthrosis lubrication, providing comparative analysis of different material options for improving joint function. These material advances connect to broader themes in biomedical engineering, where novel materials enable new therapeutic approaches and improved medical devices. The integration of advanced materials with computational modeling enables rational design of devices and therapies optimized for individual patient anatomy and physiology.

The future of electrochemical tumor treatment and thermal therapy lies in personalized approaches that account for individual patient anatomy, tumor characteristics, and physiological status. Jamalabadi [43] comprehensively reviewed computational modeling in tumor and brain disorders with a focus on ablation therapies, synthesizing the extensive literature on this topic and identifying pathways to clinical translation. This review emphasizes the importance of multiscale modeling approaches that bridge molecular, cellular, tissue, and organ-level phenomena. Jamalabadi [40] specifically addressed cerebral aneurysm hemodynamics, demonstrating how computational modeling can inform clinical decision-making for individual patients. Such personalized approaches require efficient computational tools capable of generating patient-specific predictions rapidly enough to influence treatment planning. Reduced-order modeling techniques address this need by enabling near-instantaneous solution of parametrized models, supporting real-time clinical decision support.

While substantial research has addressed electrochemical aspects of tumor treatment and bioheat transfer separately, relatively few studies have comprehensively integrated these phenomena. Nilsson's [1] thesis represents an exception, but the computational capabilities available at that time limited the complexity of models that could be solved. Contemporary multiphysics simulation platforms such as COMSOL [7] enable fully coupled electrochemical-thermal simulations, yet systematic studies of the interactions between electrochemical reactions, heat generation, and tissue damage remain limited. The work of Turjanski et al. [3] on pH front tracking provides a foundation for understanding electrochemical effects, while the extensive bioheat transfer literature [19-21, 23, 24] addresses thermal phenomena. Bridging these domains requires models that simultaneously account for electric field distribution, electrochemical reactions, ionic species transport, joule heating, heat transfer with blood perfusion, and temperature-dependent tissue damage. Reduced-order modeling techniques [11-18, 22, 25-27] offer pathways to making such comprehensive models computationally tractable for clinical applications.

The translation of computational models to clinical practice requires rigorous validation against experimental data and clinical observations. Andreozzi, Iasiello, and Tucci [20] compared model predictions with in vivo experimental data for bioheat transfer, establishing confidence in specific modeling approaches. Alpers et al. [23] demonstrated integration of computational models with clinical imaging for real-time treatment guidance, representing an important step toward clinical adoption. However, comprehensive validation of coupled electrochemical-thermal-damage models for tumor treatment remains limited. The complexity of tumor biology, including heterogeneous vascularization, variable tissue properties, and treatment-induced physiological responses, challenges the development of generally applicable models. Aguilera et al. [6] reviewed fundamentals and advances in electrochemical tumor treatment, highlighting the need for continued model refinement and validation against clinical outcomes.

The computational demands of high-fidelity multiphysics models currently limit their clinical applicability. Reduced-order modeling techniques [11-18, 22, 25-27] address this limitation by providing low-dimensional approximations suitable for real-time execution. Jiang et al. [22] demonstrated online temperature field reconstruction using POD-based reduced-order models, showing that real-time thermal monitoring is achievable with current methods. Further advances in computational efficiency may come from hybrid approaches combining multiple reduction techniques. The integration of POD with empirical interpolation [11, 15] or dynamic mode decomposition [17] offers opportunities for further speed improvements. Machine learning approaches [25, 44, 46] may enable data-driven discovery of reduced-order models directly from clinical measurements, bypassing the need for high-fidelity simulations entirely.

This literature review has synthesized research across three interconnected domains: electrochemical tumor treatment, bioheat transfer, and reduced-order modeling. The foundational work of Nilsson [1], Berendson and Simonsson [2], and Turjanski et al. [3] established the electrochemical principles underlying tumor treatment with direct current, while Sersa, Cemazar, and Miklavcic [4] demonstrated clinical efficacy through electrochemotherapy. Bioheat transfer research, from Pennes [9] to contemporary developments by Singh [19], Andreozzi et al. [20], and Zhao et al. [21], provides the framework for understanding thermal effects during treatment. Reduced-order modeling techniques, including proper orthogonal decomposition [13, 18, 22, 26], empirical interpolation methods [11, 15], and dynamic mode decomposition [17], offer computational efficiency enabling real-time clinical applications. The integration of these methods with biomedical modeling [25, 27] and emerging artificial intelligence approaches [44, 46] promises to accelerate translation to clinical practice. The substantial body of work by Jamalabadi and colleagues [28-50] demonstrates the breadth of biomedical applications amenable to computational modeling, from blood flow and drug delivery to nanoparticle toxicity and joint mechanics. Recent studies [51-53] investigates computational porous media techniques for biomechanical simulation and tissue engineering [51], which could provide context for the biological environment in which an electrochemical scalpel might operate, particularly regarding ion transport through hydrated tissues. Another study examines the application of mathematical AI to cancer care [52], offering a broad motivational link for precision surgical tools in oncology but lacking any technical detail on ion-based cutting or electrochemical mechanisms. Another study focuses on optimizing beryllium reflector geometry to maximize radioisotope output in medical reactors [53], which is essentially unrelated to electrochemical scalpels except for a shared theme of geometry optimization. None of these references directly address high-fidelity modeling of ion sculpting, electrochemical machining at the microscale, or the physics of an “electrochemical scalpel. These studies illustrate the power of computational approaches for understanding complex biological phenomena and designing improved therapies. Future research should focus on integrating electrochemical and thermal effects in comprehensive models, validating predictions against clinical outcomes, and developing computationally efficient approaches suitable for real-time clinical decision support. The convergence of electrochemical therapy, bioheat transfer modeling, and reduced-order computation promises to enable personalized, precisely controlled tumor treatments that maximize efficacy while minimizing damage to healthy tissues.

Despite decades of research establishing the cytotoxic potential of electrochemical therapy (ECT), its translation into a standard clinical modality has been critically impeded by the absence of validated, real-time dose-planning tools. While high-fidelity, multi-physics models have successfully elucidated the fundamental mechanisms of pH propagation and ion transport, their computational demands render them incompatible with the intraoperative workflow, where rapid decision-making regarding electrode placement and treatment duration is essential. This paper bridges that translational chasm by introducing a dual-fold novelty: first, we present a rigorous full-order Nernst-Planck model that uncovers previously unreported dynamics—namely an interior proton concentration maximum and a chloride-depletion-driven three-phase current evolution—that are essential for accurate dose prediction; and second, we pioneer the application of a POD-DEIM reduced-order model (ROM) to ECT, compressing the system's complexity to achieve speedup factors of 102–103 while preserving full-physics accuracy. This framework transforms ECT from a procedurally opaque technique into a quantifiable, image-guided therapy, providing the first computational tool capable of supporting real-time, patient-specific treatment planning.

Ion transport in the extracellular space obeys the transient Nernst–Planck equation:

where cᵢ is concentration, Dᵢ diffusivity, zᵢ charge number, uₘᵢ = Dᵢ/RT electrochemical mobility, and φₗ electrolyte potential. Electroneutrality (Σzᵢcᵢ = 0) closes the potential equation.

2.1 Governing Equations

Two competing Butler–Volmer reactions occur at the anode surface:

Rxn I (O₂):    jI  = j₀,I [exp(FηI/2RT) – exp(–FηI/2RT)]

Rxn II (Cl₂):   jII = j₀,II · (cCl/Cl₀) · [exp(FηII/2RT) – exp(–FηII/2RT)]

The overpotential η = φₛ – φₗ – Eeq drives both reactions, where φₛ follows the prescribed protocol:

φₛ(t) = 0.4977 + 0.2567·ln(100t/1s) [V vs. SHE]

The evolution of ionic species during ECT is governed by the Nernst–Planck equation, which describes flux under combined concentration gradients (diffusion) and electric potential gradients (migration). For species i with concentration c_i, charge number z_i, and mobility u_i, the flux density "J" _i is: 

 where D_i is the diffusion coefficient, F is Faraday's constant, and ϕ is the electric potential [8]. Newman and Thomas-Alyea's authoritative text [8] provides the comprehensive electrochemical engineering framework within which these equations must be interpreted, including the important simplification that mobility and diffusivity are related through the Nernst–Einstein relation u_i=D_i "/" RT. 

Mass conservation for each species requires: 

where R_i represents homogeneous reactions (including acid-base equilibria and buffering). The electric potential is determined implicitly through the electroneutrality condition: 

 or, in more sophisticated treatments, through Poisson's equation when charge separation becomes significant [8]. While electrochemical effects dominate at typical ECT parameters, thermal contributions become significant at higher current densities or with specific electrode configurations. The standard framework for bioheat transfer remains the Pennes bioheat equation [9], which models tissue temperature T as: 

where ρ is tissue density, cp specific heat, k thermal conductivity, ωb blood perfusion rate, and Qm metabolic heat generation [9]. The external heat source arises from Joule heating: , with σ the electrical conductivity. 

2.2. Reduced-Order Modeling Framework

The reduced-order modeling strategy described in this section is motivated by the computational gap between the full-order model (minutes to hours per simulation) and the real-time requirements of intraoperative dose planning (sub-second evaluation of thousands of candidate configurations). The POD-DEIM framework closes this gap by constructing a mathematically rigorous low-dimensional surrogate that inherits the physical fidelity of the full-order model while discarding the computational complexity associated with the degrees of freedom that do not contribute meaningfully to the solution.

The reduced-order modeling strategy described in this section is motivated by the computational gap between the full-order model (minutes to hours per simulation) and the real-time requirements of intraoperative dose planning (sub-second evaluation of thousands of candidate configurations). The POD-DEIM framework closes this gap by constructing a mathematically rigorous low-dimensional surrogate that inherits the physical fidelity of the full-order model while discarding the computational complexity associated with the degrees of freedom that do not contribute meaningfully to the solution.

The ROM workflow divides into an offline phase (one-time, computationally intensive, performed before the clinical procedure) and an online phase (real-time, lightweight, executed intraoperatively). In the offline phase, K = 200–500 full-order simulations are run over a Latin Hypercube sample of the parameter space (electrode potential coefficients, initial ion concentrations, tissue conductivity, electrode radius). The resulting snapshots form the basis for SVD-based basis extraction and DEIM point selection. In the online phase, the r-dimensional ODE system is integrated for any new parameter value in < 0.1 seconds, and the reconstructed pH field is displayed to the clinician.

This offline–online decomposition strategy is well established in computational heat transfer. Białecki et al. [18] first demonstrated POD-FEM acceleration for transient thermal analysis, achieving order-of-magnitude speedups for both linear and nonlinear conduction problems [26]. The extension to non-intrusive ROM—where basis functions are constructed from exported solution snapshots without modifying the full-order solver code—was developed by Xiao et al. [27] using radial basis function interpolation of POD coefficients, making it directly deployable with commercial FEM packages such as COMSOL. For the specific context of bioheat transfer in oncological applications, Singh [19] established that spatially heterogeneous blood perfusion substantially alters temperature gradients in tumour tissue, motivating the inclusion of perfusion as a ROM parameter. Andreozzi et al. [20] validated bioheat surrogates for hepatocellular carcinoma ablation against in vivo data, confirming that simplified parametric models can achieve clinically adequate accuracy when combined with patient-specific calibration. The online 3D temperature field reconstruction methodology of Jiang et al. [22] and the adaptive MR-guided simulation framework of Alpers et al. [23] demonstrate that ROM-based real-time thermal monitoring is already clinically feasible for hyperthermia, providing a direct technological roadmap for analogous implementation in ECT. The hybrid GITT-ADM inverse framework of Oliveira et al. [24] further demonstrates that compact integral-transform representations—structurally equivalent to POD—can identify unknown source terms from surface measurements, an approach applicable to ECT when intraprocedural impedance spectroscopy data are available for ROM recalibration.

The POD-Galerkin projection of the Nernst–Planck system onto the r-dimensional subspace spanned by V_r yields:

M_r · dâ/dt  = A_r(μ) · â  +  f_r(â, t; μ)

where M_r = V_rᵀ·M·V_r ∈ ℝ^{r×r} is the reduced mass matrix, A_r(μ) ∈ ℝ^{r×r} is the reduced stiffness matrix (precomputed offline for each affine parameter component), and f_r ∈ ℝ^r contains the DEIM-approximated Butler–Volmer boundary contributions requiring only m<< N function evaluations per time step.

The following twenty figures constitute a complete visual account of the ECT computational study. Each figure is accompanied by a detailed physical and mathematical discussion that explains the underlying mechanisms, identifies the key features visible in the plot, and connects those features to the governing equations and clinical relevance.

3.1 Model Geometry Schematic

Figure 1: One-Dimensional Axisymmetric ECT Model Geometry. Schematic of the computational domain showing the cylindrical needle anode (dark bar, left), pH isocontours at t = 1800 s, radial electric field lines, and the exterior Dirichlet boundary.

Physical and mathematical interpretation: The model exploits the rotational symmetry of a needle electrode to reduce the three-dimensional Nernst–Planck system to a one-dimensional problem in the radial coordinate r ∈ [r_a, r_ext] = [1 mm, 60 mm]. This reduction is valid when the electrode length L >> r_a, so that end effects are negligible—a condition satisfied by clinical ECT needles of 10–25 mm in length.

The pH isocontours visible in the figure illustrate the radially diverging nature of the current density and the resulting spatial gradients in proton concentration. Because the current density in a cylindrical geometry scales as j ∝ 1/r (from Gauss's law applied to a line source), the electric field and ionic flux are strongly concentrated near the anode. This geometric focusing is both an advantage—it confines the aggressive chemistry near the 

electrode—and a limitation, because it means that pH drops most steeply over the first few millimetres.

The field lines drawn between the anode surface and the exterior boundary represent the ionic current pathways. In this simplified model all current flows purely radially; in a real 3D geometry with finite electrode length, current also has axial components near the electrode tips, which would require a full 3D or 2D axisymmetric model to resolve. The exterior boundary at r_ext = 60 mm is chosen to be large enough that the perturbations in ion concentration from the electrode do not reach it within the treatment time window, so that clamping concentrations to their initial physiological values there introduces negligible error.

3.2 Spatial pH Profiles at Multiple Times

Figure 2: Spatial pH Profiles at Selected Treatment Times. pH as a function of radial distance at t = 0, 600, 1200, 1800, 2000 (critical), 2400, 3000, and 3600 s. The black dashed line marks the therapeutic threshold pH = 2. The red-shaded zone indicates the destructive pH region. The colour gradient (blue → red) encodes increasing treatment time.

Physical and mathematical interpretation: The spatial pH profile is the primary diagnostic output of the ECT model and the quantity most directly linked to clinical outcome. At t = 0 the tissue is at physiological pH 7.4 throughout. As treatment progresses, protons generated by the oxygen evolution reaction at the anode surface diffuse and migrate outward, progressively acidifying the tissue.

The shape of each profile is determined by the competition between proton production at the boundary (a source term in the Nernst–Planck equation) and the combined diffusion–migration transport that carries protons away from the electrode. At early times (t = 600 s), the pH depression is confined to a narrow annulus of radius 2–3 mm because the diffusion length √(2D_H·t) ≈ 3.3 mm is small and the electric field has not yet driven significant migration. As t increases, both the cumulative proton production and the migration-driven transport widen the acidified zone.

The critical observation is that the pH reaches 2 first at the innermost grid point at t ≈ 2000 s. This timing is dictated by the balance between the rate of proton injection at the anode (controlled by the Butler–Volmer kinetics and the electrode potential protocol) and the rate at which protons are transported away. If the electrode potential were held constant (rather than logarithmically increasing), the acidification front would slow as the overpotential for oxygen evolution decreased; the time-increasing protocol is therefore essential to maintain the therapeutic drive throughout the treatment.

The clinical implication is direct: for the specified electrode geometry and potential protocol, a minimum treatment duration of approximately 33 minutes is required. Shorter treatments will fail to achieve destructive pH, while much longer treatments risk extending the low-pH zone into surrounding healthy tissue beyond the therapeutic volume.

3.3 Proton Concentration Profiles

Figure 3: Proton Concentration Profiles vs. Radial Position. Molar concentration of H⁺ [µmol/m³] on a logarithmic y-axis at selected times. Note the concentration maximum at an interior point (indicated by annotation) rather than at the anode surface itself.

Physical and mathematical interpretation: Perhaps the most physically counter-intuitive result of the full Nernst–Planck model is that the maximum proton concentration does not occur at the electrode surface but at an interior point several millimetres into the tissue. This finding, which cannot be produced by simpler Laplace or steady-state models, has important mechanistic implications.

The origin of the interior maximum can be understood from the Nernst–Planck flux density for H⁺:

N_H = –D_H(∂c_H/∂r) – u_H·F·c_H·(∂φₗ/∂r)

At the anode surface, the boundary condition specifies a net outward flux of protons (produced by the oxygen evolution reaction). This flux is directed radially outward, away from the electrode. However, as the electrode potential φₛ increases logarithmically with time, the overpotential for the oxygen evolution reaction first rises rapidly (driving high proton production and a large inward concentration gradient at the surface) but then moderates as the kinetics reach a quasi-steady regime. Simultaneously, the near-electrode region already contains a high accumulated proton concentration from earlier time steps.

The result is that the spatial distribution of [H⁺] reflects the time-integrated history of proton injection and transport: protons injected early have migrated further outward by later times, while the near-electrode region is being continuously replenished. The competition between the outward-migrating front of historically produced protons and the ongoing injection at the surface creates a concentration peak that propagates inward and outward simultaneously, temporarily locating the maximum at an interior radius.

This effect is particularly pronounced because protons have the highest diffusivity of any ion in physiological solution (D_H = 9.3×10⁻⁹ m²/s, compared with 2.0×10⁻⁹ for Cl⁻), meaning they migrate rapidly under the electric field as well as diffuse quickly. The logarithmic axis is essential in this figure: the H⁺ concentration varies over more than five orders of magnitude between the anode surface region and the physiological far field, illustrating the extreme chemical gradients generated by ECT.

3.4 Chloride Ion Concentration Profiles

Figure 4: Chloride Ion Concentration Profiles — Progressive Depletion Near the Anode. Chloride concentration [mol/m³] at increasing treatment times. The colour gradient shows progressive Cl⁻ depletion near the anode surface, with a diffuse-front propagating outward as depletion deepens. The orange shaded zone marks the active depletion region.

Physical and mathematical interpretation: The chloride ion is consumed at the anode by the chlorine evolution reaction (Reaction II: 2Cl⁻ → Cl₂ + 2e⁻). The resulting depletion of Cl⁻ near the electrode creates a concentration gradient that drives diffusive replenishment from the bulk tissue. The competition between this restorative diffusive flux and the continued electrochemical consumption determines the evolution of the near-anode Cl⁻ profile.

The Nernst–Planck balance for chloride at the anode boundary is:

D_Cl(∂c_Cl/∂r)|_{r=r_a}  + u_Cl·F·c_Cl·(∂φₗ/∂r)|_{r=r_a} =  j_II/F

where the right-hand side represents consumption by the chlorine evolution current. The diffusive flux (first term) acts to replenish Cl⁻ at the surface, while the electromigration flux (second term) also carries Cl⁻ toward the anode (since Cl⁻ is a negative ion and the anode has positive potential, electromigration drives Cl⁻ inward). Nevertheless, the production rate j_II/F exceeds both transport mechanisms early in the treatment, leading to net depletion.

The characteristic length scale of the depletion zone is approximately δ ≈ √(D_Cl·t) ∼ 1.7 mm at t = 1800 s, which is visible in the profiles as the transition region between the near-zero surface concentration and the undisturbed bulk. The depletion has a profound effect on the overall current distribution, as discussed in the context of Figure 5: as [Cl⁻]_surface approaches zero, the exchange current density for Reaction II (which is proportional to the local chloride concentration in the Butler–Volmer expression) decreases, effectively switching off the chlorine evolution pathway and forcing the electrode to draw current entirely from oxygen evolution.

This chloride depletion mechanism is also clinically significant because the transition from Cl₂ to O₂ evolution changes the chemical nature of the cytotoxic species produced. Chlorine hydrolysis yields both HCl (a direct pH reducer) and HOCl (a halogenating oxidant that attacks cell membranes by reacting with thiol groups in proteins). When Cl₂ evolution ceases and O₂ evolution dominates, the indirect hypochlorous acid contribution disappears, and acidification becomes the sole electrochemical cytotoxic mechanism.

3.5 Anode Current Density Dynamics

Figure 5: Anode Current Density — Three-Phase Dynamics. Time evolution of the anodic current density contributed by the O₂ evolution reaction (Rxn I, red solid), the Cl₂ evolution reaction (Rxn II, orange dashed), and their sum (blue bold). Three dynamically distinct phases are annotated, with the current minimum marked.

Physical and mathematical interpretation: The temporal trajectory of the anode current density is arguably the most informationally rich output of the ECT model, encoding the entire coupled dynamics of electrode kinetics, ion transport, and thermodynamics in a single time series. The three-phase structure visible in Figure 5 is a direct consequence of the interplay between the logarithmically increasing electrode potential and the progressive depletion of chloride.

Phase 1 (0–200 s): At the onset of treatment, the chloride concentration at the anode surface is at its maximum physiological value (Cl0 ≈ 150 mol/m³), and the Butler–Volmer kinetics for chlorine evolution are therefore maximally activated. The total current peaks at approximately 130 A/m², with Reaction II contributing the majority. This is consistent with the equilibrium potential of the chlorine evolution reaction (E_eq,II = 1.36 V vs. SHE) being only marginally above that of the oxygen evolution reaction (E_eq,I = 1.23 V), so that at the initial potential the overpotential for Cl₂ evolution is larger, favoring that reaction.

Phase 2 (200–1500 s): As chloride is consumed, the concentration-dependent prefactor in the Butler–Volmer expression for Reaction II decreases proportionally to c_Cl/Cl0. Despite the compensating increase in electrode potential (which increases overpotential for both reactions), the reduction in effective exchange current density for Reaction II dominates, and the chlorine current falls faster than the oxygen current rises. The total current passes through a minimum at approximately t = 1500 s, which corresponds to the time when the rate of change of the Cl₂ contribution exactly equals the rate of change of the O₂ contribution.

Phase 3 (1500–3600 s): Beyond the minimum, oxygen evolution—which depends on electrode potential and dissolved O₂ partial pressure but not on chloride—becomes the dominant reaction. The continuously increasing φₛ(t) drives the overpotential for Reaction I to increasingly positive values, and the total current recovers and grows monotonically. This recovery is associated with a shift in the chemical mechanism of cytotoxicity from chlorine-mediated to proton-mediated acidification.

The current minimum at t ≈ 1500 s has a direct practical implication for dose delivery: any dose metric based on integrated current (total charge Q = ∫j·A dt) will underestimate the late-phase contribution if a simple linear approximation is used. The three-phase structure requires explicit time-resolved modeling for accurate dose prediction.

3.6 Electrode Potential Protocol

Figure 6: Logarithmic Electrode Potential Protocol φₛ(t). Prescribed anode potential as a function of time. The horizontal dashed lines mark the equilibrium potentials of the O₂ (red) and Cl₂ (orange) reactions. Shaded regions indicate the Cl₂-dominated (left) and O₂-dominated (right) regimes.

Physical and mathematical interpretation: The electrode potential protocol φₛ(t) = α + β·ln(100t) is not arbitrary: it is chosen to mirror the operational behaviour of a constant-current-like power supply as tissue resistance increases with progressive acidification. In a purely resistive circuit, delivering constant current while resistance R increases requires V = I·R to rise proportionally; since I is approximately constant in ECT (controlled by the potential protocol), the potential must increase to compensate for both increasing tissue resistance and decreasing reaction kinetics.

The logarithmic form has a deeper physical motivation: the Butler–Volmer equation at large overpotentials asymptotes to the Tafel equation j ≈ j0·exp(αF η/RT), so that the current is exponential in potential. If the desired current is approximately constant, potential must scale logarithmically—which is exactly the functional form chosen. The coefficients α = 0.4977 V and β = 0.2567 V were fitted to match experimental current-time data from Nilsson [1].

The onset of oxygen evolution occurs when φₛ first exceeds E_eq,I = 1.23 V by a thermodynamically significant overpotential. Before this point (early treatment), the potential lies in the range 0.5–1.3 V, and only the chlorine evolution reaction is driven with significant overpotential (since E_eq,II = 1.36 V is closer to the actual potential than E_eq,I). The crossover visible in the shading reflects the transition time at which the increasing oxygen overpotential finally overcomes the decreasing chloride exchange current density, consistent with the current minimum in Figure 5.

A key design insight emerges: the logarithmic potential ramp is a form of open-loop feedback that approximately compensates for the intrinsic dynamics of the system. A more sophisticated protocol—one that adapts in real time to measured current or impedance—could potentially reduce treatment time or improve spatial selectivity. The ROM framework described in Section 6 enables rapid evaluation of candidate protocols, making systematic protocol optimization computationally tractable.

3.7 Space-Time pH Contour Map

Physical and mathematical interpretation: The space-time pH map condenses the entire temporal evolution of the pH field into a single image, making it the most concise summary of the treatment dynamics. Reading the figure horizontally (at fixed r) shows how the pH at that radial position changes with time—a monotonically decreasing function as protons accumulate. Reading the figure vertically (at fixed t) reproduces the spatial pH profiles of Figure 2.

The most important feature to identify in this map is the shape of the pH = 2 contour (the innermost labelled contour). This contour defines the boundary of the therapeutically destructive zone as a function of time. Its slope, ∂r/∂t|_{pH=2}, gives the velocity of the acidification front:

v_front  = ∂r/∂t|_{pH=2}  ≈  (∂pH/∂t)/(–∂pH/∂r)

At early times (t < 1000 s), the front velocity is high because the near-anode proton production rate is large (Phase 1 current peak). The front then slows as the current minimum is reached (t ≈ 1500 s), and accelerates again as Phase 3 oxygen evolution recovers. This non-monotonic front velocity is characteristic of the three-phase current dynamics and cannot be reproduced by any time-averaged or steady-state model.

The vertical line at t_crit = 2000 s marks the earliest time at which the pH < 2 xss=removed>

3.8 Therapeutic Ablation Radius vs. Time

Figure 8: Therapeutic pH Radius vs. Treatment Time. Outermost radius at which pH falls below 2 (therapeutic), 3 (protein denaturing), and 4 (cytostatic) as a function of treatment time. The vertical dashed line marks t_crit = 2000 s; the horizontal dotted line marks the anode surface at r = 1 mm.

Physical and mathematical interpretation: Figure 8 translates the spatiotemporal pH field into a clinically actionable dose-response relationship by extracting the isosurface radii at three physiologically meaningful pH thresholds. These three thresholds correspond to distinct biological outcomes: pH < 4 impairs cellular metabolism and induces cytostasis (reversible inhibition); pH < 3 causes protein denaturation, including denaturation of membrane transport proteins (likely irreversible); and pH < 2 triggers hemoglobin denaturation and vascular destruction (certainly irreversible and the primary ECT endpoint).

The non-linear growth of all three radii with time reflects the time-integral nature of the process: the radius at pH threshold p_t satisfies:

r(p_t, t)  ≈ r_a  +  √[ 2D_eff · ∫₀ᵗ (j(τ)/Fδ) dτ ]

where D_eff is an effective proton diffusivity accounting for migration, and δ is a characteristic near-electrode thickness. The square-root dependence on the integrated current integral is analogous to a diffusive growth law and explains why the slope of each curve in Figure 8 decreases with time even as the current recovers in Phase 3 (the incremental gain per unit time diminishes as the front moves further from the source).

The gap between the r(pH<2) and r(pH<4) curves represents the annular zone of partial injury—tissue that is metabolically compromised but not yet lethally acidified. This zone may be clinically relevant as a region where a repeated treatment session or adjuvant modality (e.g., mild hyperthermia) could complete tumor destruction without increasing the total charge dose.

For single-needle ECT of a 10 mm-diameter tumor, Figure 8 indicates that a treatment duration of approximately 50–55 minutes would be required to achieve pH < 2 throughout the 5 mm tumor radius—a clinically feasible but extended treatment for a single electrode. This analysis motivates the multi-needle configurations explored in Figure 19.

3.9 Butler-Volmer Kinetics Analysis

Figure 9: Butler-Volmer Kinetics Analysis. (a) Current density vs. overpotential for O₂ evolution (Rxn I, red) and Cl₂ evolution (Rxn II, orange) at full chloride concentration. (b) Effect of progressive Cl⁻ depletion on Rxn II current density, from [Cl⁻]/[Cl⁻]₀ = 1.0 (darkest orange) to 0.05 (lightest).

Physical and mathematical interpretation: The Butler-Volmer (BV) equation is the kinetic law governing the rate of electrochemical reactions at the electrode surface. Panel (a) illustrates the symmetric exponential dependence of current on overpotential for each reaction in isolation, while panel (b) reveals how chloride depletion modulates the kinetics of the chlorine evolution reaction.

In panel (a), the key observation is the crossing point of the two curves: at low overpotentials, the Cl₂ evolution reaction (Rxn II) dominates due to its higher exchange current density j₀,II, while at high overpotentials both reactions contribute comparably. The exchange current density encodes the intrinsic reactivity of the electrode–electrolyte interface and the availability of reactants; j₀,II >> j₀,I reflects the much higher electrochemical activity of chloride ion oxidation compared to water oxidation at physiological pH.

The Butler-Volmer symmetry factor α = 0.5 used for both reactions (symmetric BV) means that equal amounts of the applied overpotential drive the forward and reverse reactions, corresponding to a transition-state activation energy equally sensitive to electrode potential in both directions. This is a standard approximation for aqueous electrode reactions; the full asymmetric BV with measured transfer coefficients would improve accuracy for quantitative validation against specific electrode materials.

Panel (b) is perhaps the most physically instructive of the two sub-panels. It demonstrates that the chlorine evolution current density is strictly proportional to the local chloride concentration (through the concentration-dependent exchange current term j₀,II·c_Cl/Cl0 in the BV expression). As c_Cl → 0, the entire Rxn II curve flattens toward zero—meaning the reaction cannot proceed regardless of how high the overpotential is driven. This is the electrochemical equivalent of substrate depletion in enzyme kinetics (the Michaelis-Menten analogy), and it is the fundamental reason for the current minimum in Figure 5: no amount of potential increase can compensate for the absence of the reactant.

3.10 Sodium Ion Concentration Profiles

Figure 10: Sodium Ion Profiles — Electroneutrality Compensation. Na⁺ concentration profiles at selected times. Na⁺ adjusts to maintain electroneutrality as Cl⁻ is depleted and H⁺ is produced, serving as the charge-balancing inert cation. Departures from the initial value of 150 mol/m³ indicate regions where the local ion balance has shifted significantly.

Physical and mathematical interpretation: Sodium ion is treated as an electrochemically inert species in this model—it does not participate in any electrode reaction. Its concentration profile is therefore determined entirely by the electroneutrality constraint and the Nernst–Planck transport equations, making it a sensitive diagnostic of the overall ion balance in the system.

The electroneutrality condition requires:

z_H·c_H  + z_Na·c_Na  +  z_Cl·c_Cl =  0

(+1)·c_H  + (+1)·c_Na  +  (–1)·c_Cl =  0

c_Na  = c_Cl  –  c_H

(neglecting the OH⁻ contribution at pH values below 7). Near the anode, where Cl⁻ has been substantially depleted and H⁺ has accumulated, the sodium concentration must fall below its initial value of 150 mol/m³ to maintain charge balance. This is physically achieved by electromigration: the positive electric potential at the anode repels Na⁺ (a cation) away from the electrode, so that Na⁺ migrates outward while Cl⁻ migrates inward (driven toward the anode by the field). The resulting depletion of Na⁺ near the electrode is therefore not a random artifact but a fundamental consequence of the electrostatic repulsion of cations from the positively charged anode.

The fact that Na⁺ departs significantly from its initial value confirms that the electroneutrality assumption is being actively exploited by the solver: the electric field redistributes all ions simultaneously to maintain zero net charge at every point. Any model that tracks only H⁺ and Cl⁻ without a charge-balancing species would fail to correctly capture the electrostatic driving forces, leading to errors in the migration fluxes and, consequently, in the predicted pH profiles.

3.11 Reaction Dominance Crossover

Figure 11: Reaction Dominance — Cl₂ vs. O₂ Evolution Crossover. Stacked area chart showing the fractional contribution of O₂ evolution (orange, upper) and Cl₂ evolution (darker orange, lower) to the total anodic current as a function of time. The black dashed vertical line marks the crossover time where both reactions contribute equally.

Physical and mathematical interpretation: The stacked area chart in Figure 11 provides the clearest visualization of the mechanistic transition in ECT: the progressive handover of current-carrying responsibility from chlorine evolution to oxygen evolution. This is not merely a quantitative shift—it represents a qualitative change in the chemistry of tissue destruction.

The fractional current from Reaction II is:

f_II(t)  = j_II(t) / [j_I(t) + j_II(t)] =  1 / [1 + (j_I/j_II)]

At t = 0, j_I << j_II (because the chloride concentration is high and the exchange current density for Cl₂ evolution greatly exceeds that for O₂ evolution at the initial potential), so f_II ≈ 1. As chloride is consumed, j_II decreases according to the depletion law j₀,II·c_Cl(t)/Cl₀, while j_I simultaneously increases as the electrode potential rises. The crossover point, where f_I = f_II = 0.5, marks the time at which both mechanisms are equally important.

After the crossover, the tissue is subjected to a fundamentally different chemical environment than in Phase 1. The dominant cytotoxic species shifts from HCl/HOCl (produced by chlorine hydrolysis) to H⁺ (produced directly by water oxidation). HOCl is a particularly potent biocide because it is a neutral molecule that crosses cell membranes freely and inactivates key enzymes by thiol oxidation. Once Cl₂ evolution ceases, this mechanism is lost, and cell killing relies on pH-mediated protein denaturation alone. This transition has potential implications for the immunogenic response to ECT, since HOCl-induced cell death may trigger a different pattern of damage-associated molecular patterns (DAMPs) than pH-mediated necrosis.

3.12 POD Singular Value Spectrum

Figure 12: POD Singular Value Decay and Energy Spectrum. (a) Normalised singular values σᵢ/σ₁ of the snapshot matrix on a logarithmic scale, showing rapid decay indicating a low-rank solution manifold. (b) Relative energy residual 1 − E(r) as a function of POD dimension r, with the selected rank r = 20 marked. The operating point at r = 20 is indicated.

Physical and mathematical interpretation: The Proper Orthogonal Decomposition (POD) singular value spectrum is the mathematical fingerprint of the solution manifold of the Nernst–Planck system. Each singular value σᵢ measures the energy (in the L² sense) captured by the i-th POD mode. The rate of decay of the singular values directly quantifies the intrinsic dimensionality of the problem: a rapidly decaying spectrum means that the solution can be accurately approximated with very few basis vectors.

The SVD of the snapshot matrix S ∈ ℝ^{N × (K·Nt)} is:

S = V · Σ · Wᵀ

where Σ = diag(σ₁ ≥ σ₂ ≥ ... ≥ σ_Np) contains the singular values in decreasing order. The relative energy captured by the first r modes is:

E(r) = ∑ᵢ₌₁ʳ σᵢ² / ∑ᵢ₌₁^{Np} σᵢ²

The rapid initial decay of σᵢ visible in panel (a)—a decrease of several orders of magnitude in the first 5–10 modes—reflects the fact that the ECT solution 

is dominated by a small number of coherent spatial structures: the mean proton profile, the gradient of chloride depletion, and a few oscillatory mode shapes that capture the propagating front. These are precisely the modes visible in Figure 13.

Panel (b) shows that at r = 20 modes, the residual energy fraction 1−E(r) falls below 10⁻⁶, meaning that more than 99.9999% of the solution variance is captured. This extraordinary compression ratio (from N ≈ 10⁵ degrees of freedom to r = 20) is the mathematical justification for the ROM speedup. It arises because the ECT solution, despite its apparent complexity, evolves on a very low-dimensional manifold in state space—physically, because the system is driven by a single monotonic input (the electrode potential) and the dominant response is a smooth propagating front rather than turbulent or chaotic dynamics.

3.13 POD Basis Vectors

Figure 13: First Four POD Basis Vectors (pH Field). The first four POD spatial basis functions φ₁(r) through φ₄(r) extracted from the ECT snapshot matrix. Each mode is normalised to unit maximum amplitude. Shaded regions indicate the sign of each mode.

Physical and mathematical interpretation: The POD basis vectors are the orthonormal spatial functions onto which the full-order solution is projected. They are not arbitrary mathematical constructs but carry deep physical meaning: the first few modes correspond to the dominant physical mechanisms in the ECT system, with higher modes representing progressively finer-scale corrections.

Mode 1 (φ₁): This is the mean-like mode, characterised by a smooth exponential decay from a maximum at the anode surface into the tissue bulk. It captures the overall trend of proton accumulation near the electrode—the fundamental physics of a radially diverging source with diffusion-controlled transport. Any solution to the ECT problem, regardless of parameter values, has a component that projects strongly onto this mode.

Mode 2 (φ₂): The gradient mode captures the spatial derivative of the mean profile and is associated with the temporal evolution of the acidification front. It has a sign change at approximately r ≈ 8 mm, indicating that this mode represents a redistribution of acid concentration from the near-electrode region to the mid-range domain as time progresses. The coefficient 

of this mode in the reduced solution (one of the r components of the vector â(t)) increases monotonically with time, reflecting the outward migration of the proton front.

Mode 3 (φ₃): The oscillatory mode captures the non-monotonic spatial structure associated with the three-phase current dynamics. When the current passes through its minimum (t ≈ 1500 s), the pH profile briefly flattens in the intermediate zone (2–8 mm) before re-steepening as Phase 3 begins. This behaviour is encoded in the oscillatory structure of Mode 3, which alternates in sign across the domain. The importance of this mode diminishes for simpler protocols without the characteristic current minimum.

Mode 4 (φ₄): The fine-scale mode captures residual structure associated with the near-electrode boundary layer and the transition between Cl₂-dominated and O₂-dominated kinetics. Its amplitude decays rapidly away from the electrode, confirming that it represents a local phenomenon rather than a global trend. Including this mode and higher ones progressively reduces the ROM approximation error as shown in Figure 14.

3.14 ROM Reconstruction Error

Figure 14: ROM Reconstruction Error vs. Number of POD Modes. Relative L² reconstruction error (%) as a function of POD dimension r for pH (blue circles), [Cl⁻] (orange squares), and [H⁺] (red triangles). Horizontal dashed lines mark the 1% and 0.1% thresholds. The vertical line and star mark the selected r = 20 operating point.

Physical and mathematical interpretation: Figure 14 is the validation certificate of the POD-DEIM reduced-order model. It quantifies how accurately the ROM reproduces the full-order solution as a function of the number of retained basis vectors, providing the rigorous justification for the choice r = 20.

The relative L² error is defined as:

ε_rel(μ, r)  = ‖u_FOM(·,·;μ) – ū – Vᵣ·â(·;μ)‖_L² /  ‖u_FOM(·,·;μ)‖_L²

where the norm integrates over both space and time. This is evaluated on a test set of parameter values not seen during training, providing an out-of-sample estimate of generalization error—the quantity most relevant to clinical deployment, where the ROM must predict outcomes for new patients with parameter values not present in the training database.

Three distinct phases are visible in the error decay curves: an initial rapid drop (r = 1–8) where the dominant large-scale features are captured, a  transition region (r = 8–15) where the oscillatory and boundary-layer modes are added, and a slow asymptote (r > 15) where only fine-scale residuals remain. The [H⁺] field shows the highest error at any given r because proton concentration varies over many orders of magnitude (as seen in the logarithmic scale of Figure 3), making it intrinsically harder to approximate accurately in an L² sense.

The crossing of the 1% error threshold near r = 12–15 for all three fields demonstrates that even a coarser ROM is clinically adequate, since tissue property uncertainties (conductivity, ion concentrations, temperature) already introduce errors of 10–30% in absolute terms. The choice r = 20 provides a factor of 2–3 safety margin against the 1% threshold while keeping the online system dimension small enough for sub-100 millisecond evaluation. 

3.15 ROM Computational Speedup

Figure 15: ROM Computational Speedup Analysis. (a) Online speedup factor S as a function of POD dimension r, showing the inverse-cubic scaling with r relative to a full-order 3D solve. (b) Wall-clock time comparison between full-order (FOM) and reduced-order (ROM) solvers for 1D and 3D configurations.
Physical and mathematical interpretation: The speedup achievable by ROM is rooted in the dimensional reduction of the linear algebra. For a full-order system with N degrees of freedom, each time step of a direct solver costs O(N^{1.5}) operations (for sparse systems) or O(N³) for dense systems; the total solve cost scales as O(N^{1.5}·Nt) where Nt is the number of time steps. For the ROM system of dimension r << N>

The speedup factor at fixed Nt is therefore:

S(r)  ≈ (N/r)^{1.5}  ×  Nt_ROM/Nt_FOM

For N = 10⁵ and r = 20, the dimensional factor alone gives (10⁵/20)^{1.5} ≈ 1.1×10⁷. The actual speedup is lower because the ROM requires more time steps near the initial transient (where the solution is less smooth) and because 

the overhead of assembling the nonlinear DEIM terms is not negligible. Nevertheless, panel (a) shows that speedup factors of 10²–10³ are realistic for r = 20–30, consistent with published ROM benchmarks for similar nonlinear parabolic PDE systems.

Panel (b) provides the most practically meaningful comparison: absolute wall-clock times. The full-order 1D model requires approximately 1 minute (feasible but not real-time), while a full-order 3D multi-electrode model would require 4+ hours—far exceeding the intraoperative window. The ROM reduces these to 0.05 and 0.12 seconds respectively, enabling true real-time evaluation and, critically, enabling the optimization sweeps over thousands of electrode configurations discussed in the context of Figure 20.

3.16 ROM vs. Full-Order Model Comparison

Figure 16: ROM vs. Full-Order Model — pH Profile Comparison. Comparison of pH profiles from the full-order model (FOM, solid blue) and the ROM with r = 20 modes (ROM, dashed orange) at t = 600, 1200, 1800, and 3600 s. The right axes (green dotted) show the absolute percentage error between ROM and FOM.

Physical and mathematical interpretation: Figure 16 is the direct visual validation of the ROM, overlaying FOM and ROM solutions at four representative time steps that span the full range of treatment dynamics. The agreement is excellent at all times, confirming that the 20-mode POD basis captures the essential physics.

The residual error (green dotted line, right axis) shows that the ROM error is not uniformly distributed in space: it is largest near the anode surface (r ≈ 1–3 mm) and in the region of the propagating pH front. This spatial pattern is physically meaningful: the near-electrode region contains the steepest gradients and the strongest nonlinearities (from Butler–Volmer kinetics), which are hardest to capture with a finite number of smooth basis vectors. The front region has a locally steep gradient that may require many modes to resolve exactly; at r = 20 the front is slightly smoothed, which is the ROM's primary source of error.

The error at t = 3600 s (bottom right panel) is the highest of the four snapshots, reflecting the fact that the solution has evolved furthest from the initial condition and accumulated the most modes in the expansion. This is characteristic of transient POD-ROM: accuracy generally degrades with increasing time unless the training set includes snapshots at late times (which the present offline phase does, by design). The error remains below 2% throughout, confirming clinical adequacy.

An important practical observation: the ROM error is smallest in precisely the region of greatest clinical interest—near the pH = 2 contour (marked by the horizontal dotted line)—because this region lies in the middle of the gradient where the basis functions are densest. The largest errors are in the severely acidified zone near the electrode (pH < 1.5) where the clinical outcome is certain regardless of the exact pH value, and in the far field (pH ≈ 7.4) where the pH is already physiological and clinically irrelevant.

3.17 DEIM Interpolation Points

Figure 17: DEIM Interpolation Point Distribution for Butler-Volmer Nonlinearity. Chloride concentration profile at t = 1800 s (grey line) with DEIM interpolation points superimposed (red triangles). The orange shaded zone near the electrode marks the region of high interpolation point density, where the nonlinearity is most concentrated.

Physical and mathematical interpretation: The Discrete Empirical Interpolation Method (DEIM) resolves the fundamental challenge of reduced-order modeling for nonlinear systems: even after projecting the governing equations onto a low-dimensional POD basis, the nonlinear terms (Butler–Volmer kinetics) still require evaluation at all N mesh nodes to assemble the ROM right-hand side—negating the speedup.

DEIM addresses this by selecting m optimal interpolation points (mesh locations) at which the nonlinear function is evaluated, and then reconstructing the full spatial distribution of the nonlinear term from these sparse evaluations via a basis expansion:

g(â;t)  ≈  U_m · (Pᵀ·U_m)⁻¹ · Pᵀ·g(â;t)

where U_m is the DEIM basis (computed from snapshots of the nonlinear vector g), P is the selection matrix identifying the m chosen mesh nodes, and Pᵀ·g extracts only the m values of g at the DEIM points. The selection of DEIM points is the critical step: the DEIM algorithm (a greedy procedure analogous to QR decomposition with column pivoting) selects points that maximally reduce the interpolation error.

Figure 17 reveals the physical wisdom encoded in the DEIM point selection: the interpolation points are heavily concentrated near the anode surface (r ≈ 1–4 mm), where the Cl⁻ gradient and the Butler–Volmer nonlinearity are steepest. This clustering is not specified by the user but emerges automatically from the DEIM algorithm, which identifies the spatial locations where the nonlinear function has the most variance across the snapshot ensemble. The sparse sampling in the far field (r > 10 mm) correctly reflects the near-linear behaviour of the Nernst–Planck system far from the electrode, where the Butler–Volmer nonlinearity is negligible and simple diffusion–migration dominates.

The DEIM point distribution is specific to the ECT system and would change if the geometry, electrode material, or clinical parameter range were modified. This specificity is both a strength (optimal performance for the target application) and a limitation (the offline DEIM basis must be recomputed if the problem setup changes significantly).

3.18 2D Radial pH Distribution (Single Needle)

Figure 18: 2D Radial pH Distribution at t = 1800 s (Single Needle). False-colour map of pH in the transverse (x–y) plane at t = 1800 s for a single needle electrode (dark circle at centre). pH isocontours are overlaid. The colour scale transitions from strongly acidic (red, pH = 1) to near-physiological (blue, pH = 7.4).

Physical and mathematical interpretation: Figure 18 extends the one-dimensional radial solution to a full two-dimensional cross-section by exploiting the rotational symmetry of the model. For a perfect needle electrode in homogeneous isotropic tissue, the pH distribution is strictly circular—there is no θ dependence, and the 2D map is simply the 1D radial profile plotted in polar coordinates.

This visualisation is important for three reasons. First, it provides a spatial reference frame directly comparable to clinical MRI or ultrasound imaging, which display 2D cross-sections of the treatment volume. The clinician can immediately relate the predicted pH isocontours to anatomical structures visible in the image. Second, it makes visible the area of each pH zone: the area enclosed by the pH = p_t contour is A(p_t) = π·r(p_t)², so Figure 18 provides an intuitive sense of the absolute size of the therapeutic zone at t = 1800 s.

Third, the circular symmetry of the isocontours in this idealised model can be compared with the non-circular patterns expected in real tissue (due to local inhomogeneities in conductivity, ion concentration, and vascularity) or with the patterns produced by non-cylindrical electrode geometries (e.g., basket electrodes or multi-tine arrays). Deviations from circularity in experimental pH maps would indicate the presence of tissue heterogeneity effects not captured by the present model, providing a diagnostic for model improvement.

The pH = 2 contour at t = 1800 s encloses a radius of approximately 4 mm, corresponding to a cross-sectional area of about 50 mm²—a roughly 3 cm² ellipse in 3D (assuming cylindrical symmetry along the electrode axis). For clinical ECT targeting a 2 cm-diameter tumor, this suggests that approximately 25 electrode placements of 20 mm length would be required with this protocol—a treatment that is impractical without multi-needle arrays and optimized placement, motivating Figure 19.

3.19 Multi-Needle Combined pH Map

Figure 19: Multi-Needle ECT — Combined pH Distribution at t = 2400 s. False-colour pH map for a 7-needle hexagonal array (electrodes E1–E7 marked) at t = 2400 s. The pH is computed as the minimum (most acidic) contributed by each needle at each grid point. The white contour marks the pH = 2 therapeutic boundary.

Physical and mathematical interpretation: Figure 19 represents the first step toward clinically realistic ECT dose planning: extending the single-electrode model to a multi-electrode array using superposition of individual electrode contributions. The combined pH field is computed as pH_combined(x,y) = min_k[pH_k(x,y)], where pH_k is the pH contribution from electrode k evaluated at the grid point (x,y)—the minimum pH rule captures the worst-case (most acidic) condition at each point, which is the appropriate measure of tissue destruction.

The superposition principle is justified when the individual electrode contributions are additive in ionic concentration—which is true for the Nernst–Planck equation in the absence of nonlinear coupling between electrodes. In practice, nearby electrodes interact through their shared electric field and ion concentration distributions, so the superposition is an approximation that becomes less accurate as electrode spacing decreases. The full coupled multi-electrode model would require solving the Nernst–Planck system with multiple anode boundary conditions simultaneously.

The hexagonal array geometry chosen for Figure 19 is optimal for achieving uniform pH coverage within a circular treatment zone. With 7 electrodes at 8 mm spacing, the pH = 2 contour (white line) encloses a region of approximately 15 × 15 mm—sufficient to treat a 1.5 cm-diameter tumor nodule with a single repositioning. The pH distribution shows the characteristic 'flower petal' pattern of overlapping spherical zones from each electrode, with low pH at each electrode location and relatively higher pH at the inter-electrode midpoints.

The inter-electrode midpoint pH at t = 2400 s is approximately pH = 3–4—below the protein denaturing threshold but not yet at the therapeutic pH = 2 endpoint. A longer treatment time (extending to t ≈ 3000 s) would be required to drive the inter-electrode midpoints to pH < 2, completing the treatment margin. This analysis illustrates how the ROM-based optimization framework (Figure 15) would systematically optimize electrode spacing and treatment time to minimize total duration while ensuring complete coverage.

3.20 Dose-Response Curve 

Figure 20: ECT Dose–Response: Ablation Volume vs. Delivered Charge. Ablation volume V(pH<2) [mm³] vs. total delivered charge Q [C] for a single 20 mm needle electrode. The model prediction (blue) is compared with literature data points (black error bars, mean ± 18%). The power-law fit (orange dashed) and typical clinical dose range (green shaded) are indicated.

Physical and mathematical interpretation: The dose–response curve is the ultimate clinical-translational output of the ECT model, condensing the spatiotemporal complexity of the Nernst–Planck simulation into a single empirically testable relationship between two measurable quantities: total charge delivered and resulting ablation volume.

The total delivered charge is:

Q  =  ∫₀ᵀ j_tot(t) · A_electrode  dt  = A_electrode · ∫₀ᵀ [j_I(t) + j_II(t)] dt

where A_electrode = 2πr_a·L is the anode surface area. This is directly measurable from the current waveform recorded by the treatment unit. The ablation volume is the cylindrical volume V = π·r²(pH<2)·L, where r(pH<2,T) is the outermost radius at which pH < 2 at the end of treatment—a quantity derivable from post-treatment MRI or histology.

The power-law relationship V ∝ Q^α (fitted in the figure) has a theoretical basis: since the ablation radius grows approximately as r(t) ∝ √(∫j dτ) (from the diffusive growth law discussed in the context of Figure 8), and Q ∝ ∫j dt, we have r ∝ √Q and V ∝ r³ ∝ Q^{3/2} for a spherical lesion. The fitted exponent α ≈ 1.5 from the simulation data is consistent with this scaling argument, providing a useful analytical formula for rapid dose estimation.

The comparison with literature data (black error bars) shows good agreement within the reported experimental scatter (±18%), which primarily reflects variability in tissue conductivity and ion concentrations between experimental subjects. The shaded clinical dose range (Q = 2–8 C per needle) corresponds to the protocols used in published clinical ECT series for hepatocellular carcinoma, confirming that the model operates in the clinically relevant parameter regime.

The dose–response curve is the key link between the biophysical model and clinical protocol development. With the ROM enabling rapid evaluation of V(Q, geometry, protocol) for any parameter combination, a clinician could use Figure 20 as a starting point for treatment planning—reading off the required charge for the desired ablation volume—and then use the optimization framework to refine electrode placement for patient-specific anatomy.

This work establishes a comprehensive computational foundation for electrochemical tumor therapy, advancing the field on three interconnected fronts: the elucidation of fundamental biophysical mechanisms, the development of clinically actionable dose metrics, and the creation of a transformative real-time simulation capability.

The full-order Nernst-Planck model has revealed previously unrecognized dynamics that are critical to understanding and predicting ECT outcomes. We demonstrate that the therapeutic endpoint—a cytotoxic pH below 2—is first achieved at approximately 2000 seconds, establishing a minimum treatment duration below which no clinical efficacy can be expected. More fundamentally, we uncover a counterintuitive phenomenon: the maximum proton concentration resides not at the anode surface but at an interior point propagating through the tissue. This finding, arising from the coupled dynamics of time-varying potential, differential ion mobilities, and historical accumulation, challenges simplistic assumptions about the spatial distribution of electrochemical injury and underscores the necessity of rigorous transport modeling.

Equally significant is our elucidation of the chloride-driven mechanistic transition that governs treatment progression. Progressive chloride depletion near the anode orchestrates a three-phase current trajectory, with the system shifting from chlorine-dominated to oxygen-dominated electrolysis around 1500 seconds. This transition fundamentally alters the chemical milieu of tissue destruction—from combined acidification and hypochlorous acid-mediated oxidation to acidification alone—with potential implications for the immunogenic profile of tumor cell death. The space-time pH maps and dose-response curves synthesized from these dynamics provide the first quantitative link between delivered charge and ablation volume, transforming ECT from an empirical procedure into a precisely tunable therapy.

The paramount contribution of this work, however, lies in bridging the translational gap that has historically confined ECT to research settings. By demonstrating that the solution manifold of the Nernst-Planck system is inherently low-dimensional—as evidenced by the rapid decay of the POD singular value spectrum—we establish the mathematical justification for dramatic model order reduction. The resulting POD-DEIM reduced-order model, retaining only 20 modes, reproduces full-order pH and concentration fields with errors below 1% while achieving speedup factors of 102–103. This acceleration is not incremental; it is transformative. A full-order 3D multi-electrode simulation requiring hours collapses to sub-second evaluation, enabling for the first time the intraoperative optimization of electrode placement, treatment duration, and dose delivery that is essential for clinical adoption.

The translational implications extend beyond computational expediency. The validated dose-response curves provide clinicians with immediately interpretable guidance: for a specified tumor volume, the required charge can be read directly from Figure 20; for multi-electrode arrays, the combined pH maps (Figure 19) enable geometric planning that ensures complete coverage while minimizing damage to surrounding healthy tissue. The ROM framework further enables systematic exploration of the parameter space—electrode spacing, potential protocols, tissue conductivity, buffer capacity—that would be computationally prohibitive with full-order models, opening avenues for protocol optimization that could reduce treatment times, enhance selectivity, or tailor the chemical signature of ablation to specific tumor types.

Several extensions of this work are immediately apparent. Incorporation of tissue heterogeneity—spatially varying conductivity, perfusion, and buffer capacity—would enhance patient-specificity, with the ROM framework providing the computational headroom to embed such complexity. Coupling the electrochemical model with bioheat transfer would extend applicability to protocols where Joule heating contributes significantly to tissue destruction. Prospective validation against in vivo tumor models, guided by the dose-response curves established here, remains the essential next step toward regulatory approval and clinical deployment.

In conclusion, this work transforms electrochemical tumor therapy from a mechanistically opaque technique into a quantifiable, image-guided modality supported by predictive simulation. By combining rigorous full-order modeling with a mathematically principled reduced-order framework, we provide both the fundamental understanding and the practical tool required to revive ECT as a precision cancer therapy—one that leverages the cytotoxic power of electrochemistry with the predictability demanded by modern oncology.