Variance of finite difference methods for reaction networks with non-Lipschitz rate functions

Chaojie Yuan; David F. Anderson

REVIEW 2 major objections 4 minor 28 references

Variance of finite difference methods for reaction networks with non-Lipschitz rate functions

T0 review · 2 major / 4 minor · reviewed 2026-08-14 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗ Review JSON

Pith's one-line read This paper proves that for a broad class of stochastic reaction networks with non-Lipschitz rate functions—including binary mass-action networks and models with time-dependent parameters—the stacked-coupling finite-difference sensitivity…

desk verdict Real result, repairable proof gap: the O(ε) bound holds, but the conditioning step in §3.3 is not justified as written and the binary-systems claim in the abstract is too strong. read the letter →

arxiv 1908.06880 v2 pith:AG6F5QFJ submitted 2019-08-19 math.NA cs.NAmath.PRq-bio.MN

David F. Anderson , Chaojie Yuan This is my paper · ORCID

classification math.NAcs.NAmath.PRq-bio.MN MSC 60J2760J7565C05

keywords reactionnetworksparametricsensitivityanalysisfinitedifferenceestimatorsstackedcouplingnon-Lipschitzintensityfunctionsstochasticmass-actionkineticstime-dependentrateparametersvariancebounds

verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

Sensitivity analysis asks how model outputs change with kinetic parameters, and finite-difference estimators are the workhorse for stochastic reaction networks. Their efficiency hinges on the variance of the difference between a process run at perturbed and unperturbed parameters: if it decays linearly in the perturbation ε, the estimator is practical. Earlier theory delivered that linear decay only under globally Lipschitz rate functions, which excludes most biochemical models. This paper proves the same O(ε) variance scaling for a far wider class—including binary mass-action networks and networks with time-dependent rate parameters—under the stacked coupling, confirming that good variance behavior is not an artifact of the old restrictive assumptions.

What carries the argument

The engine is the stacked coupling, which drives both processes by one space-time Poisson point process and sizes each reaction interval by the maximum of the two rates; as a result the processes use the same reaction channel whenever they jump together, and they decouple only when one rate exceeds the other. The proof then splits the expected L1 distance according to the total number $N$ of coupled jumps, bounds the $n$-th term by $n^r P(N=n)$ times a per-event decoupling probability, and shows Assumption 1(3) makes that probability $O(ε)$. A comparison with a linear pure-birth process supplies an exponential tail bound for how long the process can remain inside large L1 balls, which in turn gives finite moments for $N$. These pieces combine to sum the series and yield the $O(ε)$ bound.

What would settle it

Simulate the intracellular viral kinetics model of Example 2 under the stacked coupling for a sequence ε→0 and a fixed time t; Theorem 2 predicts that $E[||X^{θ+ε}(t)-X^θ(t)||_1]/ε$ stays bounded. If the ratio instead diverges (for example scaling like $ε^{-1/2}$), the theorem's quantitative claim is false.

Watch

Extended reading notes

Core claim

The central result is Theorem 2: when Assumption 1 holds and the processes are coupled with the stacked coupling from the same initial condition, $E[||X^{θ+ε}(t)-X^θ(t)||_1^r] ≤ C_{r,t}||ε||_1$ for every $r ≥ 1$ and sufficiently small $||ε||_1$. Assumption 1 asks that net-molecule-increasing reactions have rate functions growing at most linearly, all rates grow at most polynomially, and the relative change in the rate vector caused by perturbing θ is uniformly bounded by a constant times $||ε||_1$. Since stochastic mass-action kinetics automatically satisfies the last two conditions (Lemma 1), the theorem applies to binary networks and, with time-dependent rates, to models like circadian transcription. Taking $r=2$ and passing through a Lipschitz output $f$ gives Corollary 1: $\mathrm{Var}(f(X^{θ+ε}(t))-f(X^θ(t))) ≤ C_{T,f}||ε||_1$, which is exactly the scaling needed for finite-difference Monte Carlo sensitivity estimates.

Load-bearing premise

The result stands on the assumption that changing a parameter by a small amount changes every reaction rate by at most a small constant times that rate itself, uniformly over all possible states and times; if a rate can be tiny while its perturbed version is not, the per-event decoupling probability is no longer proportional to ε and the O(ε) variance bound can fail.

