NeuralPlasmaODE: Physics-Informed Plasma Models

• NeuralPlasmaODE is a family of physics-informed neural ODE models that learn uncertain plasma closures while preserving explicit balance-law structures.
• Its multi-region, multi-timescale formulation has been applied to tokamak burning plasmas, achieving significant error reductions using DIII-D data and transfer learning to ITER scenarios.
• The framework supports forward simulation, control-oriented inversion, and sensitivity analysis, providing a compact surrogate for complex plasma dynamics and fusion simulation.

NeuralPlasmaODE is the name used in recent plasma-physics literature for a class of physics-informed neural ordinary differential equation formulations that retain explicit balance-law structure while learning uncertain closures from data. Its most developed form is a multi-region, multi-timescale transport model for tokamak burning plasmas, trained on DIII-D data and transferred to ITER deuterium–tritium scenarios; related work uses the same name for hybrid neural ODE models of plasma current and internal inductance, and for method-of-lines ideal-MHD solvers with learned local fluxes. Across these uses, the common principle is to keep high-confidence physics in the ODE right-hand side and to learn only the parts that are poorly known, device-specific, or computationally expensive to derive directly (Liu et al., 2024, Liu et al., 2024, Liu et al., 12 Jul 2025, Liu et al., 13 Jul 2025, Wang et al., 2023, Kim et al., 2024).

1. Scope and nomenclature

The term most commonly denotes the burning-plasma transport framework introduced for tokamak dynamics analysis and later specialized to ITER. In that line of work, NeuralPlasmaODE is a data-driven, physics-connected or physics-informed Neural ODE that partitions the plasma into lumped regions and evolves particle and energy balances with learned transport coefficients. A broader reading is also warranted, because the same name has been applied to a hybrid model for the coupled evolution of plasma current and internal inductance on Alcator C-Mod, and related papers describe “NeuralPlasmaODE methods” as a family of physics-embedded neural ODE approaches for plasma dynamics (Liu et al., 2024, Wang et al., 2023, Zhu et al., 2024).

arXiv id	Physical system	Distinctive formulation
(Liu et al., 2024, Liu et al., 2024, Liu et al., 12 Jul 2025, Liu et al., 13 Jul 2025)	DIII-D and ITER burning plasmas	Multi-region, multi-timescale transport Neural ODE
(Wang et al., 2023)	Alcator C-Mod	Hybrid ODE with learned closure for $dV/dt$
(Kim et al., 2024)	Ideal MHD	Method-of-lines ODE with learned local numerical flux
(Zhu et al., 2024)	EGAM mode structure	Physics-embedded ExpNODE in the NeuralPlasmaODE family

This usage pattern suggests that NeuralPlasmaODE is best understood not as a single fixed architecture, but as a family of differentiable plasma-dynamics models organized around the same design choice: explicit governing structure combined with learned closures.

2. Burning-plasma transport formulation

In the ITER-focused formulation, NeuralPlasmaODE is a multi-region, multi-timescale transport model formulated as a Neural Ordinary Differential Equation. The model aggregates dominant burning-plasma processes—alpha self-heating, auxiliary heating, radiation, Coulomb coupling, and internodal transport—into a compact set of nodal ODEs for core and edge, with the scrape-off layer treated through boundary timescales. The evolved state contains region-wise species densities, temperatures, and equivalent energy densities for alpha particles, electrons, and thermal ions (Liu et al., 2024).

In generic form, the state evolves according to

$$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$

with

$$\mathbf{z}(t)=\big\{n_s^{r}(t),\,T_s^{r}(t),\,U_s^{r}(t)=\tfrac{3}{2}n_s^{r}T_s^{r}\big\}_{s,r},$$

where $r\in\{\text{core }(C),\text{ edge }(E)\}$ and $s\in\{\alpha,e,i=D,T\}$. The alpha population is represented by density and an energy-delay associated with slowing-down to electrons and ions. Electrons satisfy charge neutrality, while impurities such as He, Be, and Ar enter through $Z_{\mathrm{eff}}$ and radiation terms (Liu et al., 2024).

