Equilibrium Preservation and Finite Memory Bias in Fractional Monod Chemostat Models

verification.

1. Introduction

The chemostat connects microbial growth with controlled nutrient supply and biomass removal. Monod kinetics provides a classical relation between substrate availability and specific growth [1], while the mathematical theory of the chemostat identifies the operating conditions under which a microbial population persists or washes out [2], [3]. When the model is extended to include temporal memory, the central question is which features of that theory change. Slower relaxation, an altered finite-time detection outcome, and a displaced equilibrium collision are different statements and require different evidence.

A Caputo derivative represents a weighted history of the ordinary derivative [4], [5]. Replacing time derivatives in the biomass and substrate balances therefore changes their transient response. It does not automatically change the algebraic zeros of the reaction and flow terms. A numerical implementation can nevertheless produce apparently displaced boundaries if it omits the initial-state correction or discards part of the convolution without controlling the resulting error. In that case, an observed shift can belong to the discrete recurrence rather than the proposed continuous model.

This issue is particularly important when the same recurrence generates training trajectories and evaluates fitted fractional models. A low prediction error can then reflect faithful reproduction of a numerical artifact. Evidence for a physical memory parameter instead requires a correctly specified initial-value problem, a verified solver, and a calibration experiment that distinguishes order from kinetic and observation parameters. These requirements are relevant even when the fractional model is used only as a phenomenological approximation.

We study a two-state, single-species Monod system with constant dilution, constant yield, and no inhibition or delay. The analysis has three connected parts. We first establish the equilibrium and linear stability structure, including mixed derivative orders, and identify the common-order mass balance. We then derive the exact stationary distortion caused by a raw-state finite-history Grünwald–Letnikov recurrence. Finally, we execute deterministic numerical checks with two full-history methods and provide their computational definitions and implementation. The principal result is an explicit separation between the continuous washout threshold, the stationary branch collision of a discrete map, and a finite-window detection rule.

The equilibrium algebra and the underlying fractional numerical methods are established mathematical tools. The contribution developed here is their joint use to obtain an exact chemostat-specific bias equation, its two collision roots, and a fixed-history refinement limit. The computational results test those deductions and document the accuracy of the implementations used. They are model-based calculations with declared parameters; no laboratory observations or independently recovered historical experiment files are claimed.

2. Literature and Position of the Study

Fractional chemostat research addressed several distinct dynamical settings. Aris and Jamaian [6] investigated a fractional chemostat through stability and bifurcation calculations with an Adams-type predictor–corrector method. Ma et al. [7] considered a delayed Lotka–Volterra food-chain chemostat with incommensurate orders and delay-dependent control. Zhu et al. [8] studied non-monotonic growth under bounded disturbances and established conditions for extinction and persistence. These studies supplied relevant precedents for memory-dependent dynamics, but their conclusions must be read with their own growth laws, state dimensions, delays, and forcing assumptions. A bifurcation mechanism in a delayed food chain cannot be transferred directly to the two-state, delay-free Monod balance considered here.

Recent work also separated memory specifications and physical mechanisms. Elgindy [9] formulated a chemostat with a sliding-memory Caputo operator and periodic boundary conditions. A sliding-window derivative is a model definition with its own history prescription; it differs from approximating a fixed-initial-time Caputo problem and silently deleting its oldest terms. Duan and Pilyugin [10] represented biofilm growth through its effect on the fluid volume and dilution rate. Their physical state structure illustrates why an estimated fractional order cannot by itself identify biofilm retention. Zeinadini and Namjoo [11] used a nonstandard finite-difference approach for a three-dimensional fractional chemostat with variable yields, placing discrete structural preservation within the modeling problem.

The stability tools likewise depend on the derivative structure. Matignon [12] established the familiar spectral-sector condition for a common fractional order. Brandibur et al. [13] reviewed and analyzed single-order and multi-order Caputo systems, while Brandibur and Kaslik [14] characterized stability regions for planar multi-order systems in terms of matrix entries and determinant. Diethelm et al. [15] developed solution and asymptotic results for linear multi-order systems. These works support treating the mixed-order characteristic equation directly. Substituting the largest or smallest order into a common-order criterion has no general justification. Section 3 exploits the particular sign structure of the Monod Jacobian to give a direct root-exclusion argument.