Editorial extensions

If this is right

Binary mass-action networks with only zeroth- or first-order net-producing reactions now fall under the O(ε) variance guarantee without bounded-state assumptions.
Time-dependent rate parameters are covered, so oscillatory gene-transcription models with circadian forcing admit efficient finite-difference sensitivity estimates.
For any Lipschitz output function, the finite-difference estimator's difference variance scales as $C_{T,f}||ε||_1$, the same rate previously known only for globally Lipschitz systems.
The bounded-state and global-Lipschitz assumptions that excluded most literature models can be dropped for a large class of practically used networks.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

The proof mechanism suggests that any coupling based on a shared driving Poisson process whose decoupling probability per event is bounded by the same relative rate difference should inherit the O(ε) variance scaling; the paper itself only proves this for the stacked coupling.
Because Assumption 1(3) is stated in relative rather than absolute terms, it is a natural candidate hypothesis for sensitivity estimation outside the mass-action setting, including empirical rate laws in cell biology.
A conservative network whose rate function is a power of the state (for example A⇌B with A→B rate $x_A^{1+θ}$) satisfies the non-explosiveness and growth conditions but not the relative-sensitivity bound; testing it numerically would sharpen the exact boundary of the theorem's scope.

Share X Bluesky LinkedIn Reddit HN

Signed reviews

No signed human review yet.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Desk Editor's Note

read the letter

The main theorem is genuine and the paper deserves a serious referee, but the proof as printed has a hole in its central decoupling estimate, and the abstract overstates what binary networks satisfy. Both are fixable; neither sinks the result.

What is new: Theorem 2 extends Anderson's 2012 O(ε) variance guarantee from globally Lipschitz intensities to a broad non-Lipschitz class, including time-dependent rates. That is a real step, and the proof machinery is solid. Lemma 2's exponential tail bound is a clear improvement over the polynomial bounds in earlier work, and the hypoexponential calculation in Lemma 3 is done carefully. Lemma 8's moment bound on the coupled counting process is also clean. The mass-action verification in Lemma 1 is straightforward and correct.

The soft spot is the conditioning in §3.3, after equation (33). The paper conditions on the terminal count N(t)=n, which is a future event, and then writes the per-event decoupling probability as the intensity ratio. That ratio is the unconditional hazard given the history up to the event time; conditioning on the future total count biases it, because whether the pair decouples changes the subsequent rates and hence the distribution of N(t). The stress-test note gave a concrete example (rates decreasing in population) where the bias is real. So the equality in (34) and the product bound that follows are not justified as written.

That said, the gap is repairable and the theorem likely survives: replace the conditional bound by the unconditional P(β≤n) ≤ nC‖ε‖, which follows from Assumption 1(3) by iterated expectation, then use the superexponential tail of N from Lemma 8 to sum. The argument would go through with a bit more work. A referee should ask for this repair.

The second issue is presentation: the abstract and §2.4 claim that all binary systems satisfy Assumption 1, but the body itself notes that 2A→3A has a binary source and violates condition (1). The precise statement in Corollary 2 is correct—only zeroth- and first-order reactions may produce net gain. The broader claim should be rewritten.

Who gets value: people working on sensitivity estimators for stochastic reaction networks, especially those using the stacked coupling. The subfield will want this theorem. It deserves peer review, and a competent referee will catch the conditioning gap. I would expect major revision but not rejection.

Referee Report

2 major / 4 minor

Summary. The paper analyzes the variance of finite-difference estimators for parameter sensitivities of stochastic reaction networks, focusing on the stacked (split) coupling. The main result, Theorem 2, asserts that under Assumption 1—linear growth of rates for population-increasing reactions, polynomial growth for all rates, and a uniform relative parameter-sensitivity bound—the r-th moment of the L1 difference between the θ- and (θ+ε)-processes is O(||ε||₁). Corollary 1 converts this into the variance scaling Var(f(X^{θ+ε}(t))-f(X^θ(t))) = O(ε) for Lipschitz f. The proof uses an exponential tail bound for exit times (Lemma 2), an exact hypoexponential calculation (Lemma 3), a moment bound for the total number of coupled jumps (Lemma 8), and a decoupling probability argument in Section 3.3. The paper also shows that mass-action kinetics with only zeroth- and first-order growth reactions satisfy Assumption 1, covering many binary networks.

