{"id":"369d6285-d455-4c77-976e-0051d1859089","arxiv_id":"2508.10909","paper_version":1,"verdict":"CONDITIONAL","confidence":"HIGH","novelty_score":7.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":4,"one_line_summary":"A phantom-reactant sampling scheme lets particle-based simulations reproduce non-elementary bimolecular kinetics (Michaelis-Menten and Hill) without explicitly simulating fast elementary reactions.","lead":"This paper introduces a way to simulate biochemical reactions with non-elementary (Michaelis-Menten and Hill) kinetics directly at the particle level, using an imaginary 'phantom' molecule to provide a second spatial coordinate. The method is validated against Michaelis-Menten theory and applied to a circadian clock model, enabling spatial stochastic simulation of reduced reaction networks without resolving fast enzyme-binding steps.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Central rate derivation assumes a Poisson nearest-neighbor distribution for S (Eq 28), but the paper never tests Assumption 1 at low copy number or with clustered S; the claimed accuracy is therefore conditional on an unvalidated well-mixed regime.","rationale":"The paper's contribution is a way to port non-elementary rate laws into particle simulations by introducing a phantom reactant. The mathematical construction is coherent when the background of the 'second reactant' S is spatially uniform at the scale probed by the reaction boundary, because that is the only case in which Eq 7 can be inverted to obtain F(r2). The weakest point in the argument is exactly this: the method's local reaction kernel F(r2) is not derived from first principles for arbitrary S configurations; it is calibrated to the mean-field Poisson density. The paper acknowledges this as Assumption 1 but does not test it or quantify its breakdown. The validations in Section 5.1 deliberately stay in the safe regime (NS≥48, uniform, non-consumed S), and the Goldbeter demonstration does not include an explicit elementary-particle baseline, so it cannot distinguish genuine spatial biology from errors induced by the reaction condition near boundaries. A focused numerical experiment with clustered or low-copy S would settle whether the central claim holds outside the well-mixed regime. I agree with the reader's weakest_assumption; the concern does not overturn the paper's internal logic, but it does mean the central claim is conditional on a condition that is plausible yet unvalidated. Therefore the verdict should remain CONDITIONAL, and the recommended adjustment is UNCHANGED. The computational cost claim is also unsupported by runtime measurements, but the correctness concern above is more load-bearing.","tokens_in":27750,"tokens_out":23279,"duration_ms":278235,"concrete_test":"Re-run the Section 5.1 Michaelis-Menten validation with the same dimensionless parameters and 10^5 replicates, but with two perturbations: (a) initial S molecules uniformly distributed in one half of the cubic domain, leaving the other half empty; and (b) NS ∈ {1, 2, 5, 10} with otherwise uniform placement. Compute the mean steady-state NE and the implied reaction rate, and compare the signed relative error to the standard error of the mean, as in Figures 5c and 5d. If the error exceeds the error bars for any of these cases, Assumption 1 is a genuine load-bearing limitation and the paper should state the method's domain of applicability; if the error remains within error bars, the concern is resolved and the conditional verdict can be relaxed.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central claim requires Algorithm 1 to reproduce a prescribed non-elementary rate K(s). Equation (7) derives the reaction boundary from the homogeneous Poisson nearest-neighbor density P(r2)=4πr2^2 s exp(-4πs r2^3/3), and Section 6 makes this an explicit condition (Assumption 1). Algorithm 1 evaluates F(r2) using the actual distance to the closest S, so if S is not well mixed on the scale s^{-1/3} (the 99th percentile of Eq 28), the simulated rate will deviate from K(s). This is not a hypothetical edge case: the Section 5.1 validation deliberately uses NS≥48 with uniform random placements and a fixed substrate concentration (S is not consumed), so it never exercises the assumption. The Goldbeter spatial simulation (Section 5.2) imposes reflecting boundaries and a membrane that locally break the translational invariance of Eq 28; the observed discrepancy between spatial and calibration simulations is attributed to 'spatial effects,' but the reaction probability itself was derived assuming a homogeneous S density, so part of that discrepancy may be an artifact of applying F(r2) outside its validity regime. The paper provides no error bound or sensitivity analysis for these cases, leaving the central claim valid only under an unquantified well-mixed condition.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper presents an event-driven particle-based simulation framework for directly simulating non-elementary bimolecular kinetics. The key idea is to introduce a conceptual 'phantom' reactant that mimics the third spatial degree of freedom in previous trimolecular reaction conditions (Kearney and Flegg, J. Chem. Phys. 161, 194111 (2024)). The reaction probability F(r2) is obtained by inverting Eq. (7), which relates the target rate K(s) to a proximity-based reaction boundary. The framework is validated on a Michaelis-Menten system (Section 5.1) and demonstrated on a modified Goldbeter circadian model (Section 5.2), with a calibration simulation matching SSA and a spatial simulation showing expected boundary-induced deviations. The authors state three explicit assumptions, the most important being Assumption 1: the second reactant S must be well-mixed on the scale s^{-1/3} so that the nearest-neighbor density in Eq. (28) is accurate.","tokens_in":28102,"tokens_out":12721,"duration_ms":140383,"significance":"If the framework is correct, it broadens particle-based reaction-diffusion simulation to reduced non-elementary kinetics without resolving fast underlying elementary reactions, which could bring practical efficiency gains for enzymatic and signaling systems. The paper is clearly written and provides a concrete algorithm, explicit assumptions, and a non-trivial demonstration. It also correctly identifies the RDME's failure for non-elementary propensities, citing Smith and Grima. However, the central accuracy claim is only established in the well-mixed regime, and the Goldbeter demonstration uses adjusted parameters (n=2, revised vs, vd, KI, kN, kC) rather than the original model. The MM validation is partly a consistency check because F(r2) is derived from the target rate by construction. These caveats limit the strength of the conclusions as currently stated.","major_comments":[{"comment":"The central claim that the method 'accurately reproduces the target non-elementary kinetics' (abstract, Section 7) is only tested in the well-mixed regime of Assumption 1. The MM validation uses NS≥48 with uniform random initial placement and with S never consumed, so it never exercises the clustered or low-copy-number cases where Eq. (28) fails. The Goldbeter spatial simulation (Section 5.2) breaks translational invariance via reflecting boundaries and the membrane, and the observed deviation from the calibration simulation is attributed to 'spatial effects' without confirming that Eq. (28) remains valid near those boundaries. Since the method is proposed for spatial particle-based simulation, the behavior when Assumption 1 is violated is load-bearing. Please add a sensitivity analysis for clustered S or low NS, or explicitly restrict the abstract and conclusion claims to the well-mixed condition.","section":"Section 5.1, Discussion Assumption 1"},{"comment":"The text refers to the 'classical Goldbeter model' in the abstract and Section 5.2, but the simulations use n=2 instead of the original n=4 and altered values for vs, vd, KI, kN, and kC (as acknowledged in the text). This is a modified Goldbeter-type model, not the original. Please label it as a modified model throughout. More importantly, the spatial simulation is not validated against any ground truth, such as a particle simulation with explicit enzymes or a spatially resolved stochastic simulation; therefore the observed discrepancy cannot be cleanly separated into genuine spatial effects and possible artifacts of applying F(r2) outside its validity regime. Adding such a comparison would substantially strengthen the demonstration.","section":"Section 5.2, Table 1, Abstract"},{"comment":"The reaction probability F(r2) is derived by inverting Eq. (7) from the target MM rate, so the MM validation in Section 5.1 is essentially a consistency check: it confirms that the algorithm implements the intended rate under well-mixed conditions, but it does not independently validate the spatial accuracy of the boundary. The independent grounding comes from the Goldbeter calibration simulation agreeing with SSA, which is stronger. Still, the paper should state this circularity explicitly and, if feasible, test the boundary for a target rate with a different functional form (e.g., a single Hill-type reaction with n=2) to show the inversion works beyond the case it was derived from.","section":"Section 3, Eq. (7), Section 5.1"}],"minor_comments":[{"comment":"Eq. (7) contains the factor 'V sexp(...)' which is ambiguous. It should be written with explicit parentheses, e.g., '4π D1 f(r2) V s exp(...) 4π r2^2 dr2' or with the intended grouping clarified, to avoid confusion about whether V and s multiply the exponential.","section":"Eq. (7)"},{"comment":"The sentence 'we have been careful to ensure that the domain contains enough molecules of S for the probability density assumed by Eq. (7) to be accurate' is vague. Please specify the criterion (e.g., the minimum NS used and why that is sufficient) and cite the relevant analysis from Kearney et al. [78].","section":"Section 5.1"},{"comment":"Algorithm 1 computes r2 as the distance from the reactive molecule to the closest S (line 13), but the text says r2 is the distance from the centre of diffusion of the reactive molecule and the phantom to the closest S. The algorithm as written implicitly assumes D1→∞ so that x̄1≈x0. Please state this explicitly in the algorithm description, or present the general form with the phantom position sampling.","section":"Section 4, Algorithm 1"},{"comment":"The description of the spatial simulation domain is unclear: the text says 'we do not explicitly separate the nucleus from the cytoplasm' but then uses volumes VC and VN in Eq. (27). Please clarify the geometry of the cube, the location of the membrane, and how the effective volumes are handled.","section":"Section 5.2"},{"comment":"In the abstract, 'to biomolecular reactions' should likely be 'to bimolecular reactions' to match the title and the rest of the text; 'biomolecular' refers to biological molecules and is not the intended meaning here.","section":"Abstract"},{"comment":"The discussion of the spatial simulation's discrepancy attributes the effect to reduced reaction probability near boundaries, but it does not mention that the F(r2) formula itself was derived assuming a spatially homogeneous S density. Please add a sentence noting that this is an additional possible source of deviation and that quantifying it would require a boundary-specific analysis.","section":"Section 6, Discussion"}],"recommendation":"major_revision","confidential_remarks":"The paper is technically sound in its core derivation, and the phantom approach is a clever and useful extension of the authors' prior trimolecular work. The main weakness is not an error but an overstatement: the validity of the reaction probability is conditional on a well-mixed assumption that is not tested in the spatial regime where the method is expected to be used. A sensitivity analysis or an explicit scope restriction would address this. The Goldbeter demonstration would be more convincing if the spatial simulation were checked against an explicit-enzyme particle simulation or another spatially resolved reference. I would not reject; the issues are fixable within the manuscript's scope."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"The paper delivers something real: a practical way to put non-elementary bimolecular kinetics into a particle-based simulation without resolving every fast elementary step. The phantom reactant is a genuinely useful trick. By sampling arrival times from the Smoluchowski rate and closest-approach distances from a uniform CDF, the authors recover the second spatial degree of freedom that a trimolecular reaction condition needs, and they package it into an event-driven algorithm that is simple to implement. The Michaelis-Menten validation is clean so far as it goes: over a wide range of substrate concentrations, 1e5 replicates track the expected rate well. The Goldbeter calibration simulation matching SSA is a good sign, and the discussion is honest about the n=4 Hill boundary being unusable and about the open question of negative reaction probabilities.\n\nThe soft spots are real, though. The central derivation leans on Assumption 1: the closest-S distance must be governed by the homogeneous Poisson nearest-neighbor density P(r2) = 4πr2^2 s exp(-4πs r2^3/3). That condition is load-bearing, and the paper never tests it at low copy number or with clustered S. The MM validation deliberately keeps S fixed and uniform with NS>=48, so it does not exercise the assumption. The spatial Goldbeter run uses reflecting boundaries and a membrane that break translational invariance; the authors attribute the period and amplitude shift to spatial effects, but part of that shift could be the reaction probability F(r2) being applied outside the regime where Eq (28) is accurate. Without an explicit elementary-particle baseline or a sensitivity analysis, we cannot tell how much is physics and how much is artifact.\n\nThe MM test is also partly circular: F(r2) is derived by inverting Eq (7) from the target rate, so the agreement is largely a consistency check. The Goldbeter comparison to SSA provides independent grounding, which is why I would not call the circularity fatal. Still, the abstract's claim of significantly reduced computational cost is unsupported: there are no timings, no code, no data. The key interaction radius sigma is never given a numerical value, which makes reproduction harder. The citation pattern looks fine; the authors build on their own prior work but cite alternatives and relevant limitations literature.\n\nThis is a useful contribution to spatial stochastic simulation, not a revolution. It deserves a serious referee. The referee should push for (1) tests of Assumption 1 with low copy numbers and non-uniform S, ideally with an error bound; (2) an explicit-particle simulation of the Goldbeter spatial model as baseline; (3) actual timings; (4) code/data or at least full parameter values. Address those and the paper is publishable as a solid methods paper. I would take it for peer review.","headline":"A genuinely useful phantom-reactant trick for simulating non-elementary bimolecular kinetics, but the central well-mixed assumption is undertested and the cost claim is unmeasured.","tokens_in":28540,"tokens_out":3395,"would_cite":true,"duration_ms":39702,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"This paper introduces a phantom-mediated reaction probability that lets particle-based simulations directly reproduce Michaelis-Menten and circadian non-elementary kinetics without resolving the fast elementary binding steps behind them.","keywords":["particle-based simulation","non-elementary kinetics","Michaelis-Menten kinetics","phantom reactant","reaction boundary","circadian oscillations","stochastic reaction-diffusion","event-driven simulation"],"falsifier":"Run the Michaelis-Menten validation system with substrate molecules placed in tight clusters rather than uniformly, keeping the same global concentration, and measure the steady-state enzyme consumption rate; the paper's Assumption 1 predicts that once the $99$th percentile of nearest-neighbor distance ($s^{-1/3}$) approaches the domain scale or the reaction radius, the measured rate should systematically depart from $K(s)=4\\pi\\hat D_1\\sigma s/(V(\\Gamma+s))$ and should converge back to it as the substrate is re-randomized.","tokens_in":27580,"feed_emoji":"🧬","tokens_out":8468,"duration_ms":92956,"temperature":0.7,"pith_summary":"This paper introduces a particle-based simulation method that directly reproduces non-elementary bimolecular reaction kinetics, such as saturating enzyme kinetics described by the Michaelis-Menten law, without simulating the fast elementary binding steps that produce them. The trick is to mimic a third 'phantom' reactant: reactive events are scheduled according to an effective collision rate, and each event is accepted with a distance-dependent probability derived from the target reaction rate. The authors apply the method to a minimal model of circadian oscillations in the fruit fly, obtaining stochastic dynamics comparable to the chemical-master-equation reference when spatial boundaries are periodic. If correct, the method broadens what spatially resolved particle simulations can handle and offers an efficiency gain whenever quasi-steady-state reductions are appropriate.","feed_headline":"Phantom particle lets enzyme kinetics run directly in simulations","feed_subtitle":"It reproduces Michaelis-Menten and circadian-oscillation kinetics at the particle level.","key_machinery":"The central device is the phantom molecule: a fictitious third reactant that is never tracked but whose behavior is fully sampled. Its arrival times follow the diffusion-limited collision rate $K_R=4\\pi\\hat D_1\\sigma/V$, and its distance of closest approach is drawn from the uniform cumulative distribution $\\Phi_{\\infty}(r_1^{\\rm min})=r_1^{\\rm min}/\\sigma$. This turns any reaction boundary of the form $r_1=f(r_2)$ into a direct acceptance probability $F(r_2)=f(r_2)/\\sigma$, evaluated using the distance $r_2$ from the reactive molecule to the nearest molecule of the second species. The machinery carries the argument because it gives bimolecular systems the missing second spatial coordinate that the authors' earlier trimolecular condition required, without tracking an actual third particle.","core_discovery":"The central claim is that non-elementary bimolecular kinetics can be reproduced at the particle level by replacing the explicit third molecule of a trimolecular reaction with a purely conceptual 'phantom' reactant. The phantom is never tracked; instead its arrival times are sampled from the classical diffusion-limited collision rate, and its distance of closest approach to the reactive molecule is sampled uniformly on $[0,\\sigma]$. This restores a second spatial degree of freedom, allowing a reaction boundary of the form $r_1=f(r_2)$ to be converted into a simple acceptance probability $F(r_2)=f(r_2)/\\sigma$. For Michaelis-Menten kinetics the resulting probability is $F_{\\rm MM}(r_2)=\\exp(-4\\pi\\Gamma r_2^3/3)$, and for Hill-type inhibition with $n=2$ it is $F_{H2}(r_2)=\\sin^2(4\\pi K_I r_2^3/6)$. The paper validates the method against the Michaelis-Menten rate law and against a circadian oscillation model, reporting that the simulated kinetics match the target non-elementary rates without simulating the implied fast elementary reactions.","pith_inferences":["One extension left implicit in the paper is that the same phantom-arrival construction could generate acceptance probabilities for any non-elementary rate whose Laplace transform exists, not just Michaelis-Menten and Hill forms.","A testable biological prediction of the paper's spatial circadian simulation is that physical confinement alone can lengthen the oscillation period and raise mean protein levels, because reflective boundaries reduce reaction rates near the membrane.","If enzymes are spatially clustered rather than well mixed, the paper's own Assumption 1 implies that reduced enzyme kinetics should be recomputed with a position-dependent effective rate; the paper does not explore this case."],"forward_implications":["Michaelis-Menten and Hill-type ($n=2$) bimolecular rates can be simulated directly in an event-driven particle code, with no need to resolve enzyme-substrate binding events.","The method can be combined with existing unimolecular event handling and membrane transmission rules, so a five-species circadian oscillator (mRNA, cytoplasmic PER forms, nuclear PER) can be simulated at single-molecule resolution.","In the fast-diffusion limit the scheme recovers chemical-master-equation behavior, whereas the reaction-diffusion master equation is not reliable for non-elementary propensities in that limit.","Because molecular positions are known exactly at each reactive event, the event-driven implementation avoids the missed-reaction error that finite time steps introduce in the earlier trimolecular scheme."],"supporting_citations":[{"why":"Supplies the trimolecular reaction-boundary framework that this paper adapts to bimolecular reactions.","marker":"[73]"},{"why":"Introduces the diffusive separation coordinates used to quantify relative proximity of multiple diffusing molecules.","marker":"[77]"},{"why":"Derives the leading-order flux formula (Equation 7) and the boundary-correction procedure that justify the approximate reaction rates.","marker":"[78]"},{"why":"Supplies the classical diffusion-limited collision rate used to schedule phantom arrival events.","marker":"[47]"},{"why":"Provides the standard particle-based update and reaction-correction techniques on which Algorithm 1 builds.","marker":"[55]"},{"why":"Supplies the membrane transmission probabilities used for nuclear-cytoplasmic exchange in the spatial circadian simulation.","marker":"[102]"},{"why":"Defines the circadian PER oscillation model that the method reproduces in Section 5.2.","marker":"[76]"},{"why":"Provides the chemical-master-equation comparison used to assess the stochastic circadian simulations.","marker":"[81]"},{"why":"Shows the reaction-diffusion master equation breaks down with non-elementary propensities, motivating the particle-based alternative.","marker":"[74]"}],"fun_headline_variants":["Phantom particle enables direct non-elementary kinetics in simulations","Skipping fast steps, phantom particle simulates complex kinetics","Phantom reactant converts trimolecular trick into bimolecular simulation","Particle-based simulation of non-elementary kinetics without fast steps","Phantom particle mimics third molecule for direct non-elementary kinetics"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The method assumes that the substrate molecules defining $r_2$ are well mixed on length scales comparable to the typical distance to the nearest substrate (about $s^{-1/3}$); if they are clustered or present at very low copy number, the nearest-neighbor probability density used in the derivation is no longer accurate and the simulated rate will drift from the target kinetics.","fun_headline_variants_meta":{"raw":{"variants":["Phantom particle enables direct non-elementary kinetics in simulations","Skipping fast steps, phantom particle simulates complex kinetics","Phantom reactant converts trimolecular trick into bimolecular simulation","Particle-based simulation of non-elementary kinetics without fast steps","Phantom particle mimics third molecule for direct non-elementary kinetics"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.001404,"raw_usage":{"total_tokens":5689,"prompt_tokens":974,"completion_tokens":4715,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":590,"completion_tokens_details":{"reasoning_tokens":4631}},"tokens_in":590,"tokens_out":4715,"duration_ms":34783,"temperature":1.0,"reasoning_tokens":4631,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-06T10:50:56.538636+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run the Michaelis-Menten validation system with substrate molecules placed in tight clusters rather than uniformly, keeping the same global concentration, and measure the steady-state enzyme consumption rate; the paper's Assumption 1 predicts that once the $99$th percentile of nearest-neighbor distance ($s^{-1/3}$) approaches the domain scale or the reaction radius, the measured rate should systematically depart from $K(s)=4\\pi\\hat D_1\\sigma s/(V(\\Gamma+s))$ and should converge back to it as the substrate is re-randomized.","supporting_citations":[{"cited_title":"Enzyme kinetics simulation at the scale of individual parti- cles","cited_arxiv_id":null,"evidence_quote":"Supplies the trimolecular reaction-boundary framework that this paper adapts to bimolecular reactions."},{"cited_title":"Accurate stochastic simulation of nonlinear reactions between closest particles","cited_arxiv_id":"2504.03215","evidence_quote":"Derives the leading-order flux formula (Equation 7) and the boundary-correction procedure that justify the approximate reaction rates."},{"cited_title":"Toward a detailed computational model for the mammalian circadian clock","cited_arxiv_id":null,"evidence_quote":"Supplies the membrane transmission probabilities used for nuclear-cytoplasmic exchange in the spatial circadian simulation."},{"cited_title":"Deterministic Versus Stochastic Models for Circadian Rhythms","cited_arxiv_id":null,"evidence_quote":"Provides the chemical-master-equation comparison used to assess the stochastic circadian simulations."},{"cited_title":"Breakdown of the reaction-diffusion master equation with nonelementary rates","cited_arxiv_id":null,"evidence_quote":"Shows the reaction-diffusion master equation breaks down with non-elementary propensities, motivating the particle-based alternative."}],"review_version":1}