Numerical fractional calculus developed from consistent convolution approximations rather than ordinary time stepping with an added memory term. Lubich analyzed discretized fractional calculus [16]. Diethelm [17], Ford, and Freed [18] introduced an Adams predictor–corrector method and subsequently analyzed its errors. Garrappa [19] synthesized numerical approaches and software practice, and Li and Zeng [20] organized derivative and integral discretizations within a broader numerical framework. The important implementation choices include the derivative convention, the initial correction, the location at which the right-hand side is evaluated, the solution’s regularity, and the treatment of the complete history.

Long-history cost motivated algorithms that approximate history while retaining a controlled connection to the original convolution. Schädle et al. [21] developed fast and memory-saving convolution quadrature. Jiang et al. [22] approximated the Caputo history kernel using sums of exponentials with an explicit accuracy parameter. These approaches differ in purpose from a fixed term cutoff with no tail correction. A fast method is not validated by its storage requirement alone; its kernel approximation and error control are part of its definition.

Regularity near the initial time is another independent source of numerical difficulty. Jin et al. [23] analyzed the L1 scheme for nonsmooth subdiffusion data, and Stynes et al. [24] showed how regularity and mesh grading affected finite-difference convergence. Their results concern fractional diffusion discretizations, rather than a proof for the explicit recurrence used here, but they explain why endpoint error and maximum-over-time error should be distinguished. Garrappa [25] developed reliable numerical evaluation of Mittag– Leffler functions, which provide exact scalar reference solutions. Brugnano et al. [26] subsequently developed a spectral stepwise method using graded meshes and Jacobi expansions. We use two transparent lower-order methods as independent checks, without claiming to benchmark them against these more advanced algorithms.

The inference literature clarifies a further boundary. Raue et al. [27] used profile likelihood to assess practical and structural identifiability in partially observed dynamical models. Bellman and Åström [28] formalized structural identifiability as a question about what input–output observations determine uniquely. Söderström [29] reviewed errors-in-variables identification when measured regressors are noisy, and Akaike [30] developed likelihood-based model identification criteria. These ideas imply that selected fractional orders, an improved residual, and an information criterion are complementary quantities. None can certify that a numerical generator implements the intended continuous dynamics.

Table 1 compares the strands by the mathematical question they addressed. Table 2 gives an interpretive timeline, grouping representative developments rather than asserting a complete historical periodization. The gap addressed here is the missing link between a readily checked equilibrium invariant and the stationary bias induced by a specific history implementation. Sections 3–5 supply that link through formulas and reproducible calculations.

3. Model and Continuous Structure

Let X and S denote biomass and limiting-substrate concentrations in g L⁻¹. The positive parameters are dilution D, feed concentration Sin, maximum specific growth μmax, half-saturation constant K, and biomass yield Y. We fix a physical reference time τ and use dimensionless time s = t/τ, rates d = Dτ and m(S) = τμmaxS/(K + S). For an absolutely continuous scalar function and 0 < α < 1, the fixed-initial-time Caputo derivative is [4], [5].

α 1 s −α

s z(s) = ∫ (s − u) z′(u) du. (1)

Γ(1 − α) 0 At order one, the operator is the ordinary derivative. The model is

s αXX = [m(S) − d]X,

αS m(S)X 0 < αX, αS ≤ 1. (2)

s S = d(Sin − S) − , Initial values X(0) = X0 ≥ 0 and S(0) = S0 ≥ 0 are part of the problem. In physical time, the corresponding

α −1 αi

balances have the form τ i t zi = fi(z), where fi is the conventional concentration-per-time balance. This scaling makes the units explicit. Converting hours to minutes while converting τ and the rates consistently leaves the model unchanged; choosing a different physical reference time with the same numerical coefficients defines a different fractional model. All calculations below fix τ = 1 h.

The right-hand side is locally Lipschitz on the nonnegative quadrant because K > 0. At X = 0 the biomass component vanishes, and at S = 0 the substrate component equals dSin > 0. These boundary signs are consistent with nonnegative continuous solutions under the standard Caputo comparison framework [5]. They do not establish positivity for every numerical step size. The numerical study therefore checks the computed states directly and does not clip negative values.

Proposition 1. The washout state is Ew = (0, Sin). A unique strictly positive equilibrium exists precisely for 0 < d < m(Sin) and is

∗ dK ∗ ∗

S = , X = Y(Sin − S ), dc = m(Sin). (3)