Significance. If correct, the result substantially extends Anderson (2012), which required globally Lipschitz intensities, and gives the first rigorous variance bound for the stacked coupling on a broad class of non-Lipschitz models. The assumptions are explicit and checkable, and the comparison-with-birth-process proof of the exponential tail is elegant. However, I found two load-bearing gaps in the proof: an unjustified conditioning on the future terminal count in Section 3.3 and an incorrect Poisson tail bound in Lemma 8. Both are likely fixable without changing the main claims.

major comments (2)

[Section 3.3, Eq. (34) and following] The step that condition (3) of Assumption 1 gives P(β^{θ,ε}=i | N^{θ,ε}_{qK}(t)=n, β^{θ,ε}>i-1) equal to the displayed intensity ratio is not justified. The conditioning event {N(t)=n} is a terminal condition that depends on post-μ_i behavior, and the event β=i changes the post-jump state and therefore the future intensities. Consequently the conditional probability is not in general equal to the ratio evaluated at the pre-jump state. For instance, for the linear birth reaction A→2A with λ^θ(x)=θx, a split event changes the total intensity after the jump, so the distribution of the total count N(t) differs between split and common paths; conditioning on N(t)=n biases the decoupling probability. This estimate is the central load-bearing step of Theorem 2, so the proof as written is incomplete.
[Section 3.2, Lemma 8 proof] The Poisson tail bound used for Z_n(t) is incorrect as stated. For a Poisson random variable with mean μ, the standard Chernoff bound is P(Z≥a) ≤ exp(-a log(a/μ)+a-μ); the manuscript instead uses exp(-a log(a)/μ+a-μ) in the display leading to (32). For a >> μ, the latter exponent is dominated by +a and the bound is vacuous. This is load-bearing because Lemma 8 supplies the finiteness of E[(N^{θ,ε}_{qK}(t))^{r+1}] used in the final step of Theorem 2. The conclusion is likely true with the correct tail bound, but the proof must be fixed.

minor comments (4)

[Abstract and introduction] The claim that binary systems with or without time-dependent rate parameters satisfy the assumptions is not substantiated by the lemmas. Lemma 1 and Corollary 2 cover mass-action kinetics with time-independent rates; for time-dependent rates, additional hypotheses (e.g., uniform boundedness of rate constants and condition (3)) are needed. Please clarify or prove the claim.
[Assumption 1(3)] The ratio in Assumption 1(3) is undefined when the denominator vanishes, for example at states where all rates are zero. Append a convention such as setting the ratio to 0 in that case.
[Definition 1 and Section 2.3] Definition 1 contains the grammatical error 'A reaction networks is a triple', and Section 2.3 contains the typo 'arrises' instead of 'arises'.
[Lemma 8] In Lemma 8, the integer p in condition (2) may depend on the reaction k; taking p as a single order for all k requires a brief justification, for example by taking p = max_k p_k.

Circularity Check

0 steps flagged · score 0.0 of 10

No significant circularity: the variance bound is derived from Assumption 1, not assumed or fitted.

full rationale

The derivation chain is forward: Assumption 1(1)-(2) are used in Lemmas 2 and 8 to control large states and the total coupled jump count, and Assumption 1(3) gives the per-event decoupling probability; the proof of Theorem 2 then sums these ingredients to obtain the O(ε) moment bound. None of the quantities in the conclusion, such as E[||X^{θ+ε}(t)-X^θ(t)||_1^r] or Var(f(X^{θ+ε}(t))-f(X^θ(t))), is used to define Assumption 1(3), and no parameter is fitted to data and then called a prediction. The stacked coupling from [9] and the earlier global-Lipschitz result [2] are cited as constructions and background, not as premises that already contain Theorem 2; the paper proves the needed lemmas, including Lemmas 2, 4, 5, and 8, internally. The possible conditioning issue around N_{qK}^{θ,ε}(t)=n in Section 3.3 is a proof-correctness concern, not a circularity concern, since it does not turn an output into an input. No circular step can be exhibited from the paper's own equations or citations, so the appropriate finding is no significant circularity.

