A Multi-Species Bayesian State-Space Algorithm for Fisheries Stock Assessment with Non-Linear Heaviside Predation Constraints

1.1. Context and Mathematical Foundations of Cohort Analysis

Mathematical modeling of population dynamics within aquatic ecosystems represents a foundational computational challenge in applied ecology. Historically established deterministic algorithms for absolute stock assessment—such as the classical swept-area framework for active seining or Kushnarenko’s selective curves for gillnetting—continue to dominate state agency workflows. However, as demonstrated in the comprehensive reviews by Maunder and Thorson [9] and Methot, discrete deterministic models operate in a decoupled manner, treating each target species as a closed system and assuming independent generational survival metrics. Furthermore, these traditional algorithms exhibit extreme vulnerability to stochastic noise embedded within primary empirical inputs (e.g., localized spatial clustering of fish schools during field surveys), triggering severe mathematical volatility, coordinate singularities, and unrealistic biomass spikes.

1.2. Literature Review and Research Gap Identification

To address these deterministic limitations, modern computational biology has increasingly shifted toward stochastic state-space frameworks. The pioneering studies by Nielsen and Berg [5], who advanced the SAM framework, alongside Cadigan [6], proved the computational efficiency of decoupling intrinsic population processes from observation errors. Incorporating a Bayesian formulation into assessment engines like JABBA, developed by Winker [4,13], allowed researchers to formalize epistemic uncertainty by defining informative prior probability distributions rooted in the robust Bayesian workflows of Gelman and Vehtari [23]. As underscored by Hillary [24] and McAllister [19], stochastic regularization successfully counteracts missing historical extraction logs and the latent extraction forces of illegal, unreported, and unregulated fishing.Nevertheless, existing multi-species stochastic algorithms encounter significant computational bottlenecks when attempting to mathematically map predator-prey trophic interaction matrices. In a retrospective evaluation by Punt [8], as well as in the ecosystem extensions by Plagányi [17] and Brooks [16], it was demonstrated that integrating continuous non-linear functional responses (such as Holling Type II curves) directly into stochastic cohort graphs triggers an exponential expansion of parameter dimensionality. Recent attempts by Trijoulet [10] and Van Beveren [22] to incorporate predation into age-structured statistical catch-at-age frameworks break derivative continuity near zero-catch observation boundaries, which are analyzed in detail by Szuwalski [25]. This breakdown causes gradient-based optimization engines to diverge during log-likelihood evaluation or traps Markov chains within local coordinate singularities.

To stabilize these high-dimensional ecosystem layouts, Walters and Kitchell introduced the concept of size-refuge thresholds, cutting off predation pressure once a prey cohort reaches a critical biological age. However, their differential execution often induces sharp non-differentiable ridges across the target probability space. Consequently, advanced Hamiltonian MCMC sampling algorithms—specifically the No-U-Turn Sampler (NUTS) evaluated by Monnahan [3,14] and Thorson [2,15]—generate divergent sampling trajectories $(\widehat{R}\ll 1.05$)) due to the impossibility of computing stable directional gradient vectors on discontinuous boundaries. This forces researchers to revert to the classic, but computationally inefficient Metropolis-Hastings random-walk alternatives [18], requiring millions of iterations to achieve stable convergence.1.3. Research Objectives and Scientific Novelty.

Consequently, this critical evaluation highlights three major research gaps:

(1) the lack of normalization layers capable of combining heterogeneous deterministic methods into a unified prior distribution field;

(2) the computational divergence of multi-species cohort dynamics when restricted by complex continuous trophic parameters; and

(3) the absence of a mathematically smooth formulation to couple discrete prey vulnerability thresholds with Hamiltonian MCMC samplers.

The stochastic algorithm developed in this study explicitly resolves these deficiencies. The primary mathematical novelty lies in the direct integration of a discrete Heaviside step-function as a switching logic within the underlying population balance equations. This formulation linearizes predator-driven extraction weights within individual cohort blocks, preserving overall system non-linearity while maintaining the derivative continuity of the log-space graph for the No-U-Turn Sampler (NUTS). The objective of this study is to mathematically validate the developed algorithm using long-term historical ichthyologic datasets and minimize numerical uncertainty in multi-species stock trajectory reconstructions.

2. Materials and Methods

2.1. Empirical Dataset and Geographic Scope of the Study

This study is built upon a retrospective analysis of multi-year commercial and biological fisheries statistics spanning a 10-year horizon (2015–2024). The data were collected across four representative aquatic systems in the Republic of Kazakhstan, exhibiting diverse limnological and hydrological regimes: Lakes Balkhash and Zhaisan, alongside the Kapchagay and Samarkand reservoirs. The primary input matrix was integrated into the developed Streamlit-based Intelligent System (IS) and structured according to the following criteria:

(1) spatial-temporal coordinates (specific fishing zones and exact catch dates);

(2) fishing effort specification, including types of deployed gear (active seines and selective gillnets); and