τμmax − d Its linearization is asymptotically stable for every pair (αX, αS) ∈ (0,1]2. In particular, changing the orders does not move the equilibrium collision or create a nonzero imaginary characteristic-root crossing on this positive branch.

Proof: Caputo derivatives of constants vanish. For X∗ > 0, the biomass equation requires m(S∗) = d. Monotonicity of the Monod function gives the unique value in Eq. (3), and the substrate balance then gives X∗. Positivity is equivalent to S∗ < Sin, hence to d < dc. The positive branch meets Ew at dc, independently of the orders.

In the state order (X, S), put a = m′(S∗)X∗ > 0 and c = a/Y = m′(S∗)(Sin − S∗) > 0. The Jacobian is

∗ 0 a ∗

J = ( ) , det(λI − J ) = (λ + d)(λ + c). (4)

−d/Y −d − c For a common order, the eigenvalues −d and −c satisfy the sector condition for every order in (0,1] [12], [13]. For distinct orders, the relevant Laplace-variable characteristic function on the principal branch is:

Q(p) = pαX+αS + (d + c)pαX + dc. (5)

If p = ρeiϕ, with ρ > 0, 0 < ϕ ≤ π/2, and β = αX + αS, its imaginary part is:

ImQ(p) = ρβsin(βϕ) + (d + c)ραXsin(αXϕ) > 0. (6)

The first term is nonnegative because (αX + αS)ϕ ≤ π, and the second is strictly positive. The lower half-plane follows by conjugation. On the nonnegative real axis Q(p) > 0, including Q(0) = dc. Hence no root lies in the closed right half-plane. The standard planar multi-order stability characterization [14], [15] then gives asymptotic stability of the linear system. The argument includes the imaginary axis and excludes a nonzero imaginary crossing.

At washout, the biomass perturbation is triangularly decoupled with coefficient m(Sin) − d, while the substrate self-coefficient is −d. Washout is linearly stable above dc and unstable below it; equality is a nonhyperbolic boundary. Proposition 1 concerns the stated linearizations. It is not a proof of global mixed-order attraction, and it does not cover inhibition, delay, a changing yield, or additional species.

Corollary 1. If αX = αS = α, the substrate-equivalent concentration B = S + X/Y satisfies:

s αB = d(Sin − B), B(s) = Sin + (B0 − Sin)Eα(−dsα). (7)

This follows by adding the two balances after dividing the biomass balance by Y. Here Eα(z) = ∑∞ k=0 zk /Γ(αk + 1) is the Mittag–Leffler function. Its standard scalar relaxation properties imply that B remains between B0 and Sin [5], [25]. In particular, B0 = Sin defines an invariant balance line. When the orders differ, the derivatives cannot be combined into one operator, so this scalar balance and its bound cannot be asserted by cancellation of right-hand sides alone.

Table 3 defines the benchmark. At D = 0.3 h⁻¹, the exact positive equilibrium is (X∗, S∗) = (3.76,0.30) g L⁻¹, c = 1.88, and Dc = 0.4807692308 h⁻¹. Fig. 1 shows the equilibrium branches. There is no positive-branch fold in the admissible Monod range, and the common-order Jacobian has no complex eigenvalue pair. Thus an order–dilution Hopf atlas is not an output of this two-state model’s linearization.

The vertical dotted line is Dc = 0.4807692308 h⁻¹. All admissible derivative-order pairs have these same equilibrium locations; no trajectory was used to construct the curves; a. Biomass branch, and b. Substrate branch (Fig. 1).

4. Initial Correction and Finite History

For component i, let c0 (i) = 1, cj(i) = cj(− i)1(j − 1 − αi)/j, and A( ji) = −cj(i) for j ≥ 1. The raw-state explicit recurrence with a cutoff M is:

min(n+1,M)

zi,n+1 = hαigi(zn) + ∑ A( ji) zi,n+1−j. (8)

Here g is the right-hand side of Eq. (2), and the right-hand side is evaluated at the previous time node. The Grünwald–Letnikov convolution is being applied to zi, whereas the Caputo initial-value formulation requires the convolution of zi − zi,0 [5], [16], [19]. This difference can be tested before a convergence study. Starting Eq. (8) at any nonzero continuous equilibrium gives zi,1 = αizi∗ because A( 1i) = αi. For the common order 0.92, the first state is (3.4592,0.2760); for (0.85,0.95) it is (3.1960,0.2850). These changes occur without a biological disturbance.