The nodal balances include explicit particle sources from external actuation and fusion, transport relaxation between neighboring nodes, and edge ion orbit loss. Energy balances include ohmic heating, auxiliary heating, fusion heating, collisional energy exchange, transport, radiation, and, in the edge, ion orbit loss. Radiation is represented as the sum of electron cyclotron radiation, bremsstrahlung, and impurity line radiation. Alpha slowing-down is treated by delaying ion heating relative to electron heating by a timescale $\tau_{\mathrm{se}}$, so that electrons receive prompt alpha power while ions receive slowed-down alpha power. This arrangement is central to the model’s ability to represent the nonlinear thermal feedback associated with burning plasmas (Liu et al., 2024).

Internodal coupling encodes the model’s multi-timescale character. In the earlier three-node DIII-D version, the plasma is partitioned into core, edge, and scrape-off layer, and inter-node transport is parameterized by characteristic particle and energy times derived from diffusivities and geometry. In the later ITER version, transport times $\tau_P$ and $\tau_E$ connect core to edge and edge to SOL, ensuring flux continuity across interfaces. Fast radiation and Coulomb coupling coexist with slower cross-region transport, making stiffness a fundamental feature of the dynamics (Liu et al., 2024, Liu et al., 2024).

The learned closure enters through the diffusivity law rather than through replacement of the full right-hand side by a black-box network. In log-space, the ITER formulation writes

$$\ln \boldsymbol{\chi}_{r} = \mathbf{b}_{r} + \mathbf{W}_{r}\,\ln \mathbf{x}_{r},$$

and the scalar form generalizes ELMy H-mode H98 scaling through trainable exponents on $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$0, $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$1, $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$2, $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$3, $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$4, $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$5, $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$6, $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$7, and $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$8. In the baseline, electrons and ions share $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$9, and particle diffusivity satisfies $$\mathbf{z}(t)=\big\{n_s^{r}(t),\,T_s^{r}(t),\,U_s^{r}(t)=\tfrac{3}{2}n_s^{r}T_s^{r}\big\}_{s,r},$$0. Positivity and dimensional consistency are built into the parameterization (Liu et al., 2024).

3. Calibration, training, and transfer learning

The original transport-calibration study used 25 DIII-D H-mode, ELMing, non-RMP shots. Profiles from ZIPFIT with EFIT and Thomson scattering were volume-averaged into nodal quantities over core $$\mathbf{z}(t)=\big\{n_s^{r}(t),\,T_s^{r}(t),\,U_s^{r}(t)=\tfrac{3}{2}n_s^{r}T_s^{r}\big\}_{s,r},$$1, edge $$\mathbf{z}(t)=\big\{n_s^{r}(t),\,T_s^{r}(t),\,U_s^{r}(t)=\tfrac{3}{2}n_s^{r}T_s^{r}\big\}_{s,r},$$2, and SOL $$\mathbf{z}(t)=\big\{n_s^{r}(t),\,T_s^{r}(t),\,U_s^{r}(t)=\tfrac{3}{2}n_s^{r}T_s^{r}\big\}_{s,r},$$3. Twenty shots were used for training and five for testing: 131190, 140418, 140420, 140427, and 140535. The optimization objective was a mean-squared error over nodal $$\mathbf{z}(t)=\big\{n_s^{r}(t),\,T_s^{r}(t),\,U_s^{r}(t)=\tfrac{3}{2}n_s^{r}T_s^{r}\big\}_{s,r},$$4, $$\mathbf{z}(t)=\big\{n_s^{r}(t),\,T_s^{r}(t),\,U_s^{r}(t)=\tfrac{3}{2}n_s^{r}T_s^{r}\big\}_{s,r},$$5, and $$\mathbf{z}(t)=\big\{n_s^{r}(t),\,T_s^{r}(t),\,U_s^{r}(t)=\tfrac{3}{2}n_s^{r}T_s^{r}\big\}_{s,r},$$6, with gradients obtained by adjoint sensitivity through the Neural ODE solver (Liu et al., 2024).