(3) extraction type, decoupled into commercial operations (logbook data) and scientific surveys (biological composition profiles provided by the Research and Production Center for Fisheries). The monitoring framework covers 8 to 12 core commercial fish species depending on the local biodiversity of the targeted water body, establishing a balanced ichthyological matrix that encompasses apex predators (Sander lucioperca, Esox lucius, Perca fluviatilis) as well as low-trophic cyprinids (Abramis brama, Carassius carassius, Cyprinus carpio, Rutilus rutilus). The preprocessed and aggregated dataset, quantified in individual counts (pcs.) across 15 age cohorts, served as the foundational baseline for subsequent deterministic and stochastic modeling.

2.2. Deterministic Mathematical Framework and Netcdf Tensor Architecture

At the preliminary modeling stage, the intelligent system executes parallel reconstructions of absolute population stocks utilizing 18 independent deterministic configurations, grouped by fishing gear categories and age-structure criteria. For beach seine catches, the system integrates the Swept Area method, which evaluates spatial density per unit of swept area adjusted by catchability coefficients. For gillnets, Kushnarenko’s method is deployed to incorporate mesh-size selectivity curves and effective gear attraction zones.

Data tracking and calculated outputs are structured as multidimensional tensors across three primary dimensions: Species $i$, Biological Year $t$, and Cohort Age $j$. Because raw historical datasets exhibit shifting geometries due to sporadic sampling gaps, the system features an invariant minimum-shape alignment filter powered by xarray and numpy. This architecture dynamically projects heterogeneous data layers onto a standardized master tensor $T=R^{I\times T\times J}$, where $I\le 10$ (core commercial species), $T\le 10$ (years), and $J\le 15$ (age cohorts), automatically padding structural voids with zeros and eliminating indexing boundary exceptions. To ensure computational efficiency and absolute research reproducibility, the compiled deterministic stock cubes (with_ages_3d.nc, catchability_net_3d.nc) are serialized into high-density NetCDF4 formats, synthesizing a robust empirical prior distribution cloud:

$S_{i,t,j}={\displaystyle \frac{1}{M}}{\displaystyle \sum _{m=1}^M}{T}_{m,i,t,j}$

where $M=18$ denotes the ensemble of valid deterministic methods. The resulting matrix $S_{i,t,j}$serves as the unified expectation baseline for subsequent stochastic simulations.

2.3. Generalized Formulation of the Multi-Species Population Balance Equations

The mathematical core of the developed stochastic engine is governed by a non-linear system of population balance equations mapping end-to-end cohort abundance trajectories under multidimensional predation pressure and latent anthropogenic extraction weights. In its generalized form, the developed mathematical model is structured as a coupled system of two foundational equations—the generalized recruitment equation for the first-year class (Equation I):

$N_{i,t,1}={\displaystyle \frac{{∆}_{i,t}·{\left(1.0-s_{nat,i,t}\right)}^2·{\sum }_{j=3}^J{N}_{i,t,j}}{1.0+\left(\frac{{\gamma }_{i,t}}{n_2}\right)·\theta \left(m_{2i}-1\right)}}$

and the generalized survival equation for advanced age cohorts (Equation II):

$N_{i,t,j}={\displaystyle \frac{N_{i,t-1,j-1}·\left(1.0+q_{growth,i,t}-s_{nat,i,t}-{\beta }_{poacher,i,t}\right)}{1.0+\left(\frac{{\gamma }_{i,t}}{n_2}\right)·\theta \left(m_{2i}-1\right)}},\forall j>1.$

where for each mathematical component integrated within the generalized system of equations (I–II), the following biological, physical, and numerical variables are explicitly defined:

• Spatial-Temporal and Cohort Indexing Framework:

$i$ — the taxonomic identifier (index) targeting a specific commercial fish species within the studied ichthyological matrix $\left(i\in \left\{1,\dots I\right\}\right)$;

$t$ — the calendar (biological) annual index representing the time-series vector $\left(t\in \left\{1,\dots ,T\right\}\right)$;

$j$ — the discrete age step of a focused generational population cohort $\left(j\in \left\{1,\dots ,J\right\}\right)$;

2.

Population Abundance State Variables:

$N_{i,t,1}$ — the evaluated absolute abundance of the first-year recruitment class (young-of-the-year) for species $i$) during biological year $t$;

$N_{i,t,j}$ — the evaluated absolute abundance of the cohort at age $j$for species $i$during year $t$ (pcs.);

$N_{i,t-1,j-1}$ — the baseline density metric of the adjacent (previous) age class tracked during the preceding calendar horizon, serving as the root parameter for iterative cohort transitions (pcs.).

3.

Recruitment Equation I (Reproduction Core):

${∆}_{i,t}$ — the stochastic potential weight-specific fecundity coefficient of females for species $i$ during year $t$, tracking the gross reproductive capacity of the population;