Assumptions & free parameters 0 free parameters · 5 assumptions · 0 invented entities

The theorem is a conditional statement: it holds for models satisfying Assumption 1, and the proof uses standard stochastic calculus plus the stacked coupling from prior work. No parameters are fitted to data, and no new physical entities are introduced.

assumptions (5)

domain assumption Reaction networks are modeled as continuous-time Markov chains on Z^d_{\ge0} with intensity functions λ_k(x,t), via the random time-change representation (6) and space-time Poisson representation (7).
The entire analysis operates inside this modeling framework; the theorem is about these coupled CTMC processes.
domain assumption Assumption 1: linear growth for net-increasing reactions (condition 1), polynomial growth for all reactions (condition 2), and the uniform relative parameter-sensitivity bound (condition 3).
This is the sufficient condition for Theorem 2; stated in §2.4. If any part fails, the proof's lemmas (exponential tail, finite moments of N, small decoupling probability) break.
domain assumption The stacked coupling (9) from Anderson-Yuan (2018) is the coupling mechanism; its construction is used without re-derivation.
The variance bound applies specifically to this coupling. The coupling construction is from prior work by the same authors.
standard math Standard results on Poisson processes, including the Chernoff bound for Poisson tails, and the existence of unit-rate Poisson processes and space-time Poisson point processes.
Used in Lemmas 2, 8, and the proof of Theorem 2.
standard math Basic measure-theoretic probability: integrability and moment calculations, dominated convergence, and related tools.
Used throughout, especially in Lemma 8 and Theorem 2's expectation decomposition.

how reviews work

0 comments

Cite this review

Pith. "Pith review of Variance of finite difference methods for reaction networks with non-Lipschitz rate functions." pith.science (2026). https://pith.science/paper/AG6F5QFJ

@misc{pith2026190806880,
  author       = {Pith},
  title        = {Pith review of: Variance of finite difference methods for reaction networks with non-Lipschitz rate functions},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/AG6F5QFJ}},
  note         = {Machine review of arXiv:1908.06880}
}

read the original abstract

Parametric sensitivity analysis is a critical component in the study of mathematical models of physical systems. Due to its simplicity, finite difference methods are used extensively for this analysis in the study of stochastically modeled reaction networks. Different coupling methods have been proposed to build finite difference estimators, with the "split coupling," also termed the "stacked coupling," yielding the lowest variance in the vast majority of cases. Analytical results related to this coupling are sparse, and include an analysis of the variance of the coupled processes under the assumption of globally Lipschitz intensity functions [Anderson, SIAM Numerical Analysis, Vol. 50, 2012]. Because of the global Lipschitz assumption utilized in [Anderson, SIAM Numerical Analysis, Vol. 50, 2012], the main result there is only applicable to a small percentage of the models found in the literature, and it was conjectured that similar results should hold for a much wider class of models. In this paper we demonstrate this conjecture to be true by proving the variance of the coupled processes scales in the desired manner for a large class of non-Lipschitz models. We further extend the analysis to allow for time dependence in the parameters. In particular, binary systems with or without time-dependent rate parameters, a class of models that accounts for the vast majority of systems considered in the literature, satisfy the assumptions of our theory.

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

28 extracted references · 26 canonical work pages

[1]

Anderson

David F. Anderson. A modiﬁed next reaction method for simulating chemical sys- tems with time dependent propensities and delays. The Journal of chemical physics , 127(21):214107, 2007

work page 2007
[2]

Anderson

David F. Anderson. An eﬃcient ﬁnite diﬀerence method for param eter sensitivities of continuous time Markov chains. SIAM Journal on Numerical Analysis , 50(5):2237 – 2258, 2012

work page 2012
[3]

Anderson and Masanori Koyama

David F. Anderson and Masanori Koyama. An asymptotic relationship between coupling methods for stochastically modeled population processes. IMA Journal of Numerical Analysis, 35(4):1757–1778, 2015

work page 2015
[4]

Anderson and Thomas G

David F. Anderson and Thomas G. Kurtz. Continuous time Markov chain models for chemical reaction networks, chapter 1 in Design and Analysis of Biomolecular Circuits: Engineering Approaches to Systems and Synthetic Biology. Springe r, 2011

work page 2011
[5]