For 0 < α < 1, the weights are nonnegative, but the first M do not sum to one. The exact deficit is:

Γ(M + 1 − α)

δM(α) = 1 − ∑ Aj = . (9)

Γ(1 − α)Γ(M + 1)

The identity follows from the finite generalized-binomial sum. It can also be checked inductively using δM+1/δM = (M + 1 − α)/(M + 1) anδ1 = 1 − α d. At α = 1, A1 = 1 and all later weights vanish; we define δM(1) = 0. Thus positive history coefficients form a subunit sum rather than a convex average of the raw states.

Proposition 2. For fixed h > 0 and finite M, every constant stationary state of the filled-history recurrence Eq. (8) satisfies:

δM(αi)

gi(z) = εizi, εi = α (10)

For the Monod balances, put κ = K/(Sin + K), r = dc − εX − κεS, and v = κεXεS. If r > 0 and r2 > 4v, the raw recurrence has a positive stationary branch precisely for d− < d < d+, where

d± = r ± r2 − 4v , dc − d+ = εX + κεS + d− (11)

If both orders are below one, d− > 0; if both equal one, d− = 0 and d+ = dc.

Proof: At a constant state, Eq. (8) becomes zi = hαigi(z) + (1 − δM)zi, giving Eq. (10). The raw stationary washout substrate and candidate positive branch are:

dSin Sw,M = d + εS

S∗ = K(d + εX) (12)

M τμmax − d − εX

∗ d(Sin − SM ∗ ) − εSSM ∗ XM = Y d + εX The numerator in the last expression is positive exactly when SM ∗ < Sw,M, or equivalently m(Sw,M) > d + εX. On the positive-dilution domain this inequality is equivalent to:

F(d) = d2 − (dc − εX − κεS)d + κεXεS < 0 (13)

Under the stated discriminant and sign conditions, the quadratic is negative between its two nonnegative roots, giving Eq. (11). The Monod denominator in Eq. (12) is positive throughout this interval because d + εX < m(Sw,M) < τμmax. The root sum equals r, giving the exact displacement formula. The order-one case follows directly from εX = εS = 0.

The roots are algebraic collisions between the recurrence’s stationary branches. Their calculation is not a stability theorem for the finite-history map, which has additional lagged state variables, and is not a trajectory-based washout detection result. In particular, the small positive lower root in fractional cases is an artificial low-feed boundary created by the stationary losses. It does not imply a corresponding low-dilution boundary in the continuous chemostat.

The displacement identity in Eq. (11) separates biomass loss, the attenuated substrate-loss contribution, and the lower-root correction. It is exact, rather than a fitted regression. Since all terms are nonnegative, the upper collision cannot exceed dc whenever these roots exist. At fixed h, increasing M decreases each positive deficit and moves the upper collision toward dc. Table 4 gives both roots and retains the separately reported finite-window thresholds for comparison.

There is a more consequential refinement limit. For a fixed order below one, the gamma ratio in Eq. (9) gives:

(Mh)−α ε(α, h, M) ∼ Γ(1 − α)

−α (14)

Mh = L fixed: ε → > 0 Γ(1 − α) Thus reducing h while keeping the retained time span fixed does not remove the raw stationary loss. Keeping M fixed is even more problematic for stationary comparisons: The deficit is constant while h−α increases. To remove this particular asymptotic bias, the retained horizon must grow without bound or the omitted tail must be represented consistently. Eq. (14) concerns the stationary loss of Eq. (8), not a general assertion about all short-memory formulations or controlled kernel approximations.

For the mixed pair (0.85,0.95) and L = 7.5, the fixed-history limiting upper collision is 0.451459847 h⁻¹. The continuous value is 0.480769231 h⁻¹, leaving a displacement of 0.029309384 h⁻¹ even as the time grid is refined. Fig. 3 contrasts increasing the history length at fixed step with refining the step at fixed retained duration. The latter can appear numerically settled while still approximating the wrong stationary structure. Proposition 3. The corrected full-history explicit recurrence:

zi,n+1 = zi,0 + hαigi(zn) + ∑n j= + 11 A( ji) (zi,n+1−j − zi,0). (15)

Preserves every continuous equilibrium when initialized at that equilibrium and reduces to explicit Euler when αi = 1. If the sum is instead cut off at a fixed M, any limiting stationary state must satisfy gi(z) = εi(zi − zi,0).

