\documentclass[11pt,a4paper]{article}

\usepackage[english]{babel}
\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
\usepackage{lmodern}
\usepackage{microtype}
\usepackage{amsmath,amssymb,amsthm,mathtools}
\usepackage{bm}
\usepackage{booktabs,tabularx,array,multirow,longtable}
\usepackage{graphicx}
\usepackage{float}
\usepackage{caption}
\usepackage{subcaption}
\usepackage{enumitem}
\usepackage{xcolor}
\usepackage{geometry}
\usepackage{setspace}
\usepackage{siunitx}
\usepackage{csquotes}
\usepackage[backend=biber,style=numeric,sorting=none,maxbibnames=99]{biblatex}
\addbibresource{references.bib}
\input{data/table_abm_convergence.tex}
\input{data/table_split_metrics.tex}
\input{data/table_state_metrics.tex}
\usepackage[hidelinks]{hyperref}
\usepackage[nameinlink,noabbrev]{cleveref}
\usepackage{listings}
\usepackage{url}

\geometry{margin=2.4cm}
\setstretch{1.12}
\graphicspath{{figures/}}
\sisetup{scientific-notation=true,round-mode=places,round-precision=3}

\newtheorem{theorem}{Theorem}
\newtheorem{proposition}{Proposition}
\newtheorem{remark}{Remark}
\newtheorem{definition}{Definition}
\newcommand{\CD}{\,{}^{\mathrm C}\!D_t^{\alpha}}
\newcommand{\Rzero}{\mathcal R_0}
\newcommand{\vect}[1]{\bm{#1}}
\newcommand{\norm}[1]{\left\lVert #1\right\rVert}

\title{A Two-Hidden-Layer Radial Basis Neural Surrogate for a\
Caputo Fractional-Order Schistosomiasis Model}

\author{%
Gilder Cieza Altamirano$^{1,*}$ and Alvaro H. Salas$^{2,3}$\\[0.4em]
\small $^1$Universidad Nacional Aut\'onoma de Chota, Cajamarca, Peru\\
\small $^2$Department of Mathematics and Statistics, Universidad Nacional de Colombia,\\[-0.15em]
\small Sede Manizales, Manizales, Colombia\\
\small $^3$FIZMAKO Research Group, Colombia\\[0.4em]
\small $^*$Corresponding author: \href{mailto:gciezaa@unach.edu.pe}{gciezaa@unach.edu.pe}\\
\small Alvaro H. Salas: \href{mailto:ahsalass@unal.edu.co}{ahsalass@unal.edu.co};
ORCID: 0000-0001-9343-6062
}
\date{}

\begin{document}
\maketitle

\begin{abstract}
A reproducible numerical framework is developed for a five-compartment
Caputo fractional-order schistosomiasis model coupling susceptible,
infected, and recovered humans with susceptible and infected snail
populations. The compartmental equations are written in a population-balanced
form, and their fundamental properties are established: nonnegativity of
solutions, explicit bounds for the total human and snail populations, the
disease-free equilibrium, and a threshold reproduction number. Reference
trajectories are computed by the fractional Adams--Bashforth--Moulton
predictor--corrector method. A two-hidden-layer radial basis function network,
with 15 and 25 radial units, is then trained as a supervised trajectory
surrogate. The output parameters are refined by the scaled conjugate-gradient
algorithm, while a fixed 70--15--15 split is used for training, validation,
and testing. Three Caputo orders, $\alpha=0.5$, $0.7$, and $0.9$, are studied
on $0\le t\le10$. Across the held-out test sets, the root mean squared errors
are $2.406\times10^{-4}$, $2.122\times10^{-4}$, and
$1.793\times10^{-4}$, respectively, with coefficients of determination above
$0.99996$. The results show that the surrogate reproduces the reference
trajectories accurately and inexpensively after training. The distinction
between the predictor--corrector solver and the neural surrogate is maintained
throughout: the network approximates validated numerical trajectories and is
not presented as an independent replacement for the fractional differential
equations.
\end{abstract}

\noindent\textbf{Keywords:} Caputo derivative; schistosomiasis; fractional
Adams method; radial basis function network; scaled conjugate gradient;
neural surrogate.

\section{Introduction}

Schistosomiasis is an acute and chronic parasitic disease caused by blood
flukes of the genus \emph{Schistosoma}. The transmission cycle couples human
hosts, freshwater snails, and environmental exposure to cercariae. The World
Health Organization estimated that at least 253.7 million people required
preventive treatment in 2024 and reported transmission in 79 countries
\cite{WHO2026}. This burden motivates mathematical models capable of examining
human--snail transmission, treatment, loss of immunity, and long-term control
\cite{Adenowo2015,Chen2010,Okosun2017,Mangal2008}.

Classical compartmental models are usually formulated with ordinary
derivatives. Fractional derivatives provide a nonlocal alternative in which
the current rate of change depends on the history of the state. In epidemic
models, this memory may represent heterogeneous waiting times, accumulated
exposure, delayed behavioral response, or unresolved temporal structure. The
Caputo derivative is especially convenient because it admits standard initial
conditions and reduces continuously to the first derivative when
$\alpha\to1$ \cite{Caputo1967,Podlubny1999,Baleanu2012}.

A five-state fractional schistosomiasis model and a neural approximation were
previously considered by Botmart et al. \cite{Botmart2022}. That study
established an important starting point but left several issues that must be
resolved for a reproducible numerical investigation: the printed equations
contain inconsistent transfer terms, the reference Adams method is not fully
specified, and the neural approximation is sometimes described as though it
were an independent differential-equation solver. In addition, a scientifically
useful study should report initial conditions, integration interval, step size,
data partition, normalization, error definitions, random seed, and executable
code.

The present work reconstructs the model and numerical procedure in a
population-balanced and reproducible form. The specific contributions are:

\begin{enumerate}[label=(\roman*)]
\item a coherent five-compartment Caputo model with explicit human and snail
      population balances;
\item proofs of positivity and boundedness, together with the disease-free
      equilibrium and the threshold number $\Rzero$;
\item a complete fractional Adams--Bashforth--Moulton predictor--corrector
      implementation and a step-refinement check;
\item a two-hidden-layer radial basis function (RBF) surrogate with 15 and 25
      hidden units, trained by scaled conjugate gradients (SCG);
\item transparent training, validation, and test metrics for
      $\alpha\in\{0.5,0.7,0.9\}$; and
\item a supplementary package containing the \LaTeX{} source, figures, data,
      and Python code needed to regenerate the results.
\end{enumerate}

The principal methodological distinction is that the Adams method generates
the reference solution, whereas the neural network learns a compact surrogate
of that solution. This terminology avoids attributing the discretization of
the Caputo operator to a network trained only on trajectory data.

\section{Fractional schistosomiasis model}

\subsection{Caputo derivative and integral representation}

\begin{definition}
For an absolutely continuous scalar function $x$ and $0<\alpha<1$, the Caputo
derivative is
\begin{equation}
{}^{\mathrm C}D_t^{\alpha}x(t)
 =\frac{1}{\Gamma(1-\alpha)}
  \int_0^t (t-s)^{-\alpha}x'(s)\,ds.
\label{eq:caputo}
\end{equation}
For $\alpha=1$, it coincides with $x'(t)$.
\end{definition}

The initial-value problem
\begin{equation}
{}^{\mathrm C}D_t^{\alpha}\vect{x}(t)
 =\vect{f}(t,\vect{x}(t)),
\qquad \vect{x}(0)=\vect{x}_0,
\label{eq:fivp}
\end{equation}
is equivalent to the weakly singular Volterra equation
\begin{equation}
\vect{x}(t)=\vect{x}_0+
\frac{1}{\Gamma(\alpha)}\int_0^t
(t-s)^{\alpha-1}\vect{f}(s,\vect{x}(s))\,ds.
\label{eq:volterra}
\end{equation}
This representation is the basis of the fractional Adams method used below.

\subsection{Compartment definitions and equations}

Let $W(t)$, $A(t)$, and $Q(t)$ denote susceptible, infected, and recovered
humans, respectively. Let $E(t)$ and $S(t)$ denote susceptible and infected
snails. The state vector is
\[
\vect{x}(t)=(W(t),A(t),Q(t),E(t),S(t))^{\mathsf T}.
\]
The transmission structure is shown in \cref{fig:diagram}.

\begin{figure}[H]
\centering
\includegraphics[width=0.98\textwidth]{fig01_compartment_diagram.pdf}
\caption{Five-compartment human--snail transmission structure. The infection
flow from $W$ to $A$ is modulated by infected snails $S$, while the infection
flow from $E$ to $S$ is modulated by infected humans $A$.}
\label{fig:diagram}
\end{figure}

The model is
\begin{subequations}\label{eq:model}
\begin{align}
{}^{\mathrm C}D_t^{\alpha}W
 &=\lambda_1-\alpha_1SW+\gamma Q-\beta_1W,
\label{eq:modelW}\\
{}^{\mathrm C}D_t^{\alpha}A
 &=\alpha_1SW-(\beta_1+\delta)A,
\label{eq:modelA}\\
{}^{\mathrm C}D_t^{\alpha}Q
 &=\delta A-(\gamma+\beta_1)Q,
\label{eq:modelQ}\\
{}^{\mathrm C}D_t^{\alpha}E
 &=\lambda_2-\alpha_2AE-\beta_2E,
\label{eq:modelE}\\
{}^{\mathrm C}D_t^{\alpha}S
 &=\alpha_2AE-\beta_2S.
\label{eq:modelS}
\end{align}
\end{subequations}

The recovery inflow in \cref{eq:modelQ} is $\delta A$, and the snail mortality
rate in both snail equations is $\beta_2$. These choices are required by the
compartmental flow interpretation and make the total populations satisfy exact
balance equations. They correct typographical inconsistencies in the printed
system of Ref.~\cite{Botmart2022}.

\begin{table}[H]
\centering
\caption{Model parameters and values used in the numerical experiments.}
\label{tab:parameters}
\begin{tabularx}{0.96\textwidth}{@{}clXc@{}}
\toprule
Symbol & Value & Interpretation & Units\\
\midrule
$\lambda_1$ & 0.10 & Human recruitment rate & population/time$^{\alpha}$\\
$\alpha_1$  & 0.12 & Snail-to-human transmission coefficient & population$^{-1}$/time$^{\alpha}$\\
$\beta_1$   & 0.15 & Human natural mortality/removal rate & time$^{-\alpha}$\\
$\delta$    & 0.20 & Human recovery rate & time$^{-\alpha}$\\
$\gamma$    & 0.23 & Rate of immunity loss & time$^{-\alpha}$\\
$\lambda_2$ & 0.23 & Snail recruitment rate & population/time$^{\alpha}$\\
$\alpha_2$  & 0.30 & Human-to-snail transmission coefficient & population$^{-1}$/time$^{\alpha}$\\
$\beta_2$   & 0.18 & Snail mortality rate & time$^{-\alpha}$\\
\bottomrule
\end{tabularx}
\end{table}

The initial condition is
\begin{equation}
(W(0),A(0),Q(0),E(0),S(0))=(0.1,0.2,0.3,0.4,0.5).
\label{eq:initial}
\end{equation}
The parameter values are nondimensional benchmark values inherited from the
modeling scale used in Ref.~\cite{Botmart2022}; they are not asserted to be a
country-specific epidemiological calibration.

\section{Qualitative properties}

\subsection{Positivity and boundedness}

\begin{theorem}[Positivity]
Let $0<\alpha\le1$ and let the initial data in \cref{eq:initial} be
nonnegative. Then every solution of \cref{eq:model} that exists on
$[0,T]$ remains in the nonnegative orthant.
\end{theorem}

\begin{proof}
The vector field is quasi-positive. On the coordinate hyperplanes,
\begin{align*}
\left.{}^{\mathrm C}D_t^{\alpha}W\right|_{W=0}
 &=\lambda_1+\gamma Q\ge0,\\
\left.{}^{\mathrm C}D_t^{\alpha}A\right|_{A=0}
 &=\alpha_1SW\ge0,\\
\left.{}^{\mathrm C}D_t^{\alpha}Q\right|_{Q=0}
 &=\delta A\ge0,\\
\left.{}^{\mathrm C}D_t^{\alpha}E\right|_{E=0}
 &=\lambda_2\ge0,\\
\left.{}^{\mathrm C}D_t^{\alpha}S\right|_{S=0}
 &=\alpha_2AE\ge0.
\end{align*}
Hence no component can cross from zero into negative values; the standard
fractional comparison principle gives invariance of $\mathbb R_+^5$.
\end{proof}

Define the total human and snail populations
\[
N_h=W+A+Q,\qquad N_s=E+S.
\]
Adding the corresponding equations gives
\begin{equation}
{}^{\mathrm C}D_t^{\alpha}N_h=\lambda_1-\beta_1N_h,
\qquad
{}^{\mathrm C}D_t^{\alpha}N_s=\lambda_2-\beta_2N_s.
\label{eq:totals}
\end{equation}
Thus
\begin{align}
N_h(t)&=N_h(0)E_{\alpha}(-\beta_1t^{\alpha})
 +\frac{\lambda_1}{\beta_1}
  \left[1-E_{\alpha}(-\beta_1t^{\alpha})\right],
\label{eq:Nh}\\
N_s(t)&=N_s(0)E_{\alpha}(-\beta_2t^{\alpha})
 +\frac{\lambda_2}{\beta_2}
  \left[1-E_{\alpha}(-\beta_2t^{\alpha})\right],
\label{eq:Ns}
\end{align}
where $E_{\alpha}$ is the one-parameter Mittag--Leffler function.

\begin{theorem}[Bounded feasible region]
The set
\begin{equation}
\Omega=\left\{\vect{x}\in\mathbb R_+^5:
N_h\le \max\!\left(N_h(0),\frac{\lambda_1}{\beta_1}\right),\quad
N_s\le \max\!\left(N_s(0),\frac{\lambda_2}{\beta_2}\right)
\right\}
\label{eq:omega}
\end{equation}
is positively invariant.
\end{theorem}

\begin{proof}
For $0<\alpha\le1$ and $z\ge0$, $0<E_{\alpha}(-z)\le1$.
Equations \eqref{eq:Nh}--\eqref{eq:Ns} therefore express each total as a
convex combination of its initial value and its equilibrium level. Combined
with positivity, this proves the result.
\end{proof}

\subsection{Disease-free equilibrium and reproduction number}

The disease-free equilibrium (DFE) is
\begin{equation}
\mathcal E_0=\left(\frac{\lambda_1}{\beta_1},0,0,
\frac{\lambda_2}{\beta_2},0\right).
\label{eq:dfe}
\end{equation}
Following the next-generation construction \cite{VandenDriessche2002}, the
infected variables are $(A,S)$. At the DFE, the new-infection and transition
matrices are
\begin{equation}
F=\begin{pmatrix}
0 & \alpha_1W_0\\
\alpha_2E_0 & 0
\end{pmatrix},
\qquad
V=\begin{pmatrix}
\beta_1+\delta & 0\\
0 & \beta_2
\end{pmatrix},
\label{eq:FV}
\end{equation}
where $W_0=\lambda_1/\beta_1$ and $E_0=\lambda_2/\beta_2$. Therefore
\begin{equation}
\Rzero=\rho(FV^{-1})
=\sqrt{\frac{\alpha_1\alpha_2W_0E_0}
{(\beta_1+\delta)\beta_2}}.
\label{eq:R0}
\end{equation}
For the values in \cref{tab:parameters},
\begin{equation}
\Rzero=0.69769082<1.
\label{eq:R0value}
\end{equation}

\begin{proposition}[Local stability of the DFE]
If $\Rzero<1$, the DFE is locally asymptotically stable for every
$0<\alpha\le1$. If $\Rzero>1$, it is unstable.
\end{proposition}

\begin{proof}
The infected block of the Jacobian is
\[
J_I=\begin{pmatrix}
-(\beta_1+\delta) & \alpha_1W_0\\
\alpha_2E_0 & -\beta_2
\end{pmatrix}.
\]
Its trace is negative, while
\[
\det J_I=(\beta_1+\delta)\beta_2(1-\Rzero^2).
\]
Thus, for $\Rzero<1$, both infected eigenvalues are negative. The remaining
linear modes have negative real eigenvalues associated with population
relaxation and immunity loss. Negative real eigenvalues satisfy Matignon's
criterion $|\arg(\lambda)|>\alpha\pi/2$ for $0<\alpha\le1$
\cite{Matignon1996}. If $\Rzero>1$, the infected block has one positive
eigenvalue, yielding instability.
\end{proof}

\section{Reference solver: fractional Adams predictor--corrector}

Let $t_n=nh$, $n=0,\ldots,N$. For \cref{eq:fivp}, the
Diethelm--Ford--Freed PECE method \cite{Diethelm2002,Diethelm2004} uses the
predictor
\begin{equation}
\vect{x}_{n+1}^{P}=\vect{x}_0+
\frac{h^{\alpha}}{\Gamma(\alpha+1)}
\sum_{j=0}^{n}b_{j,n+1}\vect{f}(t_j,\vect{x}_j),
\label{eq:predictor}
\end{equation}
where
\begin{equation}
b_{j,n+1}=(n+1-j)^{\alpha}-(n-j)^{\alpha},
\label{eq:bweights}
\end{equation}
and the corrector
\begin{equation}
\vect{x}_{n+1}=\vect{x}_0+
\frac{h^{\alpha}}{\Gamma(\alpha+2)}
\left[
\vect{f}(t_{n+1},\vect{x}_{n+1}^{P})+
\sum_{j=0}^{n}a_{j,n+1}\vect{f}(t_j,\vect{x}_j)
\right].
\label{eq:corrector}
\end{equation}
The weights are
\begin{equation}
a_{j,n+1}=\begin{cases}
 n^{\alpha+1}-(n-\alpha)(n+1)^{\alpha}, & j=0,\\
 (n-j+2)^{\alpha+1}-2(n-j+1)^{\alpha+1}
 +(n-j)^{\alpha+1}, & 1\le j\le n.
\end{cases}
\label{eq:aweights}
\end{equation}
For $n=0$, the first weight is interpreted as $a_{0,1}=\alpha$.

The computations use $0\le t\le10$, $h=0.02$, and one corrector evaluation
per step. The direct history summation is $O(N^2)$, which is acceptable for
the 501-point benchmark and has the advantage of transparent reproducibility.

\subsection{Step-refinement verification}

A refinement study for $\alpha=0.9$ compares solutions obtained with
$h\in\{0.08,0.04,0.02\}$ against the solution with $h=0.01$ at common grid
points. The maximum componentwise discrepancies are reported in
\cref{tab:convergence} and displayed in \cref{fig:convergence}. The monotone
reduction confirms that the chosen reference grid is within the expected
convergent regime.

\begin{table}[H]
\centering
\caption{ABM step-refinement results for $\alpha=0.9$.}
\label{tab:convergence}
\begin{tabular}{cc}
\toprule
Step size $h$ & $\ell_\infty$ difference from the $h=0.01$ solution\\
\midrule
\ABMConvergenceRows
\bottomrule
\end{tabular}
\end{table}

\begin{figure}[H]
\centering
\includegraphics[width=0.72\textwidth]{fig_abm_convergence.pdf}
\caption{Step-refinement behavior of the fractional Adams
predictor--corrector method for $\alpha=0.9$.}
\label{fig:convergence}
\end{figure}

\section{Two-hidden-layer RBF neural surrogate}

RBF networks are universal approximators and are particularly effective for
smooth low-dimensional mappings \cite{Broomhead1988,Park1991,Bishop2006}.
Here, time is the scalar input and the five compartments are the outputs. The
normalized input is
\begin{equation}
\xi=\frac{t-t_0}{T-t_0}\in[0,1].
\label{eq:normtime}
\end{equation}

\subsection{First and second radial layers}

The first hidden layer contains $m_1=15$ Gaussian units,
\begin{equation}
h_i^{(1)}(\xi)=
\exp\!\left[-\frac{(\xi-c_i^{(1)})^2}
{2(\sigma_1)^2}\right],
\qquad i=1,\ldots,m_1.
\label{eq:rbf1}
\end{equation}
The centers $c_i^{(1)}$ are equally spaced across the training input range,
and $\sigma_1$ equals 1.35 times the center spacing.

Let $\vect{h}^{(1)}(\xi)\in\mathbb R^{m_1}$. The second layer contains
$m_2=25$ radial units,
\begin{equation}
h_j^{(2)}(\xi)=
\exp\!\left[-\frac{
\norm{\vect{h}^{(1)}(\xi)-\vect{c}_j^{(2)}}_2^2}
{2\sigma_2^2}\right],
\qquad j=1,\ldots,m_2.
\label{eq:rbf2}
\end{equation}
The second-layer centers are sampled uniformly along the ordered training
embeddings, and $\sigma_2$ is 1.8 times the median distance between adjacent
centers. The five-output prediction is
\begin{equation}
\widehat{\vect{x}}(t)=
\vect{\mu}_x+\operatorname{diag}(\vect{s}_x)
\left(W^{\mathsf T}
\begin{bmatrix}\vect{h}^{(2)}(\xi)\\1\end{bmatrix}
\right),
\label{eq:output}
\end{equation}
where $\vect{\mu}_x$ and $\vect{s}_x$ are the mean and standard deviation of
the training targets.

\subsection{Scaled conjugate-gradient training}

The output matrix $W\in\mathbb R^{26\times5}$ minimizes
\begin{equation}
\mathcal L(W)=\frac{1}{2N_{\mathrm{tr}}}
\sum_{n\in\mathcal I_{\mathrm{tr}}}
\norm{W^{\mathsf T}\widetilde{\vect{h}}_n-\vect{z}_n}_2^2
+\frac{\eta}{2}\norm{W_{1:25,:}}_F^2,
\label{eq:loss}
\end{equation}
where $\vect{z}_n$ are standardized targets,
$\widetilde{\vect{h}}_n=[(\vect{h}_n^{(2)})^{\mathsf T},1]^{\mathsf T}$, and
$\eta=10^{-10}$. The SCG algorithm of M\o ller \cite{Moller1993} avoids an
explicit line search by using a finite-difference estimate of the Hessian
along each conjugate direction. A lightly perturbed ridge-regression solution
is used as a stable warm start, after which SCG refines the same regularized
objective. This design keeps the experiment deterministic while retaining the
requested SCG training mechanism.

\begin{remark}
Because the hidden centers and widths are fixed, \cref{eq:loss} is convex in
the output weights. The network is therefore a supervised surrogate of the
ABM trajectories. It does not approximate the Caputo derivative inside the
loss and should not be called a physics-informed or residual-based solver.
\end{remark}

\section{Experimental protocol and error measures}

For each $\alpha\in\{0.5,0.7,0.9\}$, the ABM method produces 501 samples on
$[0,10]$. A pseudorandom permutation generated with seed 20260716 assigns 70\%
of the points to training, 15\% to validation, and 15\% to testing. The same
indices are used for all five outputs. Since points are sampled throughout the
interval, these results measure held-out interpolation rather than strict
future-time extrapolation.

For reference values $x_{n,k}$ and predictions $\widehat{x}_{n,k}$, the
reported measures are
\begin{align}
\mathrm{MAE}&=\frac{1}{5N}\sum_{k=1}^{5}\sum_{n=1}^{N}
|\widehat{x}_{n,k}-x_{n,k}|,
\label{eq:mae}\\
\mathrm{RMSE}&=\left[\frac{1}{5N}\sum_{k=1}^{5}\sum_{n=1}^{N}
(\widehat{x}_{n,k}-x_{n,k})^2\right]^{1/2},
\label{eq:rmse}\\
E_{\infty}&=\max_{n,k}|\widehat{x}_{n,k}-x_{n,k}|,
\label{eq:linf}\\
\mathrm{NRMSE}&=\frac{1}{5}\sum_{k=1}^{5}
\frac{\left[N^{-1}\sum_n(\widehat{x}_{n,k}-x_{n,k})^2\right]^{1/2}}
{\max_n x_{n,k}-\min_n x_{n,k}},
\label{eq:nrmse}\\
R^2&=1-\frac{\sum_{n,k}(\widehat{x}_{n,k}-x_{n,k})^2}
{\sum_{n,k}(x_{n,k}-\overline{x}_k)^2}.
\label{eq:r2}
\end{align}

All calculations were performed in Python using NumPy, SciPy, pandas, and
Matplotlib. The supplied script regenerates every CSV table and every figure.

\section{Results}

\subsection{Reference dynamics and influence of the fractional order}

The reference trajectories for the three fractional orders are compared in
\cref{fig:sensitivity}. All compartments remain nonnegative. Since
$\Rzero<1$, infected humans and infected snails decrease, while the susceptible
populations relax toward the disease-free levels. Smaller fractional order
produces a stronger memory effect and slower long-time relaxation.

\begin{figure}[H]
\centering
\includegraphics[width=0.92\textwidth]{fig_fractional_order_sensitivity.pdf}
\caption{Effect of the Caputo order on the five reference trajectories.}
\label{fig:sensitivity}
\end{figure}

The total populations in \cref{fig:bounds} approach
$\lambda_1/\beta_1=2/3$ and $\lambda_2/\beta_2\approx1.2778$, in agreement
with \cref{eq:Nh,eq:Ns}. This numerical balance is also a useful diagnostic:
an implementation with the incorrect recovery or mortality terms would not
reproduce these identities.

\begin{figure}[H]
\centering
\includegraphics[width=0.78\textwidth]{fig_population_bounds.pdf}
\caption{Human and snail population totals for the three fractional orders,
together with their equilibrium levels.}
\label{fig:bounds}
\end{figure}

\subsection{Reference--surrogate agreement}

\Cref{fig:overlay05,fig:overlay07,fig:overlay09} compare the ABM trajectories
with the RBF--SCG predictions. The curves are visually indistinguishable over
most of the interval. The largest discrepancies occur close to the initial
point, where fractional solutions exhibit the strongest curvature and where a
random validation subset may contain only sparse neighboring training points.

\begin{figure}[H]
\centering
\includegraphics[width=0.90\textwidth]{fig_overlay_alpha_05.pdf}
\caption{ABM reference and RBF--SCG trajectories for $\alpha=0.5$.}
\label{fig:overlay05}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[width=0.90\textwidth]{fig_overlay_alpha_07.pdf}
\caption{ABM reference and RBF--SCG trajectories for $\alpha=0.7$.}
\label{fig:overlay07}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[width=0.90\textwidth]{fig_overlay_alpha_09.pdf}
\caption{ABM reference and RBF--SCG trajectories for $\alpha=0.9$.}
\label{fig:overlay09}
\end{figure}

\subsection{Training, validation, and testing metrics}

\Cref{tab:splitmetrics} reports aggregate metrics for each data subset. The
test RMSE decreases from $2.406\times10^{-4}$ at $\alpha=0.5$ to
$1.793\times10^{-4}$ at $\alpha=0.9$. All test $R^2$ values exceed 0.99996.
Validation errors are larger than test errors because the fixed random split
places several early high-curvature points in the validation subset. This
behavior is reported rather than hidden by selecting a favorable split.

\begin{table}[H]
\centering
\caption{Aggregate neural-surrogate metrics.}
\label{tab:splitmetrics}
\resizebox{\textwidth}{!}{%
\begin{tabular}{c l c c c c c}
\toprule
$\alpha$ & Subset & MAE & RMSE & $E_\infty$ & NRMSE & $R^2$\\
\midrule
\SplitMetricRows
\bottomrule
\end{tabular}}
\end{table}

The state-by-state metrics in \cref{tab:statemetrics} show the same trend. The
largest maximum error, $2.794\times10^{-2}$, occurs for $W$ when $\alpha=0.5$,
while the global $R^2$ for that variable remains 0.999290. For
$\alpha=0.9$, all statewise maximum errors are below $5.2\times10^{-3}$.

\begin{longtable}{c c c c c c c}
\caption{Statewise metrics evaluated on all 501 grid points.}
\label{tab:statemetrics}\\
\toprule
$\alpha$ & State & MAE & RMSE & $E_\infty$ & NRMSE & $R^2$\\
\midrule
\endfirsthead
\toprule
$\alpha$ & State & MAE & RMSE & $E_\infty$ & NRMSE & $R^2$\\
\midrule
\endhead
\StateMetricRows
\bottomrule
\end{longtable}

\subsection{Pointwise errors, regression, and SCG convergence}

The absolute-error histories are displayed in
\cref{fig:error05,fig:error07,fig:error09}. The errors remain small across the
full interval and decay rapidly after the initial high-curvature region. The
regression diagrams in \cref{fig:regression} place the held-out predictions
close to the identity line. The SCG objective and gradient histories in
\cref{fig:scg} decrease monotonically after the warm start, confirming stable
refinement of the output weights.

\begin{figure}[H]
\centering
\begin{subfigure}{0.32\textwidth}
\includegraphics[width=\textwidth]{fig_error_alpha_05.pdf}
\caption{$\alpha=0.5$}\label{fig:error05}
\end{subfigure}\hfill
\begin{subfigure}{0.32\textwidth}
\includegraphics[width=\textwidth]{fig_error_alpha_07.pdf}
\caption{$\alpha=0.7$}\label{fig:error07}
\end{subfigure}\hfill
\begin{subfigure}{0.32\textwidth}
\includegraphics[width=\textwidth]{fig_error_alpha_09.pdf}
\caption{$\alpha=0.9$}\label{fig:error09}
\end{subfigure}
\caption{Pointwise absolute errors of the five neural outputs.}
\label{fig:errors}
\end{figure}

\begin{figure}[H]
\centering
\begin{subfigure}{0.32\textwidth}
\includegraphics[width=\textwidth]{fig_regression_alpha_05.pdf}
\caption{$\alpha=0.5$}
\end{subfigure}\hfill
\begin{subfigure}{0.32\textwidth}
\includegraphics[width=\textwidth]{fig_regression_alpha_07.pdf}
\caption{$\alpha=0.7$}
\end{subfigure}\hfill
\begin{subfigure}{0.32\textwidth}
\includegraphics[width=\textwidth]{fig_regression_alpha_09.pdf}
\caption{$\alpha=0.9$}
\end{subfigure}
\caption{Held-out RBF--SCG predictions versus ABM reference values.}
\label{fig:regression}
\end{figure}

\begin{figure}[H]
\centering
\begin{subfigure}{0.32\textwidth}
\includegraphics[width=\textwidth]{fig_scg_alpha_05.pdf}
\caption{$\alpha=0.5$}
\end{subfigure}\hfill
\begin{subfigure}{0.32\textwidth}
\includegraphics[width=\textwidth]{fig_scg_alpha_07.pdf}
\caption{$\alpha=0.7$}
\end{subfigure}\hfill
\begin{subfigure}{0.32\textwidth}
\includegraphics[width=\textwidth]{fig_scg_alpha_09.pdf}
\caption{$\alpha=0.9$}
\end{subfigure}
\caption{Scaled conjugate-gradient objective and gradient histories.}
\label{fig:scg}
\end{figure}

\section{Discussion}

The reconstructed model has three advantages over the incomplete formulation.
First, every transfer is attached to a source compartment: recovery is
proportional to infected humans, and snail mortality acts consistently on both
snail classes. Second, the exact total-population equations provide analytical
bounds and implementation checks. Third, the threshold number \cref{eq:R0}
connects the observed decay of infection to a standard epidemiological
criterion.

The numerical results support the RBF network as a compact trajectory
surrogate. Once trained, evaluating the network requires only two radial
feature maps and a small matrix multiplication, whereas direct evaluation of
the classical ABM history sum becomes increasingly costly as the time grid is
refined. This makes the surrogate useful for repeated visualization,
parameter-screening prototypes, inverse-problem initialization, or embedding
within a larger optimization loop.

Several limitations must be emphasized. The parameter set is a nondimensional
benchmark rather than a fit to surveillance data. The network is trained
separately for each fractional order and does not generalize automatically to
new parameter vectors. Randomly held-out time points assess interpolation, not
long-horizon extrapolation. Moreover, the surrogate inherits any discretization
error present in the ABM reference data. A future physics-informed version
would include the discretized Caputo residual and the initial-condition
residual in the training objective, but such a method would require a distinct
validation study and should not be conflated with the supervised construction
used here.

A stronger epidemiological study should calibrate the model to a specific
endemic region, quantify parameter uncertainty, compare alternative memory
kernels, and investigate intervention controls. On the computational side,
fast convolution, sum-of-exponentials, or short-memory approximations could
replace the $O(N^2)$ history summation for long simulations.

\section{Conclusions}

A complete and reproducible analysis of a five-compartment fractional
schistosomiasis model has been presented. The corrected model preserves
nonnegativity, admits explicit bounds for the human and snail totals, and has a
disease-free threshold governed by $\Rzero$. For the benchmark parameters,
$\Rzero\approx0.6977$, and the numerical trajectories consistently approach
the disease-free state.

The fractional Adams--Bashforth--Moulton method was used as the reference
solver. A two-hidden-layer RBF network with 15 and 25 radial units was trained
as a supervised surrogate using scaled conjugate-gradient refinement. Across
three fractional orders, the held-out test RMSE remained between
$1.793\times10^{-4}$ and $2.406\times10^{-4}$, with $R^2>0.99996$. These
results demonstrate accurate approximation of the validated trajectories while
preserving the methodological distinction between a numerical solver and a
neural surrogate.

\section*{Author contributions}
G.C.A. contributed to conceptualization, epidemiological modeling, and the
initial manuscript. A.H.S. contributed to mathematical analysis, numerical
methodology, software, validation, visualization, and manuscript revision.
Both authors reviewed and approved the final manuscript.

\section*{Funding}
The authors received no specific funding for this work.

\section*{Conflicts of interest}
The authors declare no conflicts of interest.

\section*{Data and code availability}
All data are generated computationally. The supplementary ZIP package contains
the Python source code, CSV files, figure files in PDF and PNG formats, and the
complete \LaTeX{} source. Running
\texttt{python code/generate\_results.py} from the package directory
regenerates the numerical results.

\section*{Ethics statement}
This work uses no individual-level human or animal data and therefore requires
no ethics approval.

\appendix
\section{Reproducibility algorithm}

The following steps reproduce the reported experiment:
\begin{enumerate}
\item Set the parameters in \cref{tab:parameters}, the initial condition in
      \cref{eq:initial}, and the seed to 20260716.
\item For each $\alpha\in\{0.5,0.7,0.9\}$, solve \cref{eq:model} on
      $[0,10]$ with $h=0.02$ using \cref{eq:predictor,eq:corrector}.
\item Normalize time and standardize each output using training-set statistics.
\item Build 15 first-layer and 25 second-layer Gaussian radial features.
\item Partition the 501 samples into 70\% training, 15\% validation, and 15\%
      testing using the fixed permutation.
\item Refine the output weights by SCG for at most 1200 iterations.
\item Evaluate \cref{eq:mae,eq:rmse,eq:linf,eq:nrmse,eq:r2} and export all
      figures and CSV files.
\end{enumerate}

\section{Software structure}

\begin{description}[style=nextline]
\item[\texttt{main.tex}] Complete manuscript source.
\item[\texttt{references.bib}] Bibliographic database.
\item[\texttt{code/generate\_results.py}] ABM solver, RBF--SCG surrogate,
      tables, and figures.
\item[\texttt{data/}] Reference trajectories, predictions, errors, metrics,
      and metadata.
\item[\texttt{figures/}] Vector PDF and 300-dpi PNG versions of all figures.
\end{description}

\printbibliography

\end{document}