Anderson and Thomas G

David F. Anderson and Thomas G. Kurtz. Stochastic analysis of biochemical systems . Springer, 2015

work page 2015
[6]

Anderson, Timo Sepp¨ al¨ ainen, and Benedek Valk´ o.Introduction to Probability

David F. Anderson, Timo Sepp¨ al¨ ainen, and Benedek Valk´ o.Introduction to Probability. Cambridge University Press, 2017

work page 2017
[7]

Anderson and Elizabeth Skubak Wolf

David F. Anderson and Elizabeth Skubak Wolf. A ﬁnite diﬀerence method for estimating second order parameter sensitivities of discrete stochastic chem ical reaction networks. J. Chem. Phys. , 137(22):224112, 2012. 19

work page 2012
[8]

Anderson and Elizabeth Skubak Wolf

David F. Anderson and Elizabeth Skubak Wolf. Hybrid pathwise sen sitivity methods for discrete stochastic models of chemical reaction systems. J. Chem. Phys. , 142:034103, 2015

work page 2015

Show all 28 references

[9]

Low variance couplings for st ochastic models of intracellular processes with time-dependent rate functions

David F Anderson and Chaojie Yuan. Low variance couplings for st ochastic models of intracellular processes with time-dependent rate functions. Bulletin of mathematical biology, pages 1–29, 2018

2018
[10]

Notes on the sum and maximum of independent e xponen- tially distributed random variables with diﬀerent scale parameters

Markus Bibinger. Notes on the sum and maximum of independent e xponen- tially distributed random variables with diﬀerent scale parameters. arXiv preprint arXiv:1307.3945, 2013

2013 arXiv
[11]

Stochastic gene expression in a single cell

Michael B Elowitz, Arnold J Levine, Eric D Siggia, and Peter S Swain. Stochastic gene expression in a single cell. Science, 297(5584):1183–1186, 2002

2002
[12]

On the stability of stochastic jump kinetics

Stefan Engblom. On the stability of stochastic jump kinetics. Applied Mathematics , 5(14):3217–3239, 2014

2014
[13]

Optimizing genetic circuits by global sensitivity analysis

Xiao-jiang Feng, Sara Hooshangi, David Chen, Genyuan Li, Ron W eiss, and Herschel Rabitz. Optimizing genetic circuits by global sensitivity analysis. Biophysical journal, 87(4):2195–2202, 2004

2004
[14]

A general method for numerically simulating the s tochastic time evolution of coupled chemical reactions

Daniel T Gillespie. A general method for numerically simulating the s tochastic time evolution of coupled chemical reactions. Journal of computational physics , 22(4):403– 434, 1976

1976
[15]

Gradient Estimation Via Perturbation Analysis

Paul Glasserman. Gradient Estimation Via Perturbation Analysis . Kluwer Academic Publishers, 1991

1991
[16]

A scalable computational frame- work for establishing long-term behavior of stochastic reaction ne tworks

Ankit Gupta, Corentin Briat, and Mustafa Khammash. A scalable computational frame- work for establishing long-term behavior of stochastic reaction ne tworks. PLoS compu- tational biology, 10(6):e1003669, 2014

2014
[17]

Unbiased estimation of par ameter sensitivi- ties for stochastic chemical reaction networks

Ankit Gupta and Mustafa Khammash. Unbiased estimation of par ameter sensitivi- ties for stochastic chemical reaction networks. SIAM Journal on Scientiﬁc Computing , 35(6):A2598–A2620, 2013

2013
[18]

Sensitiv- ity, robustness, and identiﬁability in stochastic chemical kinetics mo dels

Micha/suppress l Komorowski, Maria J Costa, David A Rand, and Michael PHStumpf. Sensitiv- ity, robustness, and identiﬁability in stochastic chemical kinetics mo dels. Proceedings of the National Academy of Sciences , 108(21):8645–8650, 2011

2011
[19]

Thomas G. Kurtz. Representations of markov processes as multiparameter time changes. Ann. Prob., 8(4):682–715, 1980

1980
[20]

Summing up the noise in gene networks

Johan Paulsson. Summing up the noise in gene networks. Nature, 427(6973):415, 2004

2004
[21]