Proof: At an equilibrium initialization every centered history term and every gi vanish, so induction preserves the constant state. For order one, only the first history weight survives and Eq. (15) reduces to zn+1 = zn + hg(zn). For the truncated version at a limiting constant state, substitution gives zi − zi,0 = hαigi + (1 − δM)(zi − zi,0), yielding the stated condition.

Correct initialization is therefore necessary but does not alone validate a finite cutoff. For example, the scalar equation s αu = −ρu with u(0) = u0 > 0 has a true relaxation to zero, whereas the corrected fixed-cutoff stationary value, if reached, is:

u∞,M = u0 (16)

ρ + ε The artificial limit retains the initial value. An equilibrium-only test would miss it because the particular initialization u0 = 0 is preserved exactly. The same issue affects the common-order chemostat mass: corrected truncation gives the stationary balance B∞,M = (dSin + εB0)/(d + ε) rather than Sin, unless the initial balance already equals the feed balance.

5. Numerical Verification

The computations use double-precision Python 3.12.14, NumPy 2.3.5, and SciPy 1.17.0. The two integration paths are corrected full-history GL as defined in Eq. (15) and the one-correction Adams predictor–corrector method [17], [18]. Predictor evaluation, corrector evaluation, and storage of the accepted right-hand side follow the usual PECE sequence. No fitting, stochastic noise, or state clipping is used. Appendix B supplies the executable integration kernels; Table 3 and the equations specify the test configurations.

For each component of the integral form of Eq. (2), the Adams predictor is the product-rectangle quadrature and the corrector is the corresponding product-trapezoidal quadrature. Writing fi,j = gi(sj, zj), the update is:

P hαi (i) zi,n+1 = zi,0 + ∑ bj,n fi,j, Γ(αi + 1)

n (17)

hαi P (i) zi,n+1 = zi,0 + [gi(sn+1, zn+1) + ∑ aj,n fi,j] .

Γ(αi + 2)

The weights for one component are:

bj,n = (n + 1 − j)α − (n − j)α,

a0,n = nα+1 − (n − α)(n + 1)α, (18)

aj,n = (n − j + 2)α+1 + (n − j)α+1 − 2(n − j + 1)α+1.

The last line applies to 1 ≤ j ≤ n; the predictor uses 0 ≤ j ≤ n. Each component uses its own order. At order one, this implementation is the explicit trapezoidal or Heun method, providing an additional integer-order consistency check. The code evaluates these weights directly; it is intended for transparent verification, not large-scale performance claims.

First, the equilibrium test runs for 30 time units at h = 0.05 for all three order pairs. Corrected GL has zero recorded drift in the stored precision and PECE has maximum absolute drift 4.44 × 10−16 g L⁻¹. The raw recurrence produces exactly the first-step states predicted after Eq. (8). Fig. 2.a also shows how the subsequent raw transient can resemble a biological response even though the continuous initial condition is already an equilibrium.

Second, we use two scalar exact references. The half-order relaxation s 1/2u = −0.3u, u(0) = 1, has solution u(s) = erfcx(0.3√s), where erfcx(x) = ex erfc(x). This identity avoids unstable evaluation of a large exponential times a small complementary error function. The second problem has the prescribed smooth solution:

2 0.92 2s1.08 2

u(s) = 1 + s , s u = − 0.3[u − (1 + s )]. (19)

Γ(2.08)

The relaxation example contains the characteristic startup singularity, while the manufactured polynomial isolates a smooth-solution check. Table 5 reports maximum-over-time error for the smooth test and endpoint error for relaxation. For the smooth test, successive GL refinement rates are 1.005, 1.003, and 1.001; PECE rates are 1.931, 1.926, and 1.924. These are observed rates for this test, not a convergence proof for every chemostat trajectory.

The maximum-over-time relaxation errors need separate interpretation. For GL they are 0.003783, 0.004329, 0.003917, and 0.003209 on the same four grids, so they are not monotone on the two coarsest grids. PECE maxima decrease from 0.000757 to 0.000142. A monotone endpoint sequence should therefore not be presented as uniform accuracy over the entire interval. This distinction is consistent with the broader regularity issues emphasized in fractional numerical analysis [23], [24].

Third, we compare complete chemostat trajectories on [0,20] from (X0, S0) = (1,1). The reference is a full-history PECE solution at href = 0.003125, evaluated only at common nodes. The reported discrepancy is:

(h) (ref)

Eh = max |zi,n − zi (sn)|. (20)

Table 6 shows decreasing discrepancies for both methods across every tested order pair. The reference is numerical, not exact. The PECE comparison between h = 0.00625 and h = 0.003125 gives maxima 3.07 × 10−6, 3.12 × 10−6, and 4.42 × 10−6 g L⁻¹ for the three pairs, respectively. These values provide a resolution diagnostic and are not rigorous bounds on the unknown reference error. Every stored concentration in this refinement campaign remained nonnegative; the smallest was above 0.30 g L⁻¹. This empirical check does not establish positivity for arbitrary steps or parameter choices.

Finally, we test the stationary prediction in Eq. (16). With α = 0.92, ρ = 0.3, u0 = 1, h = 0.05, and M = 150, the predicted truncation floor is 0.041702228. At s = 240, corrected truncated GL gives 0.041702228, while corrected full GL and full PECE give 0.001865243 and 0.001865650. Raw truncated GL gives approximately 2.13 × 10−27. Fig. 5 shows the qualitative distinction: One truncation preserves an artificial nonzero floor, while the uncorrected recurrence suppresses the full-history algebraic tail. Agreement with an algebraic stationary prediction validates the diagnosis of the discrete scheme, not its suitability as a Caputo solver.

6. Consequences for Washout Detection and Inference

The equilibrium collision Dc is an algebraic property of the continuous balances. The raw recurrence’s upper collision D+ is an algebraic property of a different stationary equation. A finite-time detection value is a third quantity. For example, a declared endpoint rule could be:

Ddet(T, η, z0) = inf{D ∈ : X(T;D, z0) ≤ ηX0}, 0 < η < 1, (21)

with a specified dilution-search set . If the set is empty the result should be recorded as not detected within the search range. A persistence requirement over a terminal interval would define another rule. Eq. (21) is an example specification, not a reconstruction of every detail of the reported benchmark. The horizon, initial condition, threshold, retained history, and grid resolution can all affect this operational quantity.

Fractional relaxation can change a finite-window outcome even when the equilibrium collision remains fixed. Near washout, the biomass mode has a small growth or decay coefficient, and observation time may be insufficient to distinguish slow persistence from slow extinction. The appropriate experiment reports Dc alongside the chosen Ddet, then repeats the detection analysis at longer horizons and finer steps using full history or an error-controlled approximation. A shifted detection curve can be a legitimate transient result; naming it a shifted equilibrium bifurcation would answer a different question.

The reported thresholds in Table 4 illustrate why this distinction matters. The M = 150 fractional values lie near the calculated raw-recurrence collisions, establishing a concrete possible numerical confound. Proximity does not prove that one mechanism explains every difference. At order one, changing M cannot change Eq. (8), since all weights beyond the first vanish. Consequently, the two reported integer comparisons, 0.4850 and 0.4801, require differences in sweep or detection settings, output processing, or implementation beyond memory length. The original settings and per-run records are needed to resolve that discrepancy.

Parameter estimation should begin only after the forward calculation has passed the numerical checks. For observations yk at times sk and an observation map H, a trajectory-based objective can be written as

J(θ) = ∑ [ yk − Hz(sk; θ)] Σk −1[yk − Hz(sk;θ)], (22)

where θ includes the unknown kinetic parameters, orders, and any explicitly estimated initial conditions, and Σk represents a declared observation covariance. Independent Gaussian errors make this a likelihood-related weighted least-squares objective up to constants; correlated errors require their joint covariance. The equation specifies an inference design, not an estimated model in this paper.

One-step regression on noisy measured states places observation error in the regressors as well as the target. The associated errors-in-variables problem [29] survives even after the solver is corrected. Moreover, a good fit does not distinguish structural uniqueness from practical concentration of an estimate [27], [28]. Orders should be profiled jointly with growth rates, yield, reference time, and initial conditions where those quantities are uncertain. Fixing a physical reference time is especially important when interpreting a non-integer order as a timescale-dependent descriptor.

Information criteria such as AIC arise from a likelihood and a complexity argument [30]. Counting two additional orders after selecting the minimum on a finite order grid does not address all consequences of selection, boundary estimates, or a misspecified observation model. The decisive assessment is a locked multistep rollout on trajectories not used for choosing the orders or solver settings. Independently generated numerical trajectories can test computational recovery; laboratory data and plausible alternative mechanisms are required for a biological claim. No new parameter-recovery or model-selection experiment is asserted here.