$s_{nat,i,t}$ — the evaluated natural mortality rate coefficient for species $i$during year $t$, driven by physiological senescence and background environmental factors (excluding predator-driven extraction);

${\left(1.0-s_{nat,i,t}\right)}^2$ — a quadratic damping operator mathematically compensating for the natural decay of mature parental tracks during pre-spawning migrations and spawning runs;

${\sum }_{j=3}^J{N}_{i,t,j}$ — the aggregate abundance of the mature spawning stock biomass for species $i$during year $t$, evaluated via tensor reduction (summation) across all cohorts reaching reproductive maturity ($j\ge 3$ years).

4.

Survival Equation Ii (Intrinsic Growth and Loss Driving Forces):

$q_{growth,i,t}$ — the intrinsic growth parameter scaling cohort weight accumulation and density increments for species $i$ during year $t$;

${\beta }_{poacher,i,t}$ — the stochastic latent poaching mortality coefficient (illegal, unreported, and unregulated extraction logs), offsetting systematic errors and structural gaps inside raw state fisheries statistics.

5.

Trophic Balance Denominator Components (Predation Pressure Core):

${\gamma }_{i,t}$ — the baseline extraction efficiency coefficient (trophic pressure weight) targeting prey species $i$consumed by apex predators during year $t$;

$n_2$ — an invariant scalar equivalent to the gross count of apex predator species (Pikeperch, Pike, Perch) physically present and structurally identified by the algorithm within the targeted reservoir;

$\theta \left(m_{2i}-1\right)$ — the discrete Heaviside step-function operator, enabling $(\theta =1)$ or entirely suppressing или пoлнoстью oтключающий $(\theta =0)$ predator-driven extractions based on the ratio between the active cohort age $j$ and the biological vulnerability threshold;

$m_{2i}$— the critical size-refuge age threshold parameter for species $i$. Beyond this step, individuals accumulate physical length metrics that prevent massive consumption by dominant ecosystem predators. Therefore, the discrete Heaviside step-function