Sergey Plyasunov and Adam P. Arkin. Eﬃcient stochastic sensitivity analysis of discrete event systems. J. Comp. Phys. , 221:724 – 738, 2007. 20

2007
[22]

Stochastic mrna synthesis in mammalian cells

Arjun Raj, Charles S Peskin, Daniel Tranchina, Diana Y Vargas, and Sanjay Tyagi. Stochastic mrna synthesis in mammalian cells. PLoS biology, 4(10):e309, 2006

2006
[23]

Moment growth bounds on continuous time m arkov processes on non-negative integer lattices

Muruhan Rathinam. Moment growth bounds on continuous time m arkov processes on non-negative integer lattices. Quarterly of Applied Mathematics , 73(2), 2013

2013
[24]

Sheppard, and Mustafa Khamma sh

Muruhan Rathinam, Patrick W. Sheppard, and Mustafa Khamma sh. Eﬃcient com- putation of parameter sensitivities of discrete stochastic chemica l reaction networks. Journal of Chemical Physics , 132:034103, 2010

2010
[25]

Anderson, and James B

Rishi Srivastava, David F. Anderson, and James B. Rawlings. Co mparison of ﬁnite diﬀerence based methods to obtain sensitivities of stochastic chem ical kinetic models. Journal of Chemical Physics , 138:074110, 2013

2013
[26]

Robustne ss properties of circadian clock architectures

J¨ org Stelling, Ernst Dieter Gilles, and Francis J Doyle. Robustne ss properties of circadian clock architectures. Proceedings of the National Academy of Sciences , 101(36):13210–13215, 2004

2004
[27]

Eﬃcient ﬁnite diﬀerence method for computing sensitivities of biochemical reactions

Vo Hong Thanh, Roberto Zunino, and Corrado Priami. Eﬃcient ﬁnite diﬀerence method for computing sensitivities of biochemical reactions. Proceedings of the Royal Society A , 474:20180303, 2018

2018
[28]

Stochastic modelling for systems biology

Darren J Wilkinson. Stochastic modelling for systems biology . Chapman and Hall/CRC, 2006. 21

2006

Pith tools

Reviewed August 14, 2026 · model on record in the stance chip above.

[1] [1]

Anderson

David F. Anderson. A modiﬁed next reaction method for simulating chemical sys- tems with time dependent propensities and delays. The Journal of chemical physics , 127(21):214107, 2007

work page 2007

[2] [2]

Anderson

David F. Anderson. An eﬃcient ﬁnite diﬀerence method for param eter sensitivities of continuous time Markov chains. SIAM Journal on Numerical Analysis , 50(5):2237 – 2258, 2012

work page 2012

[3] [3]

Anderson and Masanori Koyama

David F. Anderson and Masanori Koyama. An asymptotic relationship between coupling methods for stochastically modeled population processes. IMA Journal of Numerical Analysis, 35(4):1757–1778, 2015

work page 2015

[4] [4]

Anderson and Thomas G

David F. Anderson and Thomas G. Kurtz. Continuous time Markov chain models for chemical reaction networks, chapter 1 in Design and Analysis of Biomolecular Circuits: Engineering Approaches to Systems and Synthetic Biology. Springe r, 2011

work page 2011

[5] [5]

Anderson and Thomas G

David F. Anderson and Thomas G. Kurtz. Stochastic analysis of biochemical systems . Springer, 2015

work page 2015

[6] [6]

Anderson, Timo Sepp¨ al¨ ainen, and Benedek Valk´ o.Introduction to Probability

David F. Anderson, Timo Sepp¨ al¨ ainen, and Benedek Valk´ o.Introduction to Probability. Cambridge University Press, 2017

work page 2017

[7] [7]

Anderson and Elizabeth Skubak Wolf

David F. Anderson and Elizabeth Skubak Wolf. A ﬁnite diﬀerence method for estimating second order parameter sensitivities of discrete stochastic chem ical reaction networks. J. Chem. Phys. , 137(22):224112, 2012. 19

work page 2012

[8] [8]

Anderson and Elizabeth Skubak Wolf

David F. Anderson and Elizabeth Skubak Wolf. Hybrid pathwise sen sitivity methods for discrete stochastic models of chemical reaction systems. J. Chem. Phys. , 142:034103, 2015

work page 2015