Against a baseline multinodal model with empirical $$\mathbf{z}(t)=\big\{n_s^{r}(t),\,T_s^{r}(t),\,U_s^{r}(t)=\tfrac{3}{2}n_s^{r}T_s^{r}\big\}_{s,r},$$7 scaling, the learned parametric diffusivity produced large test-set gains. The reported shot-by-shot MSE reductions were 11.5861 $$\mathbf{z}(t)=\big\{n_s^{r}(t),\,T_s^{r}(t),\,U_s^{r}(t)=\tfrac{3}{2}n_s^{r}T_s^{r}\big\}_{s,r},$$8 0.4075 for shot 131190, 56.6859 $$\mathbf{z}(t)=\big\{n_s^{r}(t),\,T_s^{r}(t),\,U_s^{r}(t)=\tfrac{3}{2}n_s^{r}T_s^{r}\big\}_{s,r},$$9 0.3170 for 140418, 70.3650 $r\in\{\text{core }(C),\text{ edge }(E)\}$0 0.5876 for 140420, 29.7967 $r\in\{\text{core }(C),\text{ edge }(E)\}$1 0.7105 for 140427, and 88.4208 $r\in\{\text{core }(C),\text{ edge }(E)\}$2 0.7348 for 140535, corresponding to $r\in\{\text{core }(C),\text{ edge }(E)\}$3 reduction shot-by-shot and $r\in\{\text{core }(C),\text{ edge }(E)\}$4 on average across the five test shots. The remaining discrepancies were reported to be larger when multiple power sources were applied simultaneously and at the edge where ion orbit loss and gas puffing distributions are influential (Liu et al., 2024).

The ITER study reuses this calibrated transport prior through transfer learning, motivated by the absence of experimental ITER time series. The workflow ingests 2D profiles and global parameters, preprocesses them into uniform time sequences and volume-averaged nodal signals, initializes the diffusivity parameters with DIII-D-trained values, and fine-tunes them to reproduce ITER flat-top design targets. The ITER training set consists of current flat-top conditions from the ITER Technical Basis, while the test set uses plasma start-up transients. Training uses a fixed time step of $r\in\{\text{core }(C),\text{ edge }(E)\}$5, $r\in\{\text{core }(C),\text{ edge }(E)\}$6 s windows, learning rate $r\in\{\text{core }(C),\text{ edge }(E)\}$7, and tens of epochs. Reported examples include an MSE reduction from 6.7085 to 0.0016 in 14 epochs for scenario 2 flat-top training, and from 13.2741 to 0.6333 in two fine-tuning epochs for non-inductive scenario 4 (Liu et al., 2024).

4. ITER burning-plasma dynamics and thermal stability

The transferred model was exercised on three canonical ITER scenarios: inductive scenario 2, hybrid scenario 3, and non-inductive WNS scenario 4. The reported operating points were $r\in\{\text{core }(C),\text{ edge }(E)\}$8, $r\in\{\text{core }(C),\text{ edge }(E)\}$9 MA, $s\in\{\alpha,e,i=D,T\}$0 T, $s\in\{\alpha,e,i=D,T\}$1 s for scenario 2; $s\in\{\alpha,e,i=D,T\}$2, $s\in\{\alpha,e,i=D,T\}$3 MA for scenario 3; and $s\in\{\alpha,e,i=D,T\}$4, $s\in\{\alpha,e,i=D,T\}$5 MA, $s\in\{\alpha,e,i=D,T\}$6 T for scenario 4. Representative global targets included volume-averaged $s\in\{\alpha,e,i=D,T\}$7–$s\in\{\alpha,e,i=D,T\}$8 keV, $s\in\{\alpha,e,i=D,T\}$9–$Z_{\mathrm{eff}}$0 keV, $Z_{\mathrm{eff}}$1–$Z_{\mathrm{eff}}$2, $Z_{\mathrm{eff}}$3–2.07, fusion power $Z_{\mathrm{eff}}$4–400 MW, auxiliary power $Z_{\mathrm{eff}}$5–73 MW, and radiation power $Z_{\mathrm{eff}}$6 of $Z_{\mathrm{eff}}$7–55 MW) (Liu et al., 2024).

In inductive scenario 2, core and edge temperatures reached steady state by approximately 11 s. Early-time heating was dominated by ohmic and auxiliary power until fusion became appreciable. Alpha power first heated electrons, which then transferred energy to ions through Coulomb collisions, supplemented by direct ion fusion heating. Radiation and transport removed excess heat, and no energy excursion was observed. In hybrid scenario 3, stronger auxiliary heating shortened the transient to approximately 6 s, while radiation and transport increased commensurately to offset the higher heating. In non-inductive scenario 4, $Z_{\mathrm{eff}}$8 and $Z_{\mathrm{eff}}$9 saturated by approximately 12 s with higher core temperatures and steeper gradients; because $\tau_{\mathrm{se}}$0 in the core, Coulomb exchange transferred energy from ions to electrons. No runaway was observed in any of the three cases (Liu et al., 2024).

Thermal stability was assessed through the net core power balance and its temperature sensitivity. The sufficient local criterion reported against runaway is