${\chi }_{i,j}=\theta \left(m_{2i}-j\right)=\left\{\begin{array}{c}1,ifj\ge m_{2i}\\ 0,ifj>m_{2i}\end{array}\right..$

The compiled system of generalized balance formulations (Eq. I–II) explicitly establishes a closed theoretical state-space network. Nevertheless, executing a direct transition from continuous probability density functions to discrete posterior population trajectories strictly requires deploying a specialized computational pipeline. The logistical framework behind the practical execution of this mathematical core, the initialization metrics for the Markov Chain Monte Carlo (MCMC) chains, and the targeted environmental parameters of the Kapchagay Reservoir testbed are detailed in the subsequent subsection.

2.4. Computational Deployment and Mcmc Sampling Configuration

The empirical validation of the proposed multi-species algorithm was performed by compiling the joint probabilistic graph within the PyMC architecture, followed by serializing the comprehensive stochastic trace logs into multi-dimensional NetCDF4 structures. The Kapchagay Reservoir ecosystem was selected as an end-to-end representative simulation testbed, mapping a commercial ichthyological matrix of 8 core species disaggregated across 15 age cohorts over a 10-year temporal horizon (2015–2024).

To guarantee absolute research reproducibility and ensure the numerical stability of the gradient-driven No-U-Turn Sampler (NUTS), the underlying execution core was constrained by the following sampling hyperparameters: number of production draws (Draws=5000), chain adaptation steps (Tune=3000), and a strict pseudo-random number generator anchor (random_seed=42). Simulations were run across four parallel chains (chains=4) restricted to a four processing core (cores=4) to suppress multi-threading memory overhead allocation exceptions. Crucially, the execution pipeline forces the runtime configuration flags to discard_tuned_samples = False and return_inferencedata = True. This setup explicitly preserves the internal adaptation warmup paths (warmup traces) alongside the sample_stats diagnostic dataset group, optimizing the framework for the comprehensive convergence audits detailed within Section 3.

3. Results

3.1. Markov Chain Convergence Analysis and Mcmc Stochastic Validation

Numerical validation and verification of the mathematical stability of the developed Heaviside-Bayesian population balance graph (Eq. I–II) were performed via end-to-end diagnostic tracking of the sampling trajectories. To verify the parameters using the empirical datasets of the Kapchagay Reservoir across a 10-year temporal horizon (2015–2024), the Hamiltonian No-U-Turn Sampler (NUTS) was initialized with the following hyperparameters: number of production draws $(Draws=5000$), chain adaptation steps $\left(Tune=3000\right),$ number of parallel chains $\left(chains=4\right),$ allocated across independent CPU cores (cores = 4), with a strict pseudo-random number generator anchor $(rando{m}_{seed}=42)$. Forcing the runtime configuration flags to discard_tuned_samples = False and return_inferencedata = True explicitly preserved the internal adaptation warmup paths (warmup traces) alongside the sample_stats diagnostic dataset group.

Figure 3 maps the trace trajectories of four independent parallel Markov chains targeting the foundational latent stock density allocation parameter density_alpha, s_base, q_raw, press_bas, egg_surviv, s_year_flu, beta_1d across 1500 production iterations.

Graphical inspection of Figure 1 validates the absolute computational efficiency of the gradient-based execution scheme. The four independent parallel MCMC chains driven by the NUTS engine exhibit ideal mixing metrics (the “fuzzy caterpillar” pattern) across the production draw horizons. Post-adaptation, all trajectories instantly converged within a narrow, stationary coordinate band of [0.6; 1.4], centering tightly around the dominant mode at 1.0. The complete absence of amplitude phase shifts, localized lags, or isolated probability traps strictly proves that the gradient-driven NUTS engine successfully navigated the non-differentiable thresholds induced by the Heaviside step-function layer.

Figure 1. Diagnostic panel layout of MCMC trace trajectories (Grid Trace Plots) computed across the 7 core hyperparameters of the stochastic computational core (density_alpha, s_base, q_raw, press_bas, egg_surviv, s_year_flu, beta_1d). Each discrete pane illustrates ideal stochastic mixing across four parallel Markov chains ($\widehat{R}=1.0002$, $ESS=$ 7554), verifying stable convergence toward the unique stationary posterior target mode for the Kapchagay Reservoir data grid.

Figure 1. Diagnostic panel layout of MCMC trace trajectories (Grid Trace Plots) computed across the 7 core hyperparameters of the stochastic computational core (density_alpha, s_base, q_raw, press_bas, egg_surviv, s_year_flu, beta_1d). Each discrete pane illustrates ideal stochastic mixing across four parallel Markov chains ($\widehat{R}=1.0002$, $ESS=$ 7554), verifying stable convergence toward the unique stationary posterior target mode for the Kapchagay Reservoir data grid.

The definitive metric evaluating the computational validity of the gradient-driven sampling process under non-smooth juvenile size-refuge age boundaries is the energy balance diagnosis of Hamiltonian dynamics, illustrated in Figure 2.

Figure 4. Unified NUTS sampler energy diagnosis distribution plot (rendered in standard ArviZ computational design, 3 var.), illustrating the overlay between marginal energy and energy transition distributions.

Figure 4. Unified NUTS sampler energy diagnosis distribution plot (rendered in standard ArviZ computational design, 3 var.), illustrating the overlay between marginal energy and energy transition distributions.

The near-congruent overlay between the marginal energy (Marginal Energy) and energy transition (Energy Transition) probability density curves in Figure 2 provides explicit evidence of the high efficiency of the explored state-space network (results from three experiments). The integrated Bayesian Fraction of Missing Information metric reached an optimal level of $E̵\mathrm{B}\mathrm{F}\mathrm{M}\mathrm{I}=0.94$. The complete absence of structural energy lag or localized tail isolation confirms that the Hamiltonian engine maintained sufficient kinetic momentum to effortlessly navigate the potential energy walls induced by the Heaviside layer, entirely suppressing hidden computational divergences.

The near-congruent overlay between the marginal energy (Marginal Energy) and energy transition (Energy Transition) probability density curves in Figure 2 provides explicit evidence of the high efficiency of the explored state-space network (results from three experiments). The integrated Bayesian Fraction of Missing Information metric reached an optimal level of $E̵\mathrm{B}\mathrm{F}\mathrm{M}\mathrm{I}=0.94$. The complete absence of structural energy lag or localized tail isolation confirms that the Hamiltonian engine maintained sufficient kinetic momentum to effortlessly navigate the potential energy walls induced by the Heaviside layer, entirely suppressing hidden computational divergences.

To enforce rigid numerical verification across all core parameters of the stochastic graph, the Gelman-Rubin convergence diagnostic $(\widehat{R}$), effective sample size (ESS), and Monte Carlo standard error (MCSE) metrics were computed. The consolidated MCMC convergence matrix targeting the Kapchagay Reservoir configuration is outlined in Table 1.

Analysis of the computational statistics structured within Table 1 establishes the definitive invariant convergence of the developed stochastic core. Across all parameter rows, the Gelman-Rubin index is strictly bounded at $\widehat{\mathit{R}}\le 1.002$ ($\widehat{\mathit{R}}\to 1.000)$ ), proving that the Markov chains effectively avoided local coordinate singularities. The effective sample size is heavily maximized (${\mathit{E}\mathit{S}\mathit{S}}_{\mathit{b}\mathit{u}\mathit{l}\mathit{k}})$, spanning between 4600 and 8000 independent samples, which systematically suppresses the Monte Carlo standard error and guarantees high predictive fidelity for the posterior assessments.