7. Discussion and Scope

Swain et al.|Comput. Eng. Technol. Innov. 3(3) (2026) 208-225

The analysis identifies two implementation choices with different consequences. Omitting the initial correction changes the initial-value problem immediately. Cutting off the resulting raw convolution changes the stationary equation after the history fills. Adding the correction repairs the constant-solution test, but a fixed cutoff can still leave initial-condition-dependent stationary bias. These effects should be investigated separately because success on one diagnostic does not imply success on the others.

The exact quadratic provides more information than a visual comparison of washout curves. It exposes the low-dilution root, decomposes the upper displacement, and predicts what happens under different refinement strategies. Increasing the number of stored terms and decreasing the step are not interchangeable operations. At a fixed stored duration, the raw stationary distortion converges to a nonzero value. This supplies a practical explanation for apparently stable numerical thresholds that disagree with continuous equilibrium algebra. The common-order balance offers a complementary structural check. A solver should reproduce the scalar relaxation of S + X/Y and preserve the feed-balance line when initialized on it. Mixed orders remove the ability to combine the balances through one Caputo operator; they should therefore be treated as a distinct modeling decision. They can be mathematically useful without automatically inheriting every conservation or boundedness argument from the common-order model. Whether component-specific memory is physically appropriate depends on the mechanism and measurements.

Table 7 links the literature to the contribution actually established. The goal is not to replace numerical fractional calculus with a new general solver, but to make a model-specific validity test explicit. Fast convolution quadrature and controlled exponential approximations [21], [22] offer routes to larger calculations without deleting history indiscriminately. Higher-order and graded-mesh approaches [24], [26] can improve accuracy when their hypotheses and implementation are respected. Their performance has not been measured in the present two-method comparison.

The study has clear limits. The numerical reference is a refined PECE calculation rather than an exact nonlinear solution, so the chemostat discrepancies are evidence of refinement and cross-method agreement, not certified global error bounds. The exact scalar tests cover selected orders and regularity classes. The step sizes and parameter range do not establish a universal positivity or numerical stability theorem. The stationary collision analysis does not characterize all eigenvalues of the augmented finite-history recurrence.

The model itself excludes inhibition, variable yield, delay, multiple species, and explicit biofilm dynamics. These extensions can alter equilibrium multiplicity and stability and must be analyzed in their own state space [3], [7], [8], [10], [11]. The present linear stability result should not be generalized to them. No measured bioprocess dataset is analyzed, and the previously reported noisy fitting campaign cannot be independently reconstructed without its executable and per-run outputs. The results support numerical and mathematical validity claims at the stated scope, not a validated operational safety margin for a physical reactor.

8. Conclusion

For the two-state constant-yield Monod chemostat, Caputo memory changes transient response while leaving the positive-equilibrium collision at Dc = μ(Sin). The positive-branch linearization is stable for every admissible pair of orders, and the common-order model admits an exact substrate-equivalent balance. These structural facts provide immediate reference checks for simulation.

An uncorrected finite-history GL recurrence fails those checks in a calculable way. Its stationary loss produces a quadratic branch condition, two possible artificial collisions, and an exact upper-shift decomposition. The loss persists when the retained history duration stays fixed under grid refinement. Correcting the initial state restores equilibrium preservation but does not eliminate fixed-cutoff relaxation bias. Executed full-history GL and Adams PECE calculations, exact scalar references, and refinement comparisons support these distinctions. Future washout and parameter-estimation studies can use the resulting tests to separate continuous memory effects from discretization bias before interpreting an operational boundary or a fitted order biologically.

Data and Computational Availability

This study uses declared mathematical models and deterministic calculations, with no new laboratory dataset. Appendix B provides the two integration kernels; Tables 3–6, equations, and captions specify the numerical tests and reference grids. The historical benchmark summaries in Table 4 and Appendix A remain author-reported because the named files exp_25.py, exp_25_results.json, and their per-run outputs were not supplied. The new numerical verification is a separate implementation and is not described as a recovery of that package. No external repository or persistent archive is claimed.

Funding

The research received no specific grant.

Conflicts of interest

The authors declare no conflicts of interest.

Appendix

Appendix A. Reported Fitting Summaries