$\tau_{\mathrm{se}}$1

In all three scenarios, the learned transport plus radiative losses were strong enough that this net slope remained negative as $\tau_{\mathrm{se}}$2 and $\tau_{\mathrm{se}}$3 approached steady values. The study therefore concluded that radiation and transport can passively remove alpha-generated heat and suppress thermal runaway in inductive, hybrid, and non-inductive ITER operation (Liu et al., 2024).

A frequent misconception is that alpha self-heating alone necessarily implies unstable burn. Within this model class, that conclusion is not supported: the dominant result is instead the competition between alpha-driven heating and temperature-dependent loss channels, especially radiation and transport.

5. Control-oriented inversion and sensitivity analysis

A later extension recast NeuralPlasmaODE as a control-oriented differentiable simulator. In that formulation, the multinodal plasma model remains physics-informed, but external source profiles—such as NBI power, auxiliary RF heating, and fueling rates—become optimization variables. The state obeys

$\tau_{\mathrm{se}}$4

and the inverse problem minimizes trajectory-tracking error with respect to $\tau_{\mathrm{se}}$5 using automatic differentiation through the ODE solver. The framework supports actuator bounds and ramp-rate constraints in practice, although the reported paper focuses on unconstrained tracking. No new training of transport exponents was performed there; instead, previously calibrated transport closures were reused (Liu et al., 12 Jul 2025).

This inversion viewpoint converts the forward transport model into a source-design tool. The stated workflow is: define nodes and state vector, choose target trajectories, specify the physics-informed Neural ODE with calibrated transport, set a tracking objective with optional source magnitude and smoothness penalties, integrate forward, differentiate the loss with respect to $\tau_{\mathrm{se}}$6, and update source waveforms with a gradient-based optimizer. The paper presents the method and implementation rather than detailed quantitative validation, and notes that runtime, robustness metrics, and full D-T validation remain future work (Liu et al., 12 Jul 2025).

Sensitivity analysis provides a complementary use of the same differentiable structure. For a scalar output $\tau_{\mathrm{se}}$7 and parameter $\tau_{\mathrm{se}}$8, the normalized sensitivity is

$\tau_{\mathrm{se}}$9

Around a nominal ITER inductive scenario 2 model, the reported sensitivities identified transport exponents tied to magnetic field and safety factor as dominant for ion temperatures. Examples include $\tau_P$0, $\tau_P$1, $\tau_P$2, and $\tau_P$3. Electron-temperature sensitivities to transport exponents were smaller, but impurity sensitivities were more consequential: $\tau_P$4 and $\tau_P$5, with argon more influential than beryllium. ECR effects were modest, for example $\tau_P$6, and edge IOL sensitivities were effectively negligible, such as $\tau_P$7 (Liu et al., 13 Jul 2025).

Finite perturbation studies reinforced the same ranking. Halving core–edge thermal diffusivity produced $\tau_P$8 keV and $\tau_P$9 keV, while doubling it produced $\tau_E$0 keV and $\tau_E$1 keV. The authors interpret this as self-regulation: increasing temperature also increases transport through the explicit $\tau_E$2 and $\tau_E$3 dependence of $\tau_E$4, damping further heating. This is a central physical implication of the learned closure in the ITER setting (Liu et al., 13 Jul 2025).

6. Other NeuralPlasmaODE formulations, limitations, and outlook