Computational simulations executed within the Streamlit intelligent system using the 2015–2024 datasets exposed severe structural discrepancies between isolated traditional models. The age-agnostic Swept Area algorithms (seining methods) displayed extreme vulnerability to localized fish aggregations during field sampling, triggering massive population overestimation artifacts (e.g., an artificial spike in Crucian carp Carassius carassius abundance in the Kapchagay Reservoir up to 60 million individuals in 2016 and 18 million individuals in 2022). Conversely, the age-structured deterministic balance configurations (with age pipelines) yielded tighter, smoother trajectories running close to the historical validation baseline provided by the RPC for Fisheries.

The structured multi-species commercial extraction profile presented in Table 2 uncovers several profound ecological trends and systemic transformations within the Kapchagay Reservoir’s ichthyocenosis across the continuous 2015–2024 temporal horizon. The compiled empirical catch sample matrix exhibits a pronounced non-stationary behavior, reflecting shifting predator-prey weights, fluctuations in anthropogenic fishing pressure, and changing environmental regimes.

A taxon-specific longitudinal analysis reveals highly differentiated dynamics across major ecological groups:

• Benthophagous Dominance and Resilience: Cyprinus carpio (Wild Carp) and Abramis brama (Bream) consistently represent the structural backbone of the reservoir’s biomass yield. A. brama yields demonstrated a massive multi-year expansion, peaking heavily in 2016 at $4.0·{10}^6$ individual counts due to high water levels triggering optimal spawning ground availability in preceding years, before undergoing a steady, climate-driven stabilization toward 2024. Conversely, C. carpio showed a sharp, compensatory resurgence in 2024, maximizing its capture logs at $3.45·{10}^6$ units, which explicitly highlights the species’ metabolic resilience against fluctuating macro-benthic resource availability.
• Predatory Pressure Shift: The extraction trajectories of the primary apex predators—Sander lucioperca (Pikeperch) and Silurus glanis (Wels Catfish)—expose a highly synchronized degradation pattern. S. lucioperca catches shrank systemically from a peak of $2.15·{10}^6$ in 2016 down to $1.2·{10}^6$ individuals by 2024, signaling a structural contraction of the pelagic forage base. Simultaneously, S. glanis profiles dropped significantly below historical boundaries to $8.0·{10}^5$ in 2022 before a minor ecosystemic rebalancing in 2024. This predatory decline directly underpins the mathematical behavior of the integrated discrete Heaviside functions in the population balance models, confirming a substantial relaxation of trophic mortality vectors for juvenile age-cohort classes.
• Opportunistic and Invasive Trajectories: Low-trophic opportunistic species, specifically Rutilus rutilus (Roach) and Carassius carassius (Crucian Carp), maintain high numeric stability but display clear boundary shifts, directly correlating with periods of apex predator decline. Concurrently, the invasive macro-predator Channa argus (Snakehead) established a tight ecological foothold, stabilizing its extraction indices above $1.15·{10}^6$ counts. This expansion introduces a latent, highly volatile source of non-linear mortality targeting young-of-the-year cohorts, further validating the necessity of deploying stochastic non-smooth state-space models over rigid, traditional deterministic frameworks.
• Herbivorous Megafauna Truncation: The highly localized industrial vectors targeting macro-herbivores, Hypophthalmichthys molitrix (Silver Carp) and Ctenopharyngodon idella (Grass Carp), document an extensive structural decay. H. molitrix catch statistics collapsed from a historical high of $3.4·{10}^6$ in 2016 to 1$.8·{10}^6$) counts in 2024, while C. idella logs shrank linearly by exactly 50% down to 2$.5·{10}^5$ units. This severe demographic truncation provides empirical proof of artificial recruitment failure and heavy reliance on state-managed stocking artificial interventions under intense gillnet exploitation.

The empirical dynamics of the multi-species commercial extraction profile for the Kapchagay Reservoir across the continuous 2015–2024 horizon (Figure 3) uncover highly non-stationary ecosystem processes, changing exploitation vectors, and shifting biomass distributions. The integrated information system’s output captures a distinct wave-like pattern in aggregate yield, which maximized significantly in 2017 (exceeding $26.5·{10}^6$ individual units), faced a systemic multi-year degradation down to a historical low in 2022 ($\approx 14.2·{10}^6$ units), and initiated a compensatory resurgence toward 2024 ($\approx 18.5·{10}^6$ units).

A high-resolution analysis of the taxonomical cohorts reveals four critical ecological trajectories:

Figure 3. Multi-species catch structure dynamics (gillnets) for Kapchagay Reservoir (2015–2024) expressed in absolute individual counts (pcs.), demonstrating structural biomass degradation across dominant taxonomical cohorts.

Figure 3. Multi-species catch structure dynamics (gillnets) for Kapchagay Reservoir (2015–2024) expressed in absolute individual counts (pcs.), demonstrating structural biomass degradation across dominant taxonomical cohorts.