The available record describes 45 generated series from three order settings, three noise levels, and five seeds per setting, with integer and fractional models fitted to each series. At 2% observation noise it reports the held-out errors in Table 8 and a paired Wilcoxon p-value of 0.0034 for a primary comparison. These values are retained as reported context, without recalculating the statistic or treating it as validation of the Caputo model. Their shared raw-recurrence generator and missing per-series outputs prevent a physical order-recovery conclusion. The deterministic verification in Section 5 does not reproduce these statistical comparisons.

Appendix B. Executable Integration Kernels

The following Python code implements the schemes used in Section 5. Arrays have one column per state, with a separate order for each column. Omitting M retains the complete history; setting corrected=False reproduces Eq. (8). An explicit cutoff is an experimental comparator, not the recommended implementation of the fixed-initial-time Caputo model. The functions use NumPy and SciPy, with no stochastic state and no positivity clipping.

For Table 5, repeat the four stated steps with the two exact test problems. For Table 6, generate a PECE reference at h = 0.003125 for each order pair and compare the coarser arrays with the nested reference nodes, taking the largest absolute componentwise difference. For Table 4, compute the gamma deficits, form r and v in Proposition 2, and use the quadratic roots. For small d−, evaluating d− = v/d+ avoids cancellation. A direct implementation stores O(N) state or right-hand-side values and performs O(N2) history work; no fast-memory complexity is claimed.

import numpy as np from scipy.special import gamma, gammaln, erfcx def gl(fun, z0, alpha, h, T, M=None, corrected=True):

z0 = np.atleast_1d(z0).astype(float) a = np.atleast_1d(alpha).astype(float) N = round(T / h) assert abs(N * h - T) < 1e-10 c = np.ones((N + 1, len(a))) for j in range(1, N + 1):

c[j] = c[j-1] * (j - 1 - a) / j w = -c[1:] z = np.empty((N + 1, len(z0))) z[0] = z0 base = z0 if corrected else np.zeros_like(z0) for n in range(N):

m = n + 1 if M is None else min(n + 1, M) history = z[n+1-m:n+1][::-1] - base z[n+1] = (base + h**a * fun(n*h, z[n]) + np.sum(w[:m] * history, axis=0)) return np.arange(N + 1) * h, z def pece(fun, z0, alpha, h, T):

z0 = np.atleast_1d(z0).astype(float) a = np.atleast_1d(alpha).astype(float) N = round(T / h) assert abs(N * h - T) < 1e-10 q = np.arange(N + 1, dtype=float)[:, None] b = (q + 1)**a - q**a c = (q + 2)**(a + 1) + q**(a + 1) - 2*(q + 1)**(a + 1) z = np.empty((N + 1, len(z0))) f = np.empty_like(z) z[0], f[0] = z0, fun(0., z0) for n in range(N):

pred = z0 + h**a / gamma(a + 1) * np.sum( b[:n+1][::-1] * f[:n+1], axis=0) a0 = n**(a + 1) - (n - a)*(n + 1)**a history = a0 * f[0] if n:

history += np.sum(c[:n][::-1] * f[1:n+1], axis=0) z[n+1] = z0 + h**a / gamma(a + 2) * ( fun((n + 1)*h, pred) + history) f[n+1] = fun((n + 1)*h, z[n+1]) return np.arange(n + 1) * h, z def chemostat(t, z, d=0.3):

x, s = z mu = 0.5 * s / (0.2 + s) return np.array([(mu-d)*x, d*(5-s)-mu*x/0.8]) # Equilibrium preservation at mixed orders.

eq = np.array([3.76, 0.30]) t, z = gl(chemostat, eq, [0.85, 0.95], 0.05, 30.) assert np.max(np.abs(z - eq)) < 1e-10 # Exact half-order relaxation check.

t, z = pece(lambda t, u: -0.3*u, [1.], [0.5], 0.05, 10.) error = np.max(np.abs(z[:, 0] - erfcx(0.3*np.sqrt(t)))) print("Half-order maximum PECE error:", error) # Fixed-cutoff stationary floor, Eq. (16).

alpha, h, M = 0.92, 0.05, 150 delta = np.exp(gammaln(M+1-alpha) - gammaln(1-alpha) - gammaln(M+1)) eps = delta / h**alpha t, u = gl(lambda t, y: -0.3*y, [1.], [alpha], h, 240., M) assert abs(u[-1, 0] - eps/(0.3+eps)) < 1e-9

Figures