Outside burning-plasma transport, the name has been used for a hybrid neural ODE formulation for predicting coupled plasma current $\tau_E$5 and internal inductance $\tau_E$6 on Alcator C-Mod. That model retains two exact Romero equations, $\tau_E$7 and $\tau_E$8, while learning only the low-confidence closure for $\tau_E$9 through a multilayer perceptron. It was trained on 489 plasma-producing pulses with irregular sampling, mean $$\ln \boldsymbol{\chi}_{r} = \mathbf{b}_{r} + \mathbf{W}_{r}\,\ln \mathbf{x}_{r},$$0 ms and standard deviation $$\ln \boldsymbol{\chi}_{r} = \mathbf{b}_{r} + \mathbf{W}_{r}\,\ln \mathbf{x}_{r},$$1 ms. On test data, the hybrid model achieved $$\ln \boldsymbol{\chi}_{r} = \mathbf{b}_{r} + \mathbf{W}_{r}\,\ln \mathbf{x}_{r},$$2 mean percent error for $$\ln \boldsymbol{\chi}_{r} = \mathbf{b}_{r} + \mathbf{W}_{r}\,\ln \mathbf{x}_{r},$$3 and $$\ln \boldsymbol{\chi}_{r} = \mathbf{b}_{r} + \mathbf{W}_{r}\,\ln \mathbf{x}_{r},$$4 for $$\ln \boldsymbol{\chi}_{r} = \mathbf{b}_{r} + \mathbf{W}_{r}\,\ln \mathbf{x}_{r},$$5, outperforming both the physics-only baseline at $$\ln \boldsymbol{\chi}_{r} = \mathbf{b}_{r} + \mathbf{W}_{r}\,\ln \mathbf{x}_{r},$$6 and $$\ln \boldsymbol{\chi}_{r} = \mathbf{b}_{r} + \mathbf{W}_{r}\,\ln \mathbf{x}_{r},$$7, and a pure neural ODE at $$\ln \boldsymbol{\chi}_{r} = \mathbf{b}_{r} + \mathbf{W}_{r}\,\ln \mathbf{x}_{r},$$8 and $$\ln \boldsymbol{\chi}_{r} = \mathbf{b}_{r} + \mathbf{W}_{r}\,\ln \mathbf{x}_{r},$$9. The reported interpretation is that targeted hybridization improves sample efficiency, generalization, and robustness against the overfitting seen in the pure learned ODE (Wang et al., 2023).

A distinct but conceptually related usage appears in ideal MHD, where NeuralPlasmaODE denotes a method-of-lines ODE system whose right-hand side is computed from a learned local numerical flux implemented by a Flux Fourier Neural Operator. There the state is the conservative ideal-MHD vector $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$00, and the learned operator replaces hand-designed Riemann-solver components while preserving conservation-form time stepping. In 2D tests, example relative $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$01 errors at $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$02 were approximately $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$03 for $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$04, $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$05 for $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$06, $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$07 for $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$08, and $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$09 for $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$10. The reported inference time per $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$11 was $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$12 s for Flux NO versus $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$13 s for WENO-Z, or approximately $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$14 speedup, though local instabilities appeared at higher resolution and long times (Kim et al., 2024).

The EGAM study based on ExpNODE places itself in this broader physics-embedded NeuralPlasmaODE family. It predicts 2D mode profiles using exposed latent states, symmetry-preserving decoding, and a known linear coupling between $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$15 and density. On the reported data, test loss for 20-step prediction was 0.0219 for the physics-embedded ExpNODE, compared with 0.0411 for a physics-free ExpNODE and 0.0688 for ConvLSTM; for 40-step prediction the corresponding losses were 0.08772, 0.2208, and 0.2053. The explicit lesson drawn there is that embedding even a minimal set of physical constraints can substantially improve long-horizon stability and out-of-distribution generalization (Zhu et al., 2024).

Across the literature, limitations are consistent. The burning-plasma transport models use only a small number of explicit spatial regions; SOL and divertor dynamics are reduced to timescales or sink terms; triton diffusivities are set equal to deuteron diffusivities; $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$16 in the baseline; RF and NBI deposition profiles are predetermined; atomic and molecular neutral processes, recycling, pedestal dynamics, ELMs, MHD activity, and detailed impurity transport are not resolved. The inverse-control work remains only partially validated quantitatively. The C-Mod inductance model learns $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$17 without direct supervision on $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$18 in the loss. The learned-flux MHD formulation does not explicitly enforce positivity of $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$19 or $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$20, and its stability degrades in harder regimes or at higher resolutions (Liu et al., 2024, Liu et al., 12 Jul 2025, Wang et al., 2023, Kim et al., 2024).

Future directions stated in the cited works include explicit pedestal and ELM dynamics, coupling to SOL and divertor physics, refined impurity and synchrotron or ECR radiation models, integration with frameworks such as OMFIT for self-consistent $$\frac{d\mathbf{z}}{dt} = f\big(\mathbf{z}, t; \boldsymbol{\theta}\big),$$21-profile and current-drive evolution, uncertainty quantification over transport parameters, validation on experimental discharges, actuator-constrained control, closed-loop model predictive control, richer invariants for mode-structure prediction, and extensions to 3D or resistive MHD. Taken together, these directions indicate a continuing effort to make NeuralPlasmaODE a compact but physically structured surrogate for fusion-plasma simulation, inference, and control rather than a replacement for full first-principles integrated modeling.