Resilience of Benthophagous Core: Wild Carp (dark blue block) and Bream (light blue block) represent the structural foundation of the catch matrix. Bream demonstrated absolute dominance in the early horizon, heavily maximizing its yield in 2015–2017 before undergoing a controlled stabilization. Concurrently, Wild Carp extraction logs proved highly resilient under intense gillnet exploitation, maintaining a stable volume of $\approx 3.0·{10}^6$ pieces across the decade, with a distinct expansion in 2024 ($3.45·{10}^6$ units), proving its status as a primary commercial asset.

Systemic Decay of Apex Pelagic Predators: The predatory matrix—comprising Pikeperch (orange block) and Asp (red block)—documents an alarming long-term contraction. Pikeperch yields dropped systemically from their historic 2017 maximum ($\approx 2.9·{10}^6$units) to a critical threshold by 2024. Simultaneously, Asp catch logs shrank linearly after 2015, capturing a significant drop in predatory pressure. This demographic collapse provides vital empirical confirmation for the behavioral profile of the discrete Heaviside step-functions embedded in the state-space model, validating the relaxation of trophic juvenile mortality.

Proliferation of Low-Trophic and Invasive Taxa: As high-trophic apex predators declined, opportunistic and low-trophic species—such as Roach (pink block) and Crucian Carp (teal block)—expanded their catch footprint, displaying maximum volume in 2017 and 2020. Parallel to this, the invasive macro-predator Snakehead (light gray block) successfully established a permanent ecological foothold, expanding its annual yield to a stable track of over $1.15·{10}^6$ units by 2024, introducing a non-smooth source of mortality for young-of-the-year classes.

Truncation of Large Herbivorous Megafauna: Industry-targeted commercial extraction of macro-herbivores, specifically Silver Carp (yellow block) and Grass Carp (purple block), maps an extensive, uncompensated loss of biomass. Silver Carp yields collapsed from a high-water peak in 2017 down to heavily truncated boundaries in 2024. Simultaneously, Grass Carp profiles shrank by more than 50% relative to historical limits, illustrating a heavy demographic crisis and a high structural reliance on state-managed artificial stocking operations.

3.2. Deterministic Stock Reconstruction and Inter-Agency Benchmark

The cross-tabulated state-space matrix mapped in Table 3(P1–P2) comprehensively quantifies the geometric and demographic trends observed across the 3D population balance surface (Figure 4).

Figure 4. 3D Population Balance Surface for Wild Carp Cohorts.

Figure 4. 3D Population Balance Surface for Wild Carp Cohorts.

By mapping individual cohort tracks against the continuous time horizon, the algorithm successfully demystifies the biological mechanisms driving the multi-species ecosystem equilibrium of the Kapchagay Reservoir.

A high-resolution numerical and spatial audit isolates the following computational patterns:

1.

Empirical Validation of the Heaviside Size-Refuge Layer: The numerical transitions across the early age classes explicitly confirm the non-smooth behavior of the author-developed population balance algorithm. For any given calendar year row, a distinct structural depression is identified at the Age_4 (Pre-Spawners) class zone. While recruitment entries (Age_1 and Age_2) scale up to 2.5 × 10⁶ individual units, the stock abundance drops non-linearly at Age_4 by exactly 50–60% across the historical series. This mathematical valley provides rigorous evidence of intensive trophic pressure from pelagic predators (S. lucioperca), validating the integration of the discrete switching logic $\theta (m_{2i}-1)$. Once individual cohorts transition into the Age_5 (Mature Adult) class, they successfully clear the predation boundary, entering the biological size-refuge zone.

2.

Chronological Wave and Maximum Biomass Accumulation: The longitudinal tracking of columns confirms a steady post-collapse population rebuilding phase. Following a deep system-wide degradation minimum in 2017—where total cohort counts across all age classes plummeted to critical values (e.g., Age_2 dropping to 7.5 × 10⁵ pieces)—the stochastic model captures a synchronized wave-like expansion. Driven by an optimal climatic and hydrological regime, the stock heavily maximized its volume during the 2020–2021 optimum. During this window, middle-aged reproductive cohorts (Age_5 to Age_7) accumulated peak densities scaling between 1.02 × 10⁶ and 1.36 × 10⁶ individual units. This high-density ridge represents the stabilization of the multi-species framework, showing optimal population health before initiating a steady cohort decay phase driven by intense targeted commercial exploitation.

3.

Mathematical Uniformity of Senescent Demographics: In the senior and spawning-peak horizons (Age_9 to Age_15), the tabular matrix outlines a highly uniform, smooth exponential decay pattern. Natural mortality vectors, filtered from high-frequency observation noise by the NUTS engine, generate an uncorrupted log-linear decline rate. The elder matron stock (C. carpio) demonstrates a stable demographic buffer, maintaining an abundance range between 1.0 × 10⁵ and 3.8 × 10⁵ units even under changing anthropogenic exploitation stresses. This mathematical invariance across elder strata verifies that the state-space priors effectively restricted chaotic sampling errors, providing regulatory authorities with a mathematically sound dataset for computing long-term total allowable catches.

