Beyond diagonal noise: A better predator-prey modeling framework with cross-covariance

Figures Abstract The introduction of stochasticity into continuous ecological models frequently relies on phenomenological, diagonal diffusion terms that lack a rigorous microscopic basis. We demonstrate that this standard practice fundamentally misrepresents the geometry of demographic fluctuations. By deriving a stochastic Rosenzweig–MacArthur model directly from an integer-valued, Bernoulli-coupled continuous-time Markov chain, we isolate the exact diffusion covariance structure dictated by event stoichiometry. We mathematically prove that coupled predation–conversion events inherently generate a structurally negative predator–prey cross-covariance, exposing the severe mathematical and biological limitations of standard diagonal-noise approximations. Furthermore, we resolve a persistent ambiguity in stochastic population modeling by explicitly formalizing the bifurcation between open-domain formulations (for survival-conditioned interior dynamics) and absorbed formulations (for extinction-permitting dynamics). To rigorously support this distinction, we develop a tailored two-stage Lyapunov well-posedness architecture that separates non-explosion criteria from boundary-barrier positivity invariance. By bridging microscopic event stoichiometry with macroscopic boundary-degenerate diffusions, this work replaces ad hoc noise constructs with a definitive, mathematically exact template for covariance-consistent and boundary-aware ecological modeling. Citation: Yu J, Wang LS (2026) Beyond diagonal noise: A better predator-prey modeling framework with cross-covariance. PLoS One 21(5): e0350127. https://doi.org/10.1371/journal.pone.0350127 Editor: Md. Kamrujjaman, University of Dhaka, BANGLADESH Received: February 25, 2026; Accepted: May 9, 2026; Published: May 27, 2026 Copyright: © 2026 Yu, Wang. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability: The code that supports the research findings in this work is deposited in the GitHub repository https://github.com/Louis-shuo-wang/Beyond-diagonal-noise-A-better-predator-prey-modeling-framework-with-cross-covariance.git. Funding: The author(s) received no specific funding for this work. Competing interests: The authors have declared that no competing interests exist. Introduction Predator–prey interactions constitute one of the foundational subjects of mathematical ecology. In particular, the Rosenzweig–MacArthur (R–M) model with Holling type II functional response remains a central benchmark for studying coexistence, oscillatory dynamics, and the paradox of enrichment [1–4]. In the deterministic setting, the R–M model exhibits a rich bifurcation structure. As the prey carrying capacity increases, the stable coexistence equilibrium loses stability through a Hopf bifurcation, which produces sustained limit cycles and increases the risk of low-density excursions. This mechanism, in which enrichment destabilizes coexistence, has attracted sustained attention in both ecological theory and applied population biology [5–7]. For finite populations, however, the deterministic ordinary differential equation (ODE) is merely a mean-field approximation. Demographic stochasticity, which arises from randomness in individual birth, death, and predation events, introduces fluctuations that can fundamentally alter quantitative predictions. Even when the deterministic coexistence equilibrium is locally stable, demographic noise can generate quasi-cycles and facilitate noise-induced transitions toward extinction [8,9]. These phenomena are especially pronounced near the Hopf threshold, where the deterministic restoring force weakens and fluctuation amplification becomes large. A rigorous analysis of these dynamics requires a stochastic description that is both mechanistically faithful to the underlying individual-level events and mathematically tractable for asymptotic study. In this paper, we develop precisely such a description. We begin with an integer-valued continuous-time Markov chain (CTMC) formulation of the R–M system that includes explicit Bernoulli-coupled predation–conversion events. From this formulation, we derive a chemical master equation (CME) and, subsequently, a demographic diffusion approximation in the form of a chemical Langevin equation (CLE). Crucially, the drift and covariance of this CLE are strictly dictated by event stoichiometry and channel intensities. A central mathematical consequence of this derivation is that the mechanistically correct diffusion covariance is inherently non-diagonal: coupled predation–conversion events induce a structurally negative predator–prey cross-covariance. This geometric feature is absent from the phenomenological diagonal or multiplicative noise constructions that dominate the literature. We prove that this cross-covariance is a direct algebraic consequence of the stoichiometric coupling between prey loss and predator gain within the same demographic event, rather than an ad hoc modeling embellishment. Beyond the mechanistic derivation of the noise structure, a primary aim of this work is to address a persistent source of ambiguity in stochastic predator–prey modeling. This ambiguity arises from the routine conflation of two distinct diffusion formulations built from the same local stochastic differential equation (SDE) coefficients. We establish that an open-domain formulation on the interior provides the appropriate framework for studying interior coexistence dynamics. This framework captures fluctuation geometry near equilibrium, survival-conditioned statistics, and local spectral structure. Conversely, an absorbed formulation, in which trajectories are frozen upon first contact with the boundary, provides the rigorous framework for extinction-permitting dynamics. This formulation addresses questions about extinction probabilities, mean persistence times, and quasi-stationary distributions, reflecting the irreversibility of CTMC extinction. These are not competing models for the same quantity; they answer fundamentally different biological questions. We argue that their distinction must be stated explicitly to avoid both modeling ambiguity and overreach in well-posedness claims. Building upon the general literature on several topics presented in Literature review section in Appendix, the current work is structured to build our framework systematically, and the main contributions of this paper are as follows: - (C1) Mechanistic derivation of a full-covariance demographic diffusion. Starting from an integer-valued CTMC with Bernoulli-coupled predation–conversion events, we rigorously derive the CME and its associated CLE. This establishes that the diffusion covariance,is strictly dictated by event stoichiometry and channel intensities, replacing phenomenological noise constructs with a purely mechanistic foundation. - (C2) Identification of a structurally negative predator–prey cross-covariance. We mathematically isolate the effect of coupled predation–conversion events, proving they inherently generate a strictly negative cross-covariance, a12(x) < 0, on the interior domain . We demonstrate that drift-matched split-channel closures fail to capture this geometric feature. This result formalizes the limitations of diagonal-noise approximations and identifies regimes where they do not accurately represent the local fluctuation structure. - (C3) Formalization of the open-domain versus absorbed modeling bifurcation. We address a persistent ambiguity in stochastic ecological modeling by explicitly separating two global diffusion formulations built from the same local CME-consistent coefficients. We distinguish an open-domain formulation tailored for interior, survival-conditioned coexistence dynamics from an absorbed formulation designed for extinction-permitting dynamics, ensuring the mathematical framework aligns with the targeted biological observables. - (C4) A two-stage well-posedness architecture for boundary-degenerate diffusions. We develop a tailored analytical framework for the open-domain formulation that separates non-explosion from boundary avoidance. First, we establish maximal strong solvability without interior blow-up before boundary contact via a Lyapunov anti-explosion condition. Second, we establish positivity invariance as a distinct property dependent on an additional boundary-barrier condition, preventing the tacit assumption of survival from being embedded in baseline well-posedness claims. Fig 1 presents the mechanistic roadmap and modeling bifurcation for the R–M model. The structure of this work is as follows. Methods section details the mathematical framework. Deterministic R–M backbone subsection records the deterministic R–M backbone. Mechanistic CTMC and CME subsection and CME-consistent covariance structure subsection develop the CTMC/CME formulation and the full-covariance diffusion derivation. Domain formulations subsection formalizes the open-domain and absorbed formulation distinction. Well-posedness subsection develops the two-stage well-posedness framework. Results section presents the numerical simulations that verify the boundary-aware, covariance-consistent modeling structure. Discussion section discusses implications and scope. Appendices provide full derivations of the CME-to-diffusion reduction and complete proofs of the well-posedness theorems. The diagram summarizes the paper’s logic from the deterministic R–M backbone and Hopf organization, through a mechanistic CTMC/CME derivation of a full-covariance CLE, to the structural negative predator–prey cross-covariance induced by coupled predation–conversion. It then highlights the explicit modeling bifurcation between open-domain and absorbed diffusion formulations built from the same local coefficients, and the corresponding two-stage well-posedness architecture for the open-domain SDE. Methods Deterministic R–M backbone This section records the deterministic R–M structure used as the background for the stochastic construction developed later. We consider the nondimensional R–M predator–prey system with Holling type II predation [10]: (1)where N and P denote prey and predator densities, respectively, and k > 0, m > 0, c > 0 are the (scaled) carrying-capacity, predation/assimilation, and predator-mortality parameters. The biologically relevant closed state space is the positive quadrant Q, while the open quadrant (2)will be the ambient domain for the open-domain diffusion formulation from Domain formulations subsection to Well-posedness subsection in Methods section. The vector field in (1) is smooth on Q (indeed on ). We work with the standard reduced nondimensional R–M form in which the predator conversion factor is not written explicitly in the term in the P equation in (1). In the mechanistic stochastic derivation (Mechanistic CTMC and CME subsection in Methods section), this corresponds to the normalization e = 1 when matching the reduced deterministic parameterization used throughout the analysis. System (1) has three equilibria: (3)(4)(5)The coexistence equilibrium K3 lies in D if and only if (6)The first condition ensures that the predator per-capita growth term can become nonnegative for sufficiently large prey density, while the second condition is exactly the positivity condition . In the present paper, K3 serves as the deterministic coexistence reference state for the mechanistic stochastic modeling and the boundary-aware diffusion formulations. Let The Jacobian matrix is (7)At the coexistence equilibrium K3, one has whenever , so local stability is determined by the trace. A direct calculation gives the classical enrichment-driven Hopf threshold [10] (8)Hence K3 is locally asymptotically stable for k < kH and unstable for k > kH, with a Hopf bifurcation at . We primarily use (8) to organize parameter regimes (9)Fig 2 illustrates the bifurcation diagram in the carrying capacity, and the phase portraits for the two regimes (stable coexistence) and (unstable oscillation), showing the Hopf bifurcation. (a) Bifurcation diagram over the carrying capacity k. The coexistence equilibrium K3 loses stability at the Hopf threshold kH, partitioning the parameter space into the stable regime and the oscillatory regime . Solid curves denote the extrema of the stable limit cycle, and dashed lines denote the unstable equilibrium. (b) Phase portrait in the stable regime (k < kH). The predator nullcline (vertical) intersects the prey nullcline (humped) to the right of its peak, resulting in trajectories spiraling inward to a stable coexistence state K3. (c) Phase portrait in the oscillatory regime (k > kH). Enrichment shifts the intersection to the left of the peak, destabilizing K3 and generating a stable limit cycle that drives large-amplitude, low-density excursions. Mechanistic CTMC and CME We now formulate an integer-valued CTMC for the R–M predator–prey system and record the associated CME. Let denote prey and predator counts in a well-mixed population of system size (e.g., habitat volume or scaling parameter), and define the corresponding density process Because the exact CTMC is integer-valued, all jump increments in X must lie in . In particular, a compact predation–conversion increment of the form with is not an admissible exact CTMC jump unless . We therefore formulate the CTMC using integer-valued channels (with Bernoulli conversion), and only later use an effective increment as a density-level diffusion closure in CME-consistent covariance structure subsection in Methods section. This distinction is essential for mechanistic consistency. We use the following five integer-valued event channels: This is equivalent to a predation encounter process followed by Bernoulli conversion with success probability . The stoichiometric increments are (10)and the corresponding stoichiometric matrix is (11)We adopt density-dependent (Kurtz-type) propensities (12)with intensity functions (13)(14)Let . Then the deterministic drift induced by the CTMC is (15)This is the R–M drift with explicit conversion efficiency e. In the current mechanistic CTMC/CME formulation, enters the predator gain term but not the prey loss term. To maintain exact consistency with the deterministic backbone (1) used in Deterministic R–M backbone subsection in Methods section, we perform the main stochastic analysis under the normalization The mechanistic derivations in this section are written for general e where useful (to expose the origin of covariance structure), and the specialization e = 1 is imposed when matching the reduced R–M parameterization used from Domain formulations subsection in Methods section to Results section. For later comparison, it is useful to note that one can construct an alternative integer-valued CTMC with separate channels for prey removal and predator reproduction, using increments and with rates proportional to and , respectively. This split-channel construction reproduces the same deterministic drift (15) but changes the count-level and diffusion-level covariance structure. In particular, it removes the same-event predator–prey coupling in the predation contribution and will therefore serve as a diagonal-noise comparator in CME-consistent covariance structure subsection in Methods section, not as the primary mechanistic baseline. We formulate the CME. Let , . For a finite family of channels with increments and propensities , the CME is (16)In the present Bernoulli-coupled model, K = 5, the increments are those in (10), and the propensities are given by (12)–(14). We use (16) primarily as the exact count-level probabilistic description from which moment identities and the diffusion covariance inherit their structure. The associated CTMC generator acting on test functions is (17)As recalled in From the CME to the diffusion approximation section in Appendix, this yields both the count-scale drift increment and the exact infinitesimal covariance (quadratic variation density) in count variables: (18)Equivalently, in stoichiometric form, (19)for the chosen channel representation . This identity is exact at the CTMC level and already shows that second-order structure depends on event stoichiometry, not only on the drift. The key implication of (18)–(19) is that, under density scaling with , the diffusion covariance in the demographic-noise approximation is inherited directly from the count-level event structure. In CME-consistent covariance structure subsection in Methods section, we make this explicit through the unified formula and use it to compare Bernoulli-coupled, effective coupled , and split-channel closures. This is where the structural negative predator–prey cross-covariance appears transparently. CME-consistent covariance structure Let S denote a chosen stoichiometric matrix and f(x) the associated vector of density-level intensities. Under the density-dependent scaling , the chemical Langevin diffusion approximation for the density process takes the Itô form (20)with drift (21)and diffusion covariance (22)Equation (22) is the density-scale counterpart of the exact CTMC infinitesimal covariance identity in counts (Mechanistic CTMC and CME subsection in Methods section); see From the CME to the diffusion approximation section in Appendix for the derivation. To isolate the predation contribution across different closures, we write (23)Below we compare three drift-compatible predation closures (up to parameterization): (i) the exact Bernoulli-coupled CTMC mechanism, (ii) an effective diffusion-level coupled channel , and (iii) a split-channel comparator. All three closures reproduce the same first-order drift, but they differ at the level of quadratic variation and therefore in their diffusion covariance structure. For the Bernoulli-coupled predation channels and from Mechanistic CTMC and CME subsection in Methods section, the density-level predation covariance contribution is (24)Hence, the off-diagonal predation term is (25)In particular, for and e > 0, this term is strictly negative. Biologically, it records the same-event negative correlation created by prey removal and predator gain in a predation encounter with successful conversion. For compact diffusion notation, one may use a single effective predation channel at the density level: (26)This effective closure preserves the deterministic drift, but it aggregates the discrete predation–conversion mechanism into one diffusion jump. As a result, the predator fluctuation intensity is altered at second order. This is not an exact integer-valued CTMC jump when , but it is a valid diffusion-level closure consistent with the drift. The corresponding predation covariance contribution is (27)Therefore, (28)Thus, the effective coupled closure preserves the same predation-induced cross-covariance sign and magnitude as the exact Bernoulli-coupled mechanism, but it does not preserve the predator variance. Indeed, the exact Bernoulli-coupled closure gives a predator variance contribution proportional to e, whereas the effective single-channel closure gives e2. This is a direct consequence of collapsing a discrete success/failure conversion process into one diffusion jump, since diffusion covariances scale with the square of the effective jump size. The distinction is retained for the general , while the numerical baseline e = 1 makes them coincide. Consider the split-channel predation representation with Then the predation-related covariance contribution is (29)Hence (30)This representation is drift-compatible with the coupled mechanisms but removes the same-event predator–prey fluctuation coupling. For that reason, we treat it as a diagonal-noise comparator, not as the default mechanistic closure. Proposition 0.1 (Structural predation cross-covariance under coupled versus split closures). Let and . For the predation intensity , the predation contribution to the diffusion covariance satisfies: - (i) under the exact Bernoulli-coupled closure, - (ii) under the effective coupled (−1, e) closure, - (iii) under the split-channel closure, Proposition 0.1 implies that the drift equivalence does not imply covariance equivalence. Different channel representations can produce the same deterministic drift while yielding different diffusion covariances. Indeed, the drift depends on the first stoichiometric moments whereas the covariance depends on the second stoichiometric moments Thus, channel decompositions that are indistinguishable at the ODE level may differ at the diffusion level. This is the precise mechanism behind the distinction between coupled predation closures (which produce predator–prey cross-covariance) and split-channel closures (which can suppress it). Within the class of drift-compatible predation closures considered above, the sign and presence of predator–prey cross-covariance at the diffusion level is determined by whether prey loss and predator gain occur in the same event channel. In particular, the negative off-diagonal term is not an ad hoc modeling embellishment; it is a stoichiometric consequence of coupled predation–conversion. Fig 3 shows the microscopic state transitions with arrows indicating the stoichiometric increments, and the macroscopic covariance ellipses with negative correlation tilting the fluctuation geometry. (a) Microscopic state transitions on the count lattice. The exact Bernoulli-coupled mechanism strictly requires a diagonal jump vector (crimson) to represent a single predation–conversion event. The split-channel comparator falsely separates this into two independent, orthogonal jumps. (b) Macroscopic covariance geometry at the density scale. The diagonal jump in the coupled mechanism inherently tilts the local fluctuation geometry, generating a strictly negative predator–prey cross-covariance (a12 < 0). In contrast, the phenomenological split-channel closure produces an axis-aligned covariance ellipse (a12 = 0), failing to capture the true geometric structure of the demographic noise despite maintaining mean-field drift equivalence. Remark 0.2 (Bernoulli-coupled exact CTMC vs. effective closure). The exact Bernoulli-coupled CTMC closure and the effective diffusion closure agree in deterministic drift and in the predation-induced cross-covariance magnitude but they differ in the predator-variance contribution: In this paper, the effective closure is used as a compact full-covariance diffusion representation, while the Bernoulli-coupled CTMC remains the exact integer-valued mechanistic foundation. For later reference, it is sometimes convenient to decompose the full covariance into non-predation and predation contributions: where collects prey birth, prey competition death, and predator death channels, and is one of Since the base channels in the present construction contribute only diagonal terms, the sign and magnitude of the full off-diagonal covariance a12(x) are determined entirely by the predation component. Therefore, Proposition 0.1 directly controls the sign of the full predator–prey cross-covariance in the coupled closures. Since a(x) is symmetric positive semidefinite on , it admits (non-unique) factorizations (31)Two useful choices are: - (i) Event-based factorization that uses K independent Brownian drivers (one per channel): which preserves channel semantics and is convenient for interpreting noise sources; - (ii) Minimal-dimensional square root (e.g., Cholesky): choose a measurable 2 × 2 square root satisfying , which is convenient for numerical simulation and for the formulation of the open-domain SDE. These choices are equivalent in law at the level of the diffusion process defined by (20), provided they generate the same covariance matrix a(x). Assumption 0.3 (Diagonal-noise as a controlled approximation: small relative cross-covariance). Let be a state region of interest. A diagonal-noise approximation is considered only in regimes where (32)That is, diagonalization is treated as a diagnostic approximation rather than a default mechanistic model. Remark 0.4 (Relating small cross-covariance to ecological parameters). Under the coupled closures, the cross-covariance term |a12(x)| is driven directly by the predation event rate, scaling as . Assuming a standard functional response (with maximum predation rate m), the condition in (32) translates directly into specific ecological parameter regimes where independent demographic processes dominate over coupled interaction events. A modeler might find a diagonal approximation accidentally valid in the following scenarios: - (i) Low conversion efficiency (): The predator requires consuming many prey items to produce a single offspring. This mechanically decouples the noise of a prey death from a predator birth, significantly shrinking the cross-covariance relative to the individual variances. - (ii) Weak maximum predation rate (m): If the predation rate m is small relative to the intrinsic growth or natural mortality parameters of the species, the interaction variance driven by fpred(x) is dwarfed by the diagonal demographic noise terms a11(x) and a22(x). - (iii) state restriction (k): If environmental carrying capacity k restricts N or P to remain uniformly small on , the overall predation rate fpred(x) is suppressed. When these conditions fail (e.g., highly efficient predators, strong predation pressure m, and limited resources k), diagonal-noise closures are expected to misrepresent second-order statistics, in particular Cov(N, P) and the geometry of local fluctuations. In the current work, the CTMC/CME formulation provides a microscopic description, whereas the diffusion (CLE) approximation yields a macroscopic limit. Thus, our framework establishes a principled bridge between scales, aligning conceptually with hybrid multiscale modeling approaches that couple discrete events with continuum dynamics [11,12]. Domain formulations The open-domain formulation is the diffusion model posed on the interior D, with initial condition : (33)Under the full-covariance coupled closure, the predation mechanism induces a negative off-diagonal term a12(x)<0 for , as shown in CME-consistent covariance structure subsection in Methods section. At this stage, (33) is a local SDE statement on an open set; the issue of boundary attainability is deferred to Well-posedness subsection in Methods section. To model extinction-permitting dynamics, one uses an absorbed formulation based on the same local coefficients. Let (34)denote the first exit (boundary hitting) time from the interior for a local solution of (33). The absorbed process is then defined by freezing at the first boundary hit: (35)Biologically, the boundary corresponds to loss of strict coexistence: N = 0 (prey extinction), P = 0 (predator extinction), or (0,0) (total extinction). The freezing rule in (35) is the diffusion-level analogue of the irreversibility of extinction in the underlying finite-population CTMC. Remark 0.5 (Why these are different questions (and why the distinction matters)). The open-domain and absorbed formulations differ in at least three ways. First, they target different observables: interior fluctuation structure versus extinction/permanence statistics. Second, they differ in time-horizon emphasis: open-domain analyses are often local or survival-conditioned, whereas absorbed analyses directly encode rare-event accumulation and eventual boundary contact. Third, they correspond to different analogies with the exact CTMC: the absorbed formulation mirrors CTMC extinction irreversibility, while the open-domain formulation serves as a mathematically convenient interior approximation for coexistence dynamics. The distinction between the open-domain and absorbed formulations is not merely a technical distinction about how to close an SDE at the boundary. It corresponds to two different classes of biological observables. The open-domain formulation is appropriate when the empirical target is the fluctuation structure of populations that remain extant over the observational window. Examples include the amplitude and phase of predator–prey cycles, local covariance geometry around coexistence, quasi-cycles, stochastic amplification near a stable equilibrium, and survival-conditioned summaries of population variability. Such observables are central in empirical studies of population cycles and consumer–resource time series, where the goal is to explain persistent oscillatory abundance patterns rather than to estimate the first time at which one species disappears [13–15]. By contrast, the absorbed formulation is required when the empirical or management target is a boundary-hitting observable. In conservation applications, population viability analysis commonly estimates probabilities of extinction or quasi-extinction over a fixed horizon, distributions of persistence time, and mean or median time to extinction [16–18]. These quantities depend explicitly on whether and when a trajectory reaches a low-abundance or extinction boundary. A model that is conditioned to remain in the interior cannot answer such questions without changing the target quantity. Conversely, an absorbed model is not designed to describe the stationary or quasi-stationary geometry of surviving interior fluctuations after extinction events have been removed or conditioned away. Thus, the modeling bifurcation reflects a biological difference in the question being asked. If the question concerns coexistence fluctuations among surviving populations, the relevant object is the interior law of the diffusion, often interpreted conditionally on non-extinction over the time scale of observation. If the question concerns extinction risk, persistence time, or quasi-extinction thresholds, the relevant object is the absorbed process. The same local drift and covariance coefficients may therefore support two globally different stochastic models because the empirical observables are different. Fig 4 contrasts the representative trajectories in the open domain under the open-domain formulation, and those in the first quadrant under the absorbed diffusion formulation. Both panels simulate the exact same local chemical Langevin drift and covariance parameters (, a fluctuation-amplified regime) originating from the identical initial state x0. (a) The absorbed formulation, where trajectories are permanently frozen upon first boundary contact (), marked by red stars. This is the mathematically appropriate framework for evaluating extinction probabilities and mean persistence times. (b) A numerical proxy for the open-domain formulation, where trajectories are conditioned to survive and continuously explore the interior domain . This framework captures the geometry of quasi-stationary interior fluctuations but must not be conflated with the extinction-permitting dynamics strictly captured in (a). Well-posedness We consider the Itô SDE (33) on the open domain D. For the R–M demographic diffusion considered here, the drift b and covariance a are smooth on D, while a chosen factorization may fail to be globally Lipschitz near because demographic noise intensities vanish as and/or . Nevertheless, on every precompact set , the coefficients are locally bounded and (for the factorizations used in this paper) locally Lipschitz. This is the natural setting for localization and Lyapunov arguments. Define the first exit (boundary hitting) time from the interior by (34) and freeze xt according to the rule in (35). We put the detailed proof for Theorem 0.6 and Proposition 0.7 in Proofs for well-posedness section in Appendix. Theorem 0.6 (Maximal strong solution and no interior explosion before boundary hitting). Assume: (A1) Local regularity on compacts. For every precompact set , the coefficients b and are Lipschitz and bounded on K. (A2) Lyapunov control at infinity. There exist and a constant C > 0 such that (36)where and whenever with . Then for every , the SDE (33) admits a unique maximal strong solution with lifetime , and (37)Equivalently, loss of well-posedness cannot occur in the interior of D; if , the only possible obstruction is boundary hitting (or boundary approach), i.e., (38)Proposition 0.7 (Positivity invariance under a boundary-barrier Lyapunov condition). Assume the hypotheses of Theorem 0.6. In addition, suppose there exists a function such that: (B1) Boundary blow-up (barrier property). (39)(B2) Localized generator bound. There exists such that for every n ≥ 1 and all , (40)Then (41)Consequently, the maximal strong solution is global and remains in D for all times: (42)For better understanding, we apply Theorem 0.6 to our model. We note that the computations below are applicable to all cases, while we only take e = 1 to remain consistent with the current framework. Here, We take the Lyapunov function as where and . Then, it is easy to see Here, By elementary inequalities, e.g., it is easy to find a constant that satisfies Although the choice is primarily motivated by analytical convenience, it admits a natural ecological interpretation. Since V(x) is a strictly increasing function of , it can be viewed as a radial abundance index for the predator–prey system, measuring the overall magnitude of the population state (N,P). In particular, is comparable to the squared Euclidean distance from the extinction state (0,0), and hence V(x) provides a smooth proxy for the distance from extinction. Large values of V(x) correspond to configurations where at least one population is large, while small values correspond to states close to extinction. We emphasize that V(x) is not intended to represent a conserved ecological quantity such as total biomass N + P. Instead, it plays the role of a coercive control function that penalizes large excursions of either species. On comparable density scales, it can nevertheless be interpreted as a surrogate for aggregate biomass, while retaining the regularity and growth properties required for the Lyapunov analysis. Fig 5 visualizes the mathematical mechanisms for the prevention of interior blow-up and boundary hitting that underlies Theorem 0.6 and Proposition 0.7. (a) The anti-explosion Lyapunov function V∞(x) (Theorem 0.6) acts as a geometric bowl that grows radially to infinity, strictly confining the stochastic trajectory from escaping to infinite states without restricting its boundary approach. (b) The boundary-barrier Lyapunov function Vb(x) (Proposition 0.7) acts as a steep geometric cliff. The surface is extremely flat within the interior domain , permitting natural quasi-stationary fluctuations, but spikes asymptotically to infinity precisely at the boundary faces N = 0 and P = 0, theoretically repelling the trajectory and ensuring positivity invariance almost surely. Demographic diffusions often invite comparison with one-dimensional square-root diffusions, where Feller boundary classification [19,20] gives explicit criteria for boundary attainability in terms of drift–diffusion balance. That intuition is useful but incomplete here. In the present two-dimensional predator–prey diffusion: - (i) drift and diffusion depend on both coordinates via nonlinear interaction terms; - (ii) the covariance need not be diagonal (indeed, the mechanistic full-covariance closure produces a12(x)<0 in D); - (iii) the boundary has multiple faces (N = 0, P = 0) with distinct biological meanings. Accordingly, a literal one-dimensional Feller criterion does not transfer directly to the coupled two-dimensional setting. The barrier-based Lyapunov framework above is designed precisely to handle this multivariate, boundary-degenerate structure without forcing an oversimplified 1D analogy. A step-by-step guide for applyivng the framework to other ecological systems Although the R–M predator–prey system provides the main worked example in this paper, the framework is not specific to two species or to Holling type II predation. The essential object is not the deterministic vector field alone, but the event-level representation from which the drift and covariance are derived. We therefore summarize the transferable procedure as a five-step workflow. Step 1: Identify the integer-valued demographic events. The starting point is a list of biologically meaningful events that change integer population counts. These may include births, natural deaths, competition-induced deaths, predation events, conversion or reproduction events, immigration, harvesting, or disease transmission. Each event must be formulated as an admissible jump in the count process. For a system with species counts each event changes the state by an integer vector. This step is where biological mechanism enters the model most directly. In particular, if one ecological event simultaneously decreases one species and increases another, such as predation with conversion, then this coupling should be represented as a single joint event rather than as two independent noise sources. For example, in a three-species food chain with basal prey X1, intermediate predator X2, and top predator X3, one may include prey birth, prey competition death, intermediate-predator death, top-predator death, predation of X1 by X2, and predation of X2 by X3. If predation is coupled to conversion, the corresponding count-level jumps include vectors such as These joint events automatically encode negative covariance between resource and consumer fluctuations. In a two-species competition model, by contrast, typical demographic events might include intrinsic births , , intrinsic deaths , , and competition-induced deaths , with rates depending on both species densities. Unless a single event changes both species simultaneously, the demographic covariance may remain diagonal even though the drift is dynamically coupled through the rates. Step 2: Define the stoichiometric matrix and channel intensities. Once the event list is fixed, assign to each event k a stoichiometric increment and collect the increments in the stoichiometric matrix Next, specify density-dependent propensities where fk(x) is the density-level event intensity and is the system-size parameter. This step separates the biological event mechanism from the scaling approximation. The increments describe what happens when an event occurs; the functions fk describe how often those events occur. For the three-species food-chain example, a coupled predation event of the intermediate predator on the basal prey may have or another functional response appropriate to the system. A top-predator conversion event may similarly have The specific functional response can change, but the covariance construction is unchanged. Step 3: Derive both the drift and the covariance from the same event representation. The deterministic drift and the demographic diffusion covariance are then obtained from the same stoichiometric data: This is the central methodological step. It ensures that the stochastic model is not obtained by adding phenomenological noise to an ODE, but by carrying the integer-valued event structure to the diffusion approximation. This step is also where covariance geometry is determined. A joint event with increment contributes to the covariance matrix. Hence the basal prey and intermediate predator acquire a negative cross-covariance contribution. Similarly, a joint event contributes a negative cross-covariance between the intermediate and top predator. Thus, in a food chain, the full covariance matrix is generally banded or structured according to trophic event couplings, rather than diagonal. For a competition model, if interspecific competition affects only the rates of single-species deaths, then the drift can be strongly coupled while the demographic covariance remains diagonal. This contrast illustrates the main principle of the present paper: drift coupling and noise coupling are different objects. Drift coupling is controlled by the dependence of fk(x) on the full state x, whereas covariance coupling is controlled by whether the same event channel changes multiple species simultaneously. Step 4: Specify the biologically appropriate domain formulation. After the local coefficients b(x) and a(x) have been derived, the modeler must specify the global interpretation of the diffusion. For a d-species system, the natural interior domain is The open-domain formulation is appropriate when the target observables concern surviving coexistence dynamics: local covariance near an equilibrium, quasi-cycles, stochastic amplification, spectral structure, phase relationships, or survival-conditioned fluctuations. In contrast, the absorbed formulation is appropriate when the target observables are extinction probabilities, quasi-extinction probabilities, persistence-time distributions, or mean time to extinction. In conservation and population viability contexts, such boundary-hitting quantities are often the primary empirical or management target. For a three-species food chain, different boundary faces have different biological meanings. The face x3 = 0 corresponds to loss of the top predator, x2 = 0 corresponds to collapse of the intermediate consumer, and x1 = 0 corresponds to basal-resource extinction. Absorbing all boundary faces may be appropriate when any species loss is treated as irreversible over the time scale of interest. Alternatively, one may define face-specific absorbing rules if the biological question distinguishes top-predator loss from total food-chain collapse. For a competition model, the boundary face xi = 0 represents extinction of competitor i, and the absorbed formulation is the natural framework for coexistence probability, exclusion probability, or time-to-loss of a focal species. Step 5: Verify well-posedness and boundary behavior. Finally, the modeler must verify that the resulting SDE is mathematically well posed for the intended formulation. The verification should be separated into two stages. First, local regularity on compact subsets of D and an anti-explosion Lyapunov function establish maximal strong solvability and rule out blow-up in the interior. A typical choice is a coercive abundance function such as with , followed by a generator estimate of the form This step prevents unbounded population excursions from causing loss of the solution before any boundary is reached. Second, if the intended model is an open-domain survival formulation, one must separately verify boundary non-attainment. This typically requires a boundary-barrier function Vb satisfying together with a localized generator bound This step should not be assumed automatically. It is a biological and mathematical modeling choice: if extinction-permitting dynamics are the target, then boundary hitting is not a pathology but the event of interest, and the absorbed formulation should be used instead. This five-step workflow turns the present R–M analysis into a general method. The modeler first builds the integer-valued event system, then derives the diffusion covariance from stoichiometry, then selects the appropriate global domain interpretation, and finally verifies the well-posedness conditions required by that interpretation. The same procedure applies to predator–prey systems, food chains, competition models, host–parasitoid systems, epidemic–ecological models, and higher-dimensional food webs. What changes from model to model is the event list and the biological interpretation of the boundary; the covariance construction and domain logic remain the same. Results This section provides illustrative numerical experiments supporting the structural results established from Mechanistic CTMC and CME subsection to Well-posedness subsection in Methods section. The emphasis is on visualizing the consequences of (i) the mechanistically induced diffusion covariance structure, in particular the predator–prey cross-covariance term, and (ii) the distinction between open-domain and absorbed formulations built from the same local coefficients. All simulations use the nondimensional R–M system (1) with the Bernoulli-coupled predation–conversion mechanism of Mechanistic CTMC and CME subsection in Methods section. Unless otherwise stated, figures are generated from the exact CTMC (Gillespie SSA) or from Euler–Maruyama discretizations of the density-level diffusion under the covariance closures introduced in CME-consistent covariance structure subsection in Methods section. Throughout, we fix the baseline parameters (43)and vary the carrying capacity k and system size as indicated. Under (43), the Hopf threshold (8) is Covariance ellipses: Coupled versus split-channel closures CME-consistent covariance structure subsection in Methods section established that drift-compatible closures can produce different diffusion covariances, and in particular that coupled predation–conversion induces a negative off-diagonal entry a12(x)<0, while the split-channel comparator yields (Proposition 0.1). We visualize this local geometric difference at the coexistence equilibrium K3 through covariance ellipses. Fig 6 displays and covariance ellipses at K3 under three closures: (i) the effective coupled closure (27), (ii) the exact Bernoulli-coupled closure (24), and (iii) the drift-compatible split-channel comparator (29). The two coupled closures produce tilted ellipses, reflecting the negative predator–prey cross-covariance, whereas the split-channel ellipse is axis-aligned. This provides a direct geometric signature of the structural sign result. Left: effective coupled . Center: exact Bernoulli-coupled. Right: split-channel (diagonal comparator). Solid and dashed ellipses show and contours, respectively. The negative cross-covariance tilts the coupled-closure ellipses, while the split-channel ellipse is axis-aligned. The panels also annotate the normalized cross-covariance ratio As predicted by CME-consistent covariance structure subsection in Methods section, the Bernoulli-coupled and effective closures share the same off-diagonal term (hence the same sign and similar tilt direction), while differing in the predator variance entry a22; see Remark 0.2. Cross-covariance ratio across parameter space To assess when a diagonal-noise approximation may be acceptable in practice (cf. Assumption 0.3 and Remark 0.4), we evaluate the normalized cross-covariance ratio across the (m,k) plane with c = 0.4 and e = 1 fixed. Fig 7 shows under the effective coupled closure. The ratio is uniformly non-positive across the parameter region, in agreement with Proposition 0.1, and becomes more negative when m,k increase. The Hopf threshold is overlaid to indicate the partition of (m,k) parameter space into and . The dashed curve is the Hopf threshold . White regions indicate parameter values for which coexistence is infeasible. Absorbed versus open-domain boundary behavior Domain formulations subsection in Methods section introduced two distinct global interpretations built from the same local coefficients: an absorbed process for extinction-permitting dynamics and an open-domain formulation for interior dynamics. To visualize the modeling consequences of this distinction, we consider a parameter regime with amplified fluctuations (, k = 5.5) and simulate representative sample paths over a finite horizon. Fig 8(a,d) uses the absorbed formulation, in which trajectories are frozen at first boundary contact, yielding an empirical extinction-time distribution. Fig 8(b) shows a positivity-corrected Euler–Maruyama proxy intended to visualize the interior-dynamics perspective associated with the open-domain formulation: when a discrete-time step produces a small sign violation, a numerical clipping near zero is applied to suppress step-size-induced boundary crossing behaviors. This numerical device is used only for visualization and should not be conflated with the continuous-time well-posedness statements proved in Well-posedness subsection in Methods section. (a) Absorbed formulation phase-plane paths. (b) Positivity-corrected Euler–Maruyama proxy for the open-domain interior-dynamics perspective. (c) Boundary-distance diagnostic for absorbed paths. (d) Empirical extinction-time distribution from n = 300 absorbed simulations. The contrast is nevertheless instructive. The absorbed formulation produces a nontrivial distribution of extinction times, while the open-domain proxy emphasizes continued interior fluctuation dynamics. Thus, although the two formulations share the same local drift and covariance in the interior, they answer different modeling questions and produce different observables. Parameter-regime visualization relative to the Hopf threshold Finally, we provide a brief visual comparison of stochastic trajectories across the parameter partition induced by the Hopf threshold (8). Fig 9 compares full-covariance diffusion simulations (with deterministic trajectories overlaid) for , , and . Left: (inside ). Center: (threshold). Right: (inside ). Top row: prey time series (diffusion vs. ODE). Bottom row: phase-plane trajectories. The purpose of this figure is primarily organizational: it illustrates that the partition into and corresponds to visibly different fluctuation regimes (stable-coexistence fluctuations, threshold-amplified oscillatory fluctuations, and noise-perturbed cycling, respectively). Taken together, these numerical experiments reinforce the structural conclusions developed in the analytical sections. The covariance-ellipse visualizations provide a direct geometric manifestation of the mechanistically induced negative cross-covariance; the parameter-space scan shows that this structural sign persists across the parameter sweeps; and the boundary comparisons clarify that absorbed and open-domain formulations, though locally identical, encode distinct modeling objectives. The simulations, therefore, serve not as independent evidence but as visual diagnostics of the theoretical architecture: they illustrate how event-channel design, covariance structure, and boundary interpretation propagate from microscopic assumptions to macroscopic stochastic behavior. Discussion This paper develops a mechanistically consistent stochastic framework for the R–M predator–prey model. It demonstrates that diffusion-level drift and covariance must be rigorously derived from an integer-valued CTMC/CME description rather than imposed directly at the SDE level. Ultimately, the value of this work lies not in proposing yet another stochastic variant of the R–M model, but in establishing a mathematically rigorous bridge linking: (i) microscopic event stoichiometry to macroscopic demographic covariance, (ii) open-domain survival to absorbed-boundary extinction, and (iii) boundary-degenerate local coefficients to a globally tailored well-posedness architecture. Mechanistic covariance structure as a primary result (Main result 1) Our first principal contribution isolates the macroscopic consequences of microscopic event coupling. Starting from the Bernoulli-coupled predation–conversion CTMC, the diffusion covariance is strictly inherited through the CME-consistent identity Within this framework, the predator–prey cross-covariance is not a tunable modeling embellishment: it is an unavoidable structural output of the stoichiometric coupling between prey loss and predator gain during predation–conversion events. In particular, the coupled closures systematically studied in CME-consistent covariance structure subsection in Methods section strictly produce a negative predation-induced off-diagonal term. This rigorously exposes the fallacy that drift equivalence implies covariance equivalence, demonstrating that diagonal-noise closures must be relegated to regime-specific approximations rather than treated as default mechanistic baselines. Boundary-aware formulation distinction as a modeling result (Main result 2) The second principal result is the formal bifurcation of two global diffusion formulations built from the exact same local coefficients. One is an open-domain formulation on , used strictly for interior, survival-conditioned dynamics. The other is an absorbed formulation, which allows for extinction-permitting dynamics. While routinely blurred or entirely ignored in the stochastic ecological literature, this distinction is both mathematically profound and biologically consequential. The open-domain formulation provides the precise analytical arena for coexistence fluctuations and local spectral geometry. Conversely, the absorbed formulation is the exclusive framework for evaluating extinction probabilities and mean persistence times, mirroring the irreversibility of CTMC extinction. By stating this bifurcation explicitly, our framework eradicates a pervasive source of modeling ambiguity. Two-stage well-posedness architecture as an analytical result (Main result 3) The third principal result is the customized two-stage well-posedness architecture for the open-domain SDE. Well-posedness subsection in Methods section deliberately separates: - (i) maximal strong solvability together with non-explosion in the interior before boundary hitting, and - (ii) positivity invariance under an additional boundary-barrier Lyapunov condition. This separation is not a mere technical bookkeeping device; it is a vital architectural safeguard. It explicitly aligns the mathematical theorems with the biological modeling distinction formalized in Domain formulations subsection in Methods section. Positivity invariance is a property that may hold in a specific open-domain analytical regime, but it is not a default mathematical right one can claim when the intended model is extinction-permitting. This formulation thereby avoids analytical overreach and provides a pristine mathematical template for other boundary-degenerate population diffusions. Scope and limitations We contextualize these contributions within several recognized limitations. First, the diffusion model remains a second-order approximation to the exact CTMC/CME. It is therefore most robust in regimes where the system-size parameter is sufficiently large and event counts over observational windows are strictly non-sparse. For severely small populations, integer-valued effects and immediate boundary absorption dominate, rendering exact CTMC/SSA descriptions as the mandatory baseline. Second, the present paper emphasizes exact structural derivation and well-posedness architecture rather than quantitative diffusion error bounds; no uniform-in-time approximation bound in is pursued here. Third, the open-domain analysis is formulated at the abstract level of local coefficients and Lyapunov or barrier criteria. Verifying these conditions in highly specific or extreme parameter regimes may require bespoke analytical constructions beyond the generic bounding techniques presented. Generalizations and portability of the framework Crucially, the analytical framework developed here is highly portable and extends well beyond the R–M model. The core covariance identity depends universally on event stoichiometry and intensity functions, completely agnostic to the specific Holling type II form. Thus, this exact blueprint can be deployed for alternative predator–prey mechanisms (e.g., ratio-dependent interactions) and seamlessly scaled to high-dimensional multi-species food-web models. More broadly, the integration of biological mechanisms, rigorous mathematical modeling, and threshold-driven organization reflects a portable analytical paradigm that extends to diverse biological systems beyond ecological interactions [21,22]. In every such extension, the same three modeling questions reappear. First, which event couplings mandate nontrivial cross-covariances? Second, which specific boundary faces represent biologically distinct extinction cascades? Third, which open-domain versus absorbed interpretations align with the target observables? This paper provides the foundational, reusable template for answering those questions explicitly. Conclusions The central thesis of this paper is that demographic covariance in predator–prey diffusion models must be treated as a mechanistic output, not a stylistic choice. For the R–M system, predation–conversion coupling structurally mandates a negative predator–prey cross-covariance. This feature interacts fundamentally with boundary interpretation and stochastic well-posedness. We derive the diffusion directly from an integer-valued CTMC model, explicitly separate open-domain from absorbed formulations, and formalize a bespoke two-stage well-posedness architecture. Together, these steps establish a mathematically rigorous foundation for covariance-consistent, boundary-aware stochastic ecological modeling. Appendix Literature review We present the necessary literature on which the current study is based. The present work connects several partially overlapping bodies of literature: (i) stochastic predator–prey modeling and its mechanistic foundations, (ii) demographic-noise amplified population oscillation, (iii) CME and chemical Langevin methods in population biology, and (iv) well-posedness and boundary analysis for diffusion processes arising in ecology and mathematical biology. We review each part in turn and identify the specific gaps that motivate the contributions of this paper. Stochastic predator–prey models: Noise structures and their origins. The introduction of noise into predator–prey systems has a long history. Early works by May [3] and Gard [23] studied stochastic Lotka–Volterra-type equations with environmental noise. They establish that random fluctuations can qualitatively alter or preserve persistence and coexistence depending on noise intensity. Subsequently, a large literature developed stochastic R–M models and related systems, using both Itô and Stratonovich SDE formulations with various noise structures [24–29]. A common modeling approach in this literature is to introduce multiplicative noise directly at the SDE level, typically by adding separate, independent Wiener processes to the prey and predator equations: or variants thereof [30–34]. Such constructions are mathematically tractable and allow for clear persistence and extinction analyses. However, they remain fundamentally phenomenological rather than mechanistic because the noise intensities and are free parameters not derivable from event stoichiometry, and the independence of the two noise sources is imposed rather than justified. In particular, a predation event that simultaneously removes a prey individual and possibly adds a predator individual counts as a single demographic event. Its natural diffusion-level representation must involve a coupled, non-diagonal fluctuation increment rather than two independent ones. This mechanistic structure is inherently absent from diagonal noise models, and quantifying the severe consequences of this omission is a primary motivation for the present work. A different strand of stochastic ecological modeling introduces noise through environmental (extrinsic) stochasticity rather than demographic (intrinsic) stochasticity [35–38]. Environmental noise models capture variation in parameters (e.g., carrying capacity or predation rate) driven by external fluctuations, and they also typically produce diagonal or block-diagonal noise structures. This paper focuses exclusively on demographic noise, arising from the discreteness of individual events in finite populations, and does not incorporate environmental noise. This distinction is standard in the theory of stochastic population models [39–42] and is strictly maintained throughout. Several authors have emphasized the importance of the mechanistic derivation of noise structure in population models. Renshaw [43] and Allen [44] discuss CTMC formulations of birth–death–immigration processes and their diffusion limits. Naghnaeian and Del Vecchio [45] and Wang [46] develop moment approximations for interacting population models starting from the CME. Recently, Black and McKane [9] provide an accessible account of CME-based stochastic modeling in ecology, while Rogers et al. [47] use CME methods to study demographic fluctuations in multi-species systems. In these works, the importance of deriving the diffusion covariance from event stoichiometry rather than imposing it phenomenologically is recognized. However, the specific structural consequences for predator–prey covariance geometry and their implications for boundary behavior have not been systematically unraveled. Quasi-cycles, stochastic amplification, and resonance near Hopf bifurcations. The phenomenon of quasi-cycles in predator–prey and other ecological systems, which are stochastic oscillations in the stable regime driven by the amplification of demographic fluctuations, was clarified theoretically by McKane and Newman [8]. They used an individual-level model for a predator–prey system to demonstrate this behavior. Their analysis showed that demographic noise can produce power spectral peaks at a characteristic frequency, even when the deterministic equilibrium is stable. This result provides a mechanistic explanation for observed population oscillations without invoking nonlinear limit cycles. This result was extended to R–M-type systems by Barraquand [25], to spatial models by Butler and Goldenfeld [48,49], and to systems with multiple interacting species by Biancalani et al. [50]. Theoretical tools designed to capture these phenomena, such as the linear noise approximation (LNA), rely explicitly on the evaluation of the diffusion covariance matrix at equilibrium, , to quantify these spectral peaks. However, without a bottom-up derivation, D* is often misspecified as diagonal. On the other hand, stochastic sensitivity analysis for nonlinear stochastic systems has been developed by Bashkirtseva and colleagues [51–53], including applications to mathematical ecology models. Their stochastic sensitivity function (SSF) formalism provides a geometric description of fluctuation amplification via Lyapunov-equation-based covariance analysis, and quantifies the local propensity for noise-induced transitions. Consequently, noisy precursors and early-warning signals near tipping points and bifurcations have attracted substantial interest as indicators of regime shifts in ecological systems [54–59]. Our rigorous identification of the negative cross-covariance provides the exact geometric input needed for these sensitivity and early-warning frameworks. This examines whether risk assessments of noise-induced extinction is fundamentally biased by diagonal-noise assumptions. CME and Langevin methods in population biology. The CME and the associated CLE derived via van Kampen system-size expansions are classical [60–62]. Gillespie [62] gives a rigorous derivation of the CLE from the CME under a “many firings” condition, yielding a diffusion approximation whose covariance is the density-rescaled version of the exact CTMC infinitesimal covariance. The key formula appears in this context as a direct mathematical consequence of the event stoichiometry and propensities; see also Kampen [60] and Elf and Ehrenberg [63]. Applications of the CME/CLE framework to ecological population models have been developed by Allen [44], Black and McKane [9], Ovaskainen and Meerson [40], and others. In predator–prey contexts specifically, Pineda-Krch et al. [64] apply Gillespie SSA to predator–prey systems; Alonso et al. [65] study demographic stochasticity and critical transitions using CME-based approximations. While the present paper builds upon this robust CLE/CME framework, it pivots from general moment closures to focus specifically on: (i) the mechanistic isolation of the non-diagonal covariance structure, (ii) the explicit comparison between coupled and split-channel CTMC closures, and (iii) the rigorous treatment of boundary behavior in the resulting diffusion. Well-posedness and boundary analysis for population diffusions. The well-posedness of SDEs arising in mathematical biology has been studied extensively. This is especially true for models with boundary-degenerate or state-constrained diffusion. For one-dimensional diffusions, Feller’s boundary classification [19,66] provides explicit conditions for the accessibility and attainability of boundaries in terms of drift–diffusion balance. Applications to population genetics (Wright–Fisher diffusion) are classical [67–69]. For square-root diffusions (CIR processes) [70], the Feller condition governs whether the boundary is accessible, and is widely used in financial mathematics and population modeling. In two or higher dimensions, direct analogues of the Feller criterion are generally unavailable because boundary attainability depends on the full nonlinear drift–diffusion interaction and may differ across distinct boundary faces. Lyapunov function methods provide the standard replacement: non-explosion and positivity are established by finding a Lyapunov function satisfying a linear growth bound of the generator [71–73]. These methods have been used for stochastic Lotka–Volterra systems and related models [74–77], typically in the context of global well-posedness and long-time behavior under various noise conditions. A relevant model is the mathematical Leslie–Gower system, where a stochastic extension driven by demographic noise is established [78], but the well-posedness analysis is missing. Our work adopts this Lyapunov perspective for the boundary-degenerate demographic diffusion arising from the R–M CTMC, but introduces a critical structural departure from the existing literature: a deliberate, two-stage separation of well-posedness claims. We first establish a maximal strong solution with no interior blow-up before boundary hitting (Theorem 0.6), and only then establish positivity invariance under an additional barrier condition (Proposition 0.7). This separation is explicitly motivated by the modeling bifurcation between open-domain and absorbed formulations. It also prevents positivity invariance, which is biologically relevant only in the open-domain formulation, from being inadvertently embedded as a baseline mathematical assumption in extinction-permitting modeling. Harvesting and interaction. Recent studies have further clarified how harvesting, heterogeneity, and transport mechanisms shape the fate of interacting populations. In a heterogeneous environment with harvesting, the competition of two species was shown to admit a well-posed PDE formulation, and the analysis revealed regimes in which coexistence, competitive exclusion, or even simultaneous extinction may occur, with numerical experiments supporting the theoretical conclusions [79]. A closely related work on two competing species in a nonhomogeneous habitat with harvesting established existence and uniqueness of solutions, identified conditions for coexistence and extinction, and compared the effects of different harvesting levels through both analysis and simulation [80]. More recently, a time-dependent advection-reaction-diffusion competition model with stocking and harvesting was analyzed for existence, uniqueness, and positivity, and the authors further developed fully discrete decoupled linearized schemes and showed numerically that appropriate stocking or harvesting can promote coexistence [81]. In a heterogeneous advective setting, another study investigated a two-species reaction–diffusion–advection model with unequal diffusion and advection rates, deriving local and global stability results and showing that the outcome of competition can depend on the ratio of transport rates [82]. Finally, a spatially distributed harvesting model with resource-based diffusion demonstrated that harvesting intensity and migration strategy can determine whether coexistence or competitive exclusion prevails, and it also provided estimates for harvesting regimes under which coexistence is possible [83]. Together, these works highlight that harvesting and transport effects can fundamentally alter persistence outcomes in structured ecological systems, thereby motivating mechanistically faithful stochastic formulations that track not only mean-field drift but also the associated covariance and boundary behavior. Recent literature has continued to develop stochastic predator–prey theory from both general and mechanism-specific perspectives. A comprehensive treatment of deterministic, stochastic, and thermodynamic modeling of interacting species emphasized how birth and death, competition, consumption, and environmental fluctuations shape the stability and dynamics of ecological systems, thereby providing a broad foundation for stochastic ecological analysis [84]. Within predator–prey systems, one study incorporating additional food for the predator perturbed the prey growth rate and the predator death rate by Gaussian white noise, and established existence and uniqueness of a global positive solution together with extinction and persistence criteria [85]. Related work on a stochastic predator–prey model with additional food and a fear effect proved existence, uniqueness, boundedness, and uniform continuity of the global positive solution, while also deriving conditions for extinction and persistence under environmental noise [86]. Another study combined prey refuge, additional food for the predator, mutual interference, and white-noise perturbations, again proving positivity, boundedness, existence and uniqueness of global positive solutions, and identifying extinction and persistence regimes [87]. These works show that stochastic perturbations, additional food, fear, refuge, and environmental heterogeneity can substantially alter persistence and extinction outcomes, which motivates our mechanistic treatment of covariance structure and boundary-aware stochastic dynamics. From the CME to the diffusion approximation This appendix records a standard but explicit route from a mechanistic CTMC, via the CME, to the diffusion approximation used in the main text. It also clarifies the origin of the covariance identity (44)which underlies the full-covariance chemical Langevin formulation. At the count level, the corresponding infinitesimal covariance is (45)For the predator–prey application in Mechanistic CTMC and CME subsection in Methods section, the key message is that the diffusion matrix is not an ad hoc modeling choice: it is inherited directly from the stoichiometry and propensities of the underlying CTMC/CME description. CTMC formulation and the CME. Let denote the vector of population counts (in the main text, d = 2 for prey and predator). Consider K reaction/event channels indexed by , with stoichiometric increments . Collect them in the stoichiometric matrix When channel k fires at state X, the process jumps according to Let be the propensity of channel k at state X. By definition, for a short interval [t,t + dt), (46)and, up to o(dt), at most one event occurs. Let . Conditioning on what occurs during [t,t + dt), there are (up to o(dt)) two contributions to p(X,t + dt): - (i) no event occurs and the system was already in state X, contributing(47) - (ii) exactly one event k occurs and moves the system from into X, contributing(48) Summing these terms, subtracting p(X,t), dividing by dt, and letting yields the CME: (49)Generator and moment identities. Let be a test function of suitable growth. Multiplying (49) by , summing over states, and using index shifts gives (50)where the CTMC generator is (51)We compute the first moment (drift). Taking gives hence (52)In vector form, this is the exact CTMC drift identity (53)where . We compute the second moment and quadratic variation structure. Taking and using we obtain (54)The term is the key source of infinitesimal covariance; in matrix form it appears as . Infinitesimal covariance at the CTMC level. Let . Conditioning on X(t)=X, equation (46) implies and Therefore, (55)since the outer product of the conditional mean is O(dt2). Hence, the exact infinitesimal covariance per unit time at the count scale is (56)This identity is exact for the CTMC and is the source of the diffusion covariance used in the chemical Langevin approximation. Density-dependent scaling and the diffusion matrix. Introduce a system-size parameter and define density variables Assume density-dependent propensities of Kurtz type: (57)for smooth rate functions . Let . Under the standard chemical Langevin scaling, the density process is approximated by a diffusion with drift and covariance (58)The same formula follows directly by rescaling the CTMC short-time covariance (55): since and , (59)Thus, the diffusion matrix a(x) is precisely the density-scale counterpart of the exact count-scale infinitesimal covariance G(X). Proofs for well-posedness This appendix provides detailed proofs of the main well-posedness statements in Well-posedness subsection in Methods section. We work on the open state space and consider the Itô SDE (60)where W is an r-dimensional Brownian motion and . For convenience, we restate the generator acting on : (61)We also recall the first boundary-hitting time (62)7.3.1. Preliminaries: exhaustion, localization, and a geometric lemma. We begin with standard localization objects and a geometric fact used in the proof of Theorem 0.6. We detail exhaustion of D by precompact sets. Fix an increasing exhaustion of D by compact sets with nonempty interior such that (63)A concrete choice is For a continuous adapted process x, define the exit times (64)Then is increasing in n, and we set (65)Lemma 0.8 (Leaving every compact subset of D while norm-bounded implies boundary approach). Let . Suppose satisfies Furthermore, for every compact , only finitely many xn lie in K. Then Equivalently, there exists a subsequence with after passing to a further subsequence. Proof. Since (xn) is norm-bounded in , Bolzano–Weierstrass theorem yields a convergent subsequence . Suppose, for contradiction, that . Then there exists such that . By convergence, for sufficiently large j we have , so the infinite tail is contained in the compact set . This contradicts the standing assumption that for every compact only finitely many terms of belong to K. Hence , which implies . □ 7.3.2. Complete proof of Theorem 0.6. Proof. The proof is divided into five steps. Step 1: Local strong solutions on each localization domain. Fix n ≥ 1. By assumption (A1), the coefficients are Lipschitz and bounded on . Extend to globally Lipschitz, bounded coefficients on such that (For example, use a cutoff and McShane-type extension componentwise.) Then the globally defined SDE (66)admits a unique global strong solution by standard Itô SDE theory. Define the exit time from Kn for this solution: On , , so b(n) = b and . Therefore X(n) solves the original SDE (60) up to . Step 2: Consistency and pathwise uniqueness patching. Let m > n. Consider the extended solutions X(n) and X(m) driven by the same Brownian motion and the same initial condition. On Kn, both sets of extended coefficients agree with the original coefficients , and (A1) yields local Lipschitz continuity on Kn. By pathwise uniqueness for locally Lipschitz SDEs up to the exit time from Kn, we have Hence the stopped processes are consistent across n, and we may define a process xt on by patching: where is the (common) exit time from Kn under the patched process. This yields a strong solution of (60) on , with By the same local uniqueness argument, this maximal solution is pathwise unique up to . Step 3: A Lyapunov estimate up to bounded-radius localization. Fix R > 0 and define (67)For each R, the stopped process remains in the bounded set , so Itô’s formula applies to . Using (61) and assumption (A2), (68)where is a martingale with zero expectation. Taking expectations gives (69)By Grönwall’s inequality, (70)The right-hand side is independent of R. Step 4: No explosion to infinity before boundary hitting. We show that the process cannot lose well-posedness in the interior by escaping to infinity before hitting . Suppose, for contradiction, that there exists T > 0 such that On the event inside the probability, for all sufficiently large R, and . Since , the process remains in D up to , so within D. By assumption (A2), Thus on a set of positive probability. By Fatou’s lemma, , contradicting the uniform estimate (70). Therefore, (71)Step 5: If every compact set is eventually left while the norm remains bounded, the path approaches . It remains to rule out the possibility that On this event, since with and increasing to D, the trajectory leaves each compact Kn before time . Hence we may choose times such that Because , the sequence is norm-bounded. Moreover, for any compact there exists N such that ; therefore for all n ≥ N. Thus is eventually outside every compact subset of D. By Lemma 0.8, any norm-bounded sequence in D that eventually leaves every compact subset of D admits a subsequence whose distance to tends to zero. Consequently, some subsequence of approaches . Since xt is continuous on , it follows that that is, the path approaches the boundary as . This contradicts the event , which asserts that no boundary hit or exit occurs before time . Hence any loss of compactness in D while the norm remains bounded must correspond to approach to the boundary rather than interior breakdown. Combining this with (71), we conclude This completes the proof. 7.3.3. Complete proof of Proposition 0.7. Proof. Let (Kn) and be the exhaustion and exit times defined in (63)–(64). By Theorem 0.6, a unique maximal strong solution exists up to , and any finite lifetime obstruction can only occur through boundary hitting/approach. Step 1: Localized Itô estimate for the barrier Lyapunov function. Fix n ≥ 1 and t ≥ 0. Since for all s, Itô’s formula for yields (72)where is a true martingale. Taking expectations, (73)Applying Grönwall, (74)Hence (75)Step 2: Contradiction argument if boundary is hit with positive probability. Assume for contradiction that . Then there exists T > 0 such that On the event , continuity of sample paths implies as from below (or equivalently within D). By the barrier blow-up property (B1), Because and Theorem 0.6 implies a.s., the stopping sequence approaches from below on . In particular, (76)Applying Fatou’s lemma to the nonnegative random variables , we obtain (77)contradicting the uniform bound (75) at time T. Therefore for every T > 0, hence Step 3: Globality and positivity invariance. Since boundary hitting is a.s. impossible and Theorem 0.6 excludes interior breakdown before boundary hitting, the maximal solution extends for all times: Moreover, This proves positivity invariance. 7.3.4. Remarks on verifiability of the assumptions in the present model. Remark 0.9 (Local regularity of b and a on D). For the R–M demographic diffusion constructed from Mechanistic CTMC and CME subsection to CME-consistent covariance structure subsection in Methods section, the drift b(x) and covariance a(x) are polynomial/rational combinations of N,P with denominator 1 + N. Since N > 0 on , the denominator never vanishes on D, and therefore b and a are smooth on D (in particular, locally Lipschitz on every ). Remark 0.10 (Local regularity of common factorizations ). Two coefficient factorizations used in the paper fit assumption (A1) on precompact subsets : - (i) Event-based factorization. Ifthen on , each intensity fj(x) is smooth and bounded away from singularities; in the present model the relevant fj are nonnegative smooth functions on D. Hence is locally Lipschitz on K (in particular when fj is bounded away from 0, and more generally for the explicit polynomial/rational intensities used here on compacts). - (ii) Matrix square-root factorization. If is chosen as a measurable square root of a(x) (e.g., Cholesky where a(x) is positive definite), then one typically works locally on compact subsets on which rank/signature does not change. On such sets, can be chosen locally Lipschitz. For the well-posedness statements in Well-posedness subsection in Methods section, only local Lipschitz regularity on each is required. Remark 0.11 (Why boundary degeneracy does not conflict with the open-domain theorem). Demographic diffusions are often degenerate at because reaction intensities vanish as and/or . This does not contradict assumption (A1), which is imposed only on precompact subsets , i.e., strictly away from the boundary. The boundary behavior is handled separately through the maximal-solution formulation and, when needed, the barrier-Lyapunov condition of Proposition 0.7. 7.3.5. Abstract lemma: Two-stage Lyapunov architecture on open domains. The arguments above can be summarized in a reusable form for boundary-degenerate ecological diffusions on open domains. Lemma 0.12 (Two-stage Lyapunov template on an open domain). Let be open, and consider with coefficients locally Lipschitz and locally bounded on every precompact . Assume: - (a) there exists with as in D, and - (b) optionally, there exists with as within D, and (locally on an exhaustion) Then: - (i) under (a), there exists a unique maximal strong solution, and no loss of well-posedness occurs in the interior before boundary hitting; - (ii) under (a)+(b), the boundary is a.s. not hit in finite time, hence the solution is global and remains in D for all t ≥ 0. Proof. Part (i) follows the localization/patching and V∞-Lyapunov estimate used in the proof of Theorem 0.6; part (ii) follows by the barrier contradiction argument used in the proof of Proposition 0.7. Lemma 0.12 is not needed for any statement in the main text, but it clarifies the general structure of the two-stage Lyapunov architecture and may be useful for extensions to other demographic diffusions and food-web models. References - 1. Rosenzweig ML, MacArthur RH. Graphical representation and stability conditions of predator-prey interactions. The American Naturalist. 1963;97(895):209–23. - 2. Rosenzweig ML. Paradox of enrichment: destabilization of exploitation ecosystems in ecological time. Science. 1971;171(3969):385–7. pmid:5538935 - 3. May RM. Stability and complexity in model ecosystems. 1st ed. Princeton: Princeton University Press. 2001. - 4. Beay LK, Suryanto A, Darti I, Trisilowati. Stability of a stage-structure Rosenzweig-MacArthur model incorporating Holling type-II functional response. IOP Conf Ser: Mater Sci Eng. 2019;546(5):052017. - 5. Gilpin ME, Rosenzweig ML. Enriched predator-prey systems: theoretical stability. Science. 1972;177(4052):902–4. pmid:17780992 - 6. Arditi R, Ginzburg LR. Coupling in predator-prey dynamics: Ratio-Dependence. Journal of Theoretical Biology. 1989;139(3):311–26. - 7. Roy S, Chattopadhyay J. The stability of ecosystems: a brief overview of the paradox of enrichment. J Biosci. 2007;32(2):421–8. pmid:17435332 - 8. McKane AJ, Newman TJ. Predator-prey cycles from resonant amplification of demographic stochasticity. Phys Rev Lett. 2005;94(21):218102. pmid:16090353 - 9. Black AJ, McKane AJ. Stochastic formulation of ecological models and their applications. Trends Ecol Evol. 2012;27(6):337–45. pmid:22406194 - 10. Grunert K, Holden H, Jakobsen ER, Stenseth NC. Evolutionarily stable strategies in stable and periodically fluctuating populations: The Rosenzweig-MacArthur predator-prey model. Proc Natl Acad Sci U S A. 2021;118(4):e2017463118. pmid:33479183 - 11. Yu J, Wang LS, Liu Z, Liu J. Pattern suppression and recovery under one-way versus two-way chemotactic coupling in hybrid partial differential equation–ordinary differential equation models. Transport Phenomena. 2026. https://doi.org/10.1515/tp-2026-0023 - 12. Liu Z, Wang LS, Yu J, Zhang J, Martel E, Li S. Bidirectional Endothelial Feedback Drives Turing-Vascular Patterning and Drug-Resistance Niches: A Hybrid PDE-Agent-Based Study. Bioengineering (Basel). 2025;12(10):1097. pmid:41155096 - 13. Turchin P, Taylor AD. Complex Dynamics in Ecological Time Series. Ecology. 1992;73(1):289–305. - 14. Kendall BE, Briggs CJ, Murdoch WW, Turchin P, Ellner SP, McCauley E, et al. Why do populations cycle? a synthesis of statistical and mechanistic modeling approaches. Ecology. 1999;80(6):1789–805. - 15. Barraquand F, Louca S, Abbott KC, Cobbold CA, Cordoleani F, DeAngelis DL, et al. Moving forward in circles: challenges and opportunities in modelling population cycles. Ecol Lett. 2017;20(8):1074–92. pmid:28633194 - 16. Brook BW, O’Grady JJ, Chapman AP, Burgman MA, Akçakaya HR, Frankham R. Predictive accuracy of population viability analysis in conservation biology. Nature. 2000;404(6776):385–7. pmid:10746724 - 17. Henle K, Sarre S, Wiegand K. The Role of Density Regulation in Extinction Processes and Population Viability Analysis. Biodiversity and Conservation. 2004;13(1):9–52. - 18. Steifetten Ø, Dale S. Viability of an endangered population of ortolan buntings: The effect of a skewed operational sex ratio. Biological Conservation. 2006;132(1):88–97. - 19. Feller W. Diffusion processes in one dimension. Transactions of the American Mathematical Society. 1954;77(1):1–31. - 20. Peskir G. On boundary behaviour of one-dimensional diffusions: from brown to feller and beyond. Selected Papers II. Cham: Springer International Publishing. 2015. 77–93. - 21. Liang Y, Wang LS, Yu J, Liu Z. Global Well-Posedness and Stability of Nonlocal Damage-Structured Lineage Model with Feedback and Dedifferentiation. Mathematics. 2025;13(22):3583. - 22. Wang LS, Yu J, Liu Z. A damage-structured PDE model of stem cell hierarchies: The dual role of dedifferentiation in tissue homeostasis and aging. PLoS One. 2026;21(2):e0335163. pmid:41734247 - 23. Gard T. Persistence in stochastic food web models. Bulletin of Mathematical Biology. 1984;46(3):357–70. - 24. Yuan S, Wang Z. Bifurcation and chaotic behavior in stochastic Rosenzweig–MacArthur prey–predator model with non-Gaussian stable Lévy noise. International Journal of Non-Linear Mechanics. 2023;150:104339. - 25. Barraquand F. No sensitivity to functional forms in the Rosenzweig-MacArthur model with strong environmental stochasticity. J Theor Biol. 2023;572:111566. pmid:37422068 - 26. Braumann CA. Variable effort harvesting models in random environments: generalization to density-dependent noise intensities. Math Biosci. 2002;177–178:229–45. https://doi.org/10.1016/s0025-5564(01)00110-9 pmid:11965257 - 27. Braumann CA. Itô versus Stratonovich calculus in random population growth. Math Biosci. 2007;206(1):81–107. pmid:16214183 - 28. Liu M, Wang K, Wu Q. Survival analysis of stochastic competitive models in a polluted environment and stochastic competitive exclusion principle. Bull Math Biol. 2011;73(9):1969–2012. pmid:20725797 - 29. Zhao Y, Yuan S, Ma J. Survival and Stationary Distribution Analysis of a Stochastic Competitive Model of Three Species in a Polluted Environment. Bull Math Biol. 2015;77(7):1285–326. pmid:25990043 - 30. Huang Y, Shi W, Wei C, Zhang S. A stochastic predator-prey model with Holling II increasing function in the predator. J Biol Dyn. 2021;15(1):1–18. pmid:33357105 - 31. Du NH, Nguyen DH, Yin GG. Conditions for permanence and ergodicity of certain stochastic predator–prey models. J Appl Probab. 2016;53(1):187–202. - 32. Qi H, Meng X. Threshold behavior of a stochastic predator–prey system with prey refuge and fear effect. Applied Mathematics Letters. 2021;113:106846. - 33. Yasin MW, Ahmed N, Iqbal MS, Raza A, Rafiq M, Eldin EMT, et al. Spatio-temporal numerical modeling of stochastic predator-prey model. Sci Rep. 2023;13(1):1990. pmid:36737648 - 34. Liu Q, Jiang D, Hayat T, Alsaedi A. Dynamics of a Stochastic Predator–Prey Model with Stage Structure for Predator and Holling Type II Functional Response. J Nonlinear Sci. 2018;28(3):1151–87. - 35. Wu DM, Wang H, Yuan SL. Stochastic sensitivity analysis of noise-induced transitions in a predator-prey model with environmental toxins. Math Biosci Eng. 2019;16(4):2141–53. pmid:31137203 - 36. Chatterjee A, Abbasi MA, Venturino E, Zhen J, Haque M. A predator–prey model with prey refuge: under a stochastic and deterministic environment. Nonlinear Dyn. 2024;112(15):13667–93. - 37. Das A, Samanta GP. Modelling the effect of resource subsidy on a two-species predator-prey system under the influence of environmental noises. Int J Dynam Control. 2021;9(4):1800–17. - 38. Belabbas M, Ouahab A, Souna F. Rich dynamics in a stochastic predator-prey model with protection zone for the prey and multiplicative noise applied on both species. Nonlinear Dyn. 2021;106(3):2761–80. - 39. Newman TJ, Ferdy J-B, Quince C. Extinction times and moment closure in the stochastic logistic process. Theor Popul Biol. 2004;65(2):115–26. pmid:14766186 - 40. Ovaskainen O, Meerson B. Stochastic models of population extinction. Trends Ecol Evol. 2010;25(11):643–52. pmid:20810188 - 41. Wang LS, Yu J, Li S, Liu Z. Analysis and Mean-Field Limit of a Hybrid PDE-ABM Modeling Angiogenesis-Regulated Resistance Evolution. Mathematics. 2025;13(17):2898. - 42. Dobramysl U, Täuber UC. Environmental versus demographic variability in stochastic predator–prey models. J Stat Mech. 2013;2013(10):P10001. - 43. Renshaw E. Modelling biological populations in space and time. 1st ed. Cambridge: Cambridge Univ. Press. 1995. - 44. Allen LJS. An introduction to stochastic processes with applications to biology. 2nd ed. Boca Raton, FL: Chapman and Hall/CRC. 2010. https://doi.org/10.1201/b12537 - 45. Naghnaeian M, Del Vecchio D. Robust moment closure method for the chemical master equation. In: 2017 IEEE Conference on Control Technology and Applications (CCTA), 2017. 967–72. https://doi.org/10.1109/ccta.2017.8062585 - 46. Wang LS, Yu J. Analysis framework for stochastic predator–prey model with demographic noise. Results in Applied Mathematics. 2025;27:100621. - 47. Rogers T, McKane AJ, Rossberg AG. Demographic noise can lead to the spontaneous formation of species. EPL. 2012;97(4):40008. - 48. Butler T, Goldenfeld N. Robust ecological pattern formation induced by demographic noise. Phys Rev E Stat Nonlin Soft Matter Phys. 2009;80(3 Pt 1):030902. pmid:19905053 - 49. Butler T, Goldenfeld N. Fluctuation-driven Turing patterns. Phys Rev E Stat Nonlin Soft Matter Phys. 2011;84(1 Pt 1):011112. pmid:21867118 - 50. Biancalani T, Rogers T, McKane AJ. Noise-induced metastability in biochemical networks. Phys Rev E Stat Nonlin Soft Matter Phys. 2012;86(1 Pt 1):010106. pmid:23005355 - 51. Bashkirtseva IA, Ryashko LB. Stochastic sensitivity of 3D-cycles. Mathematics and Computers in Simulation. 2004;66(1):55–67. - 52. Bashkirtseva I, Ryashko L. Stochastic sensitivity of the closed invariant curves for discrete-time systems. Physica A: Statistical Mechanics and its Applications. 2014;410:236–43. - 53. Bashkirtseva I, Ryashko L. Stochastic Sensitivity Analysis of Noise-Induced Extinction in the Ricker Model with Delay and Allee Effect. Bull Math Biol. 2018;80(6):1596–614. pmid:29611109 - 54. Scheffer M, Bascompte J, Brock WA, Brovkin V, Carpenter SR, Dakos V, et al. Early-warning signals for critical transitions. Nature. 2009;461(7260):53–9. pmid:19727193 - 55. Dakos V, Scheffer M, van Nes EH, Brovkin V, Petoukhov V, Held H. Slowing down as an early warning signal for abrupt climate change. Proc Natl Acad Sci U S A. 2008;105(38):14308–12. pmid:18787119 - 56. Boettiger C, Hastings A. Quantifying limits to detection of early warning for critical transitions. J R Soc Interface. 2012;9(75):2527–39. pmid:22593100 - 57. Carpenter SR, Cole JJ, Pace ML, Batt R, Brock WA, Cline T, et al. Early warnings of regime shifts: a whole-ecosystem experiment. Science. 2011;332(6033):1079–82. pmid:21527677 - 58. Lenton TM. Early warning of climate tipping points. Nature Clim Change. 2011;1(4):201–9. - 59. Boettiger C, Hastings A. No early warning signals for stochastic transitions: insights from large deviation theory. Proc Biol Sci. 2013;280(1766):20131372. pmid:23843397 - 60. van Kampen NG. Stochastic processes in physics and chemistry. 3rd ed. Amsterdam ; Boston: Elsevier. 2007. - 61. Gardiner CW. Stochastic methods: a handbook for the natural and social sciences. 4th ed. Berlin Heidelberg: Springer. 2009. - 62. Gillespie DT. The chemical Langevin equation. The Journal of Chemical Physics. 2000;113(1):297–306. - 63. Elf J, Ehrenberg M. Fast evaluation of fluctuations in biochemical networks with the linear noise approximation. Genome Res. 2003;13(11):2475–84. pmid:14597656 - 64. Pineda-Krch M. GillespieSSA: A user-friendly stochastic simulation package for R. Nat Prec. 2009. https://doi.org/10.1038/npre.2009.3673.1 - 65. Alonso D, McKane AJ, Pascual M. Stochastic amplification in epidemics. J R Soc Interface. 2007;4(14):575–82. pmid:17251128 - 66. Feller W. The Parabolic Differential Equations and the Associated Semi-Groups of Transformations. The Annals of Mathematics. 1952;55(3):468. - 67. Ethier SN, Kurtz TG. Markov processes: characterization and convergence. 1st ed. New York: Wiley. 1986. https://doi.org/10.1002/9780470316658 - 68. Hofrichter J. On the diffusion approximation of Wright–Fisher models with several alleles and loci and its geometry. 2014. - 69. Ewens WJ. Mathematical population genetics. New York, NY: Springer New York. 2004. https://doi.org/10.1007/978-0-387-21822-9 - 70. Cox JC, Ingersoll JE, Ross SA. A Theory of the Term Structure of Interest Rates. Econometrica. 1985;53(2):385. - 71. Khasminskii R. Stochastic stability of differential equations. vol. 66 of Stochastic Modelling and Applied Probability. Berlin, Heidelberg: Springer Berlin Heidelberg. 2012. https://doi.org/10.1007/978-3-642-23280-0 - 72. Meyn SP, Tweedie RL. Stability of Markovian processes III: Foster–Lyapunov criteria for continuous-time processes. Advances in Applied Probability. 1993;25(3):518–48. - 73. Mao X. Stochastic differential equations and applications. 2nd ed. Chichester: Horwood Pub. 2008. - 74. Cui Y, Huang S, Li X, Shi S. The long-term dynamic behaviors of a stochastic predator-prey model with fear and group defense. DCDS-B. 2026;34(0):106–40. - 75. Zhao X, Liu L, Liu M, Fan M. Stochastic dynamics of coral reef system with stage-structure for crown-of-thorns starfish. Chaos, Solitons & Fractals. 2024;181:114629. - 76. Zhan J, Wei Y. Dynamical behavior of a stochastic non-autonomous distributed delay heroin epidemic model with regime-switching. Chaos, Solitons & Fractals. 2024;184:115024. - 77. Liu G, Qi H, Chang Z, Meng X. Asymptotic stability of a stochastic May mutualism system. Computers & Mathematics with Applications. 2020;79(3):735–45. - 78. Wang LS, Yu J. Algebraic–spectral thresholds and discrete–continuous stability transfer in Leslie–Gower systems. era. 2026;34(1):251–90. - 79. Adan MMIY, Kamrujjaman M, Molla MM, Mohebujjaman M, Buenrostro C. Interplay of harvesting and the growth rate for spatially diversified populations and the testing of a decoupled scheme. Math Biosci Eng. 2023;20(4):6374–99. pmid:37161112 - 80. Yasin Adan MMI. Evolution of harvesting in competitive population dynamics with spatial heterogeneity. GANIT: J Bangladesh Math Soc. 2025;45(1). - 81. Tisha MS, Kamrujjaman Md, Mohebujjaman M, Khan T. Decoupled algorithms and analyses for an advection-reaction-diffusion model with stocking and harvesting. Computers & Mathematics with Applications. 2025;189:24–47. - 82. Tisha MS, Kamrujjaman Md, Zahan I. Dynamics of reaction–diffusion–advection system and its impact on river ecology in the presence of spatial heterogeneity. Partial Differential Equations in Applied Mathematics. 2024;11:100802. - 83. Zahan I, Kamrujjaman M. Evolution of dispersal and the analysis of a resource flourished population model with harvesting. Heliyon. 2024;10(10):e30737. pmid:38770280 - 84. Samanta G. Deterministic, stochastic and thermodynamic modelling of some interacting species. Springer Singapore. 2021. https://doi.org/10.1007/978-981-16-6312-3 - 85. Das A, Samanta GP. Stochastic prey–predator model with additional food for predator. Physica A: Statistical Mechanics and its Applications. 2018;512:121–41. - 86. Das A, Samanta GP. Modeling the fear effect on a stochastic prey–predator system with additional food for the predator. J Phys A: Math Theor. 2018;51(46):465601. - 87. Das A, Samanta GP. A prey–predator model with refuge for prey and additional food for predator in a fluctuating environment. Physica A: Statistical Mechanics and its Applications. 2020;538:122844.