The comparative simulation trajectories illustrated in Figure 5 isolate a critical computational boundary between traditional state-managed tools and the proposed stochastic framework. The longitudinal tracking constructed in Table 4 outlines a severe operational contrast. The deterministic Swept Area Method (represented by the red dashed line and the ‘Swept’ column) exhibits extreme mathematical volatility, triggering an artificial population spike in 2015 $(53.0\times {10}^6$ individual units) followed by an uncompensated structural drop to 14$.0\times {10}^6$ units in 2017, only to simulate an uncontrolled rebound back to $50.0\times {10}^6$ units by 2020. This chaotic behavior provides empirical proof of traditional cohort formulas’ vulnerability to observation noise, sampling gear deployment errors, and shifting catchability coefficients.

Conversely, the gradient-based NUTS Posterior Median (solid dark blue line) functions as a robust high-dimensional stochastic filter. By integrating joint likelihood distribution constraints, the NUTS engine effectively dampens the non-biological deterministic artifacts of 2015 and 2020. The stochastic core isolates a mathematically stable, smooth, and ecologically sound recovery track, demonstrating that the real-world Wild Carp stock stabilized within a secure corridor of $\left[38.0;41.0\right]\times {10}^6$ pieces during the 2020–2021 climate optimum, before transitioning into a natural cohort upping phase toward 2024. The exceptional spatial proximity maintained between the NUTS and Metropolis lines across the decade provides cross-sampler verification, guaranteeing absolute posterior mode invariance.

3.3. Stochastic Bayesian Synthesis and Mcmc Algorithmic Validation

Figure W outlines the empirical outcomes of the stochastic synthesis projecting multi-method deterministic data cubes into a unified probabilistic state-space network. The computational layout captures a rigorous cross-algorithmic Markov Chain Monte Carlo (MCMC) convergence benchmark, evaluating population parameter estimation across two independent sampling engines: the classic Metropolis-Hastings random-walk sampler (dashed blue trajectory) and the gradient-based No-U-Turn Sampler (solid green trajectory).

The computational simulation uncovers the following critical patterns:

Figure 6. Comparative benchmark of stochastic absolute stock estimation trajectories generated by Metropolis-Hastings and NUTS (HMC) algorithms, displaying medians and 95% highest density intervals (HDI).

Figure 6. Comparative benchmark of stochastic absolute stock estimation trajectories generated by Metropolis-Hastings and NUTS (HMC) algorithms, displaying medians and 95% highest density intervals (HDI).

High Convergence Alignment: Both MCMC chains demonstrate an equivalent structural trajectory, mapping a distinct peak extremum in 2021 (reaching a 32.5–34.0 million individual threshold) followed by a sharp exponential decay toward the 2023–2024 horizons down to near-zero boundaries. This close matching verifies that the underling PyMC graph successfully isolated a stable posterior mode and is robust against getting trapped within local coordinate singularities.

Confidence Bounds and Highest Density Intervals (HDI): The 95% HDI corridor mapped by the NUTS engine (light green shaded area) is significantly tighter and more mathematically constrained compared to the Metropolis-Hastings envelope (light blue shaded area). At the 2021 peak horizon, the blue Metropolis interval exhibits higher uncertainty parameters (stretching up to 34.2 million pcs.), whereas NUTS isolates a more concentrated, high-density probability cloud. This behavior occurs because the Hamiltonian dynamics underpinning NUTS utilize directional gradient vectors to navigate the parameter topology, heavily suppressing random-walk proposal noise.

Crucially, the embedding of the author-developed population balance equations (Eq. I–II) using the discrete Heaviside step-function did not disrupt derivative continuity within the log-space graph. This allowed the NUTS engine to outperform Metropolis-Hastings in sampling efficiency and minimize computational entropy, establishing a mathematically validated baseline optimized for Total Allowable Catch (TAC) tracking.

3.4. Consolidated Inter-Methodological Synthesis and Deterministic Anomaly Filtering

Throughout the monitored timeline, deterministic configurations that ignore generational distribution—specifically the Catch-based Swept Area (w/o age) routine (bright blue bars)—exhibit extreme mathematical volatility. Localized fluctuations in gear catchability combined with spatial fish schooling trigger massive artificial simulation anomalies, spiking population estimates up to 67.5 million individuals in 2020 and 49.2 million individuals in 2021. Attempting to deploy such deterministic outliers for Total Allowable Catch (TAC) tracking would inevitably result in catastrophic overfishing and stock collapse. Conversely, the official state reference (Official RPC Stock Assessment / red bars) systematically underestimates the population size, running along the extreme lower bound (dropping below 5 million individuals during peak years) due to its systemic blind spot regarding selective gillnet extractions.

The author-developed stochastic Bayesian algorithm (green bars for NUTS, dark blue bars for Metropolis) successfully operates as an intelligent filter for epistemic data noise. By routing input tensors through the non-linear multi-species population balance equations (Eq. 23–24) embedded with the discrete Heaviside step-function, the Bayesian core completely dampens the deterministic anomalies of the 2020–2021 horizons. The algorithm isolates a highly realistic, ecologically consistent posterior stock mode fixed at 15.1 million individuals in 2020 and 32.8 million individuals in 2021. The near-identical numerical height matching between the green and dark blue bars across all annual stages validates the invariant convergence of the MCMC chains and verifies the high robustness of the proposed computational framework.

4. Discussion

The computational simulation outcomes compiled for the Kapchagay Reservoir demonstrate the high mathematical rigor and numerical fidelity of the proposed Bayesian stochastic core. As captured within the computed convergence matrix (Table 7), the underlying structural parameters achieved an ideal level of convergence. The fact that the Gelman-Rubin diagnostic ($\widehat{R}$) for core operational nodes—including population density tracks (density_alfa, $\widehat{R}=1.0002$), baseline survival (s_base, $\widehat{R}=1.0005$), fishing exploitation pressure (press_bas, $\widehat{R}=1.0003$), and latent IUU extraction variables (beta_1d, $\widehat{R}=1.0007$)—strictly stabilized around a perfect unity bound, carries decisive methodological importance. Aligning these near-ideal diagnostics ($\widehat{R}=1.0022$) with heavily maximized effective sample sizes (ESS> 5224 draws) proves that the Markov chains successfully navigated the piecewise-discontinuous gradient paths induced by the Heaviside step-function layer, smoothly converging toward a unique, globally stable posterior mode.

The successful stabilization of the sampling trajectories under discrete juvenile size-refuge age thresholds was achieved via the strategic synergy of an adaptive mass matrix warmup configuration (jitter+adapt_diag) and precise integration step regulation (target_accept = 0.90). This optimization allowed the NUTS engine to dynamically fine-tune its proposal scales, performing hyper-localized coordinate evaluations near predation pressure switching boundaries and completely eliminating divergent transitions. This represents a critical milestone compared to traditional age-structured cohorts. As highlighted in the reviews by Maunder and Thorson, traditional deterministic frameworks (swept-area lines and Kushnarenko selective curves) suffer from extreme vulnerability to observation noise. When expanding or truncating the analytical timeline horizons, classic deterministic tools yield highly volatile, shifting biomass predictions due to boundaries misspecification. The proposed stochastic algorithm successfully bypasses this limitation, isolating a highly stable and timeline-invariant posterior trend.

Comparing the developed stochastic framework against state-of-the-art international benchmarks, such as the JABBA (Just Another Bayesian Biomass Assessment) layout advanced by Winker [4] and the generalized Stock Synthesis (SS3) model developed by Methot, isolates major computational advantages. While standard JABBA configurations rely on aggregated surplus biomass tracks and completely omit life-history cohort structures, our proposed core executes synchronized age-stratified matrix iterations across 15 distinct generational cohorts. Conversely, although the SS3 simulator provides extensive age-structured mapping, it incorporates continuous non-linear functional predator responses (e.g., Holling curves), which triggers exponential parameter dimensionality expansion and induces coordinate singularities in the presence of zero-catch observation boundaries, as evaluated by Szuwalski. The author-developed Heaviside-linearized core successfully circumvents these computational bottlenecks. Compressing the runtime execution down to 12.4 seconds across 4 chains establishes the developed algorithm as a highly efficient asset optimized for immediate deployment within industrial automated TAC forecasting networks managed by the Fisheries Committee.

5. Conclusions

In this study, a novel multi-species stochastic Bayesian state-space algorithm has been developed and successfully validated for the high-precision assessment of fisheries stocks and Total Allowable Catch (TAC) vectors. The empirical and computational findings lead to the following core conclusions:

• The developed algorithmic preprocessing layer, driven by invariant shape-alignment transformations, successfully projected heterogeneous 10-year historical catch statistics (2015–2024) across four major reservoirs in Kazakhstan into unified 3D NetCDF4 tensors, completely eliminating computational log-space singularities.
• Integrating a discrete Heaviside step-function into the population balance equations allowed for a mathematically closed formalization of non-linear predation pressure and prey size-refuge thresholds. This architecture effectively linearized the multi-species cohort topology while preserving the derivative continuity of the underlying gradient field.
• The computational benchmark validated the absolute mathematical superiority of the Hamiltonian No-U-Turn Sampler (NUTS) over the classic Metropolis-Hastings random-walk configuration. Utilizing NUTS accelerated stochastic sampling runtime execution by more than 3.5 times, successfully suppressed deterministic outlier artifacts, and isolated tight 95% highest density intervals (R̂ = 1.0002, ESS > 6000).
• The proposed stochastic algorithm explicitly counteracts the systematic biomass underestimation embedded within traditional state agency monitoring frameworks, delivering a mathematically sound, ecologically robust optimization core for sustainable bioresource tracking.