Tau Leaping and Reaction Diffusion Model Error — Examples

By [F3], the fine count is $P_c+P_f$ and the coarse count is $P_c$, where the three count components have the independent laws in the Example. Their marginals have parameters $13.1h$ and $13h$ because independent Poisson counts add; this also follows by convolving the masses in [F1].

The coupled difference is exactly $P_f$, so its variance is $0.1h$ by [F2]. If the two counts were independent, their difference would have variance $13.1h+13h=26.1h$. Thus the shared minimum reduces the correction variance by a factor of $261$ in this one-step example.

The leap state is negative exactly when $1-K<0$, equivalently $K\ge 2$. Using the masses in [F3], $\mathbb{P}(K\ge 2)=1-\mathbb{P}(K=0)-\mathbb{P}(K=1)=1-e^{-ch}(1+ch)>0.$ The strict inequality follows from $e^{ch}>1+ch$ for $ch>0$.

For small $h$, Taylor expansion of $e^{-ch}(1+ch)$ gives $\mathbb{P}(K\ge 2)=\frac{1}{2}{c}^2{h}^2+O(h^3)$. The exact chain cannot become negative because its only positive-rate event sends $1$ to $0$, after which the rate is zero by [F1]. Thus the invalid tau-leap event is a distinct numerical exit from the physical state set.

By [F1], $N_T\sim \mathrm{Pois}(m)$, so $X_T=N_T\mathrm{/}V\ge 0$ almost surely. The point-process counts on disjoint intervals are independent by [F2]. On an interval of length $h$, put $r=V\lambda h$. By [F3], $\mathbb{P}(K=0)=e^{-r}=1-r+O(r^2)$ and $\mathbb{P}(K=1)=r{e}^{-r}=r+O(r^2)$; also $\mathbb{P}(K\ge 2)\le \mathbb{E}[K(K-1)]\mathrm{/}2=r^2\mathrm{/}2$. For bounded $f$, the multiple-count contribution to its expected increment is therefore $O(h^2)\mathrm{\parallel }f{\mathrm{\parallel }}_{\mathrm{\infty }}$. Dividing the zero/single-count expansion by $h$ gives the density-count generator $V\lambda (f(x+1\mathrm{/}V)-f(x))$, as in [F4].

Let $B$ be the Brownian motion in [F6] and define $Z_t=\lambda t+\sqrt{\lambda \mathrm{/}V}{B}_t$. By [F7] and the independent increments in [F6], $Z$ is a continuous Markov process. For any $f$ with three bounded continuous derivatives, put ${\mathrm{\Delta }}_h=\lambda h+\sqrt{\lambda \mathrm{/}V}(B_{t+h}-B_t)$. The Brownian increment is independent of the past and has variance $h$. The density in [F5] has finite third absolute moment by substituting $v=u^2\mathrm{/}2$, so $\mathbb{E}[{\mathrm{\Delta }}_h]=\lambda h,\mathbb{E}[{\mathrm{\Delta }}_h^2]=\frac{\lambda h}{V}+{\lambda }^2{h}^2,\mathbb{E}\mathrm{\mid }{\mathrm{\Delta }}_h{\mathrm{\mid }}^3=O(h^{3\mathrm{/}2})\mathrm{.}$ Taylor's formula therefore gives ${\lim }_{h\downarrow 0}\frac{\mathbb{E}[f(x+{\mathrm{\Delta }}_h)]-f(x)}{h}=\lambda f^{\mathrm{\prime }}(x)+\frac{\lambda }{2V}{f}^{\mathrm{\prime }\mathrm{\prime }}(x),$ the generator in the Example. Also $B_T\sim N(0,T)$ by [F6], so $G=B_T\mathrm{/}\sqrt{T}\sim N(0,1)$ and $Z_T=\mu +\sqrt{\mu \mathrm{/}V}G$.

The density $\phi$ in [F5] is even. The moments of orders one and three are absolutely integrable, as follows by substituting $v=u^2\mathrm{/}2$ in the tail integrals, so symmetry gives $\mathbb{E}G=\mathbb{E}{G}^3=0$. Since ${\phi }^{\mathrm{\prime }}(u)=-u\phi (u)$ and $u\phi (u)\to 0$ at both infinities, integration by parts gives $\mathbb{E}[G^2]={\int }_{\mathbb{R}}{u}^2\phi (u)du={\int }_{\mathbb{R}}\phi (u)du=1.$ Expanding the Gaussian marginal in step 1.2 now gives $\mathbb{E}[Z_T^3]={\mu }^3+3{\mu }^2\mathrm{/}V\mathrm{.}$

Reindexing the nonnegative Poisson mass series gives $\mathbb{E}[N_T(N_T-1)(N_T-2)]=e^{-m}{\sum }_{k=3}^{\mathrm{\infty }}\frac{k(k-1)(k-2)m^k}{k!}=m^3{e}^{-m}{\sum }_{r=0}^{\mathrm{\infty }}\frac{m^r}{r!}=m^3\mathrm{.}$ Using [F3] and $N^3=N(N-1)(N-2)+3N(N-1)+N$ yields $\mathbb{E}[X_T^3]={\mu }^3+\frac{3{\mu }^2}{V}+\frac{\mu }{V^2}\mathrm{.}$

Subtracting the two moments in steps 2.2 and 2.1 gives $\mathbb{E}[X_T^3]-\mathbb{E}[Z_T^3]=\mu \mathrm{/}{V}^2>0.$ Thus this exact marginal comparison has a weak discrepancy of order $V^{-2}$ for the cubic payoff.

The CLE value is negative exactly when $G<-\sqrt{V\mu }$. The normal density is strictly positive on every nonempty interval, so this event has positive probability. As $V\mu \downarrow 0$, continuity of the normal distribution and symmetry give $\mathbb{P}(Z_T<0)=\mathbb{P}(G<-\sqrt{V\mu })⟶\frac{1}{2}\mathrm{.}$ The jump marginal remains nonnegative by step 1.1.

Apply [F2] to $v(t,x)=x$ and use [F1]. The exact mean $m(t)=\mathbb{E}{X}_t$ satisfies $m(t)=x_0+{\int }_0^t(c_{2}M-\lambda m(s))ds\mathrm{.}$ This makes $m$ absolutely continuous; differentiating almost everywhere shows that the derivative of $e^{\lambda t}(m(t)-c_{2}M\mathrm{/}\lambda )$ is zero almost everywhere. This product is absolutely continuous, so it is constant; its value at $t=0$ is $x_0-c_{2}M\mathrm{/}\lambda$. Rearranging gives $m(t)=\frac{c_{2}M}{\lambda }+(x_0-\frac{c_{2}M}{\lambda })e^{-\lambda t}\mathrm{.}$

By [F3], one Poisson leap has $\mathbb{E}{Y}_1=(1-\lambda h)x_0+c_{2}Mh$. Its linear mean amplification factor is $1-\lambda h$, whose absolute value exceeds $1$ when $h>2\mathrm{/}\lambda$. For the symmetric case $c_1=c_2=r$ and initial state $x_0=0$, the choice $h=3\mathrm{/}(2r)$ gives $\mathbb{E}{Y}_1=3M\mathrm{/}2$, outside the physical upper bound $M$, whereas the exact expectation from step 1.1 is $\frac{M}{2}(1-e^{-3})<M\mathrm{/}2$. This is a one-step comparison; it does not assume that invalid Poisson leaps remain in the physical state set.

Since $n>n_c$, the death channel is noncritical, and $\tau =\min (T,{\tau }_{\max },m\mathrm{/}a(n))=h$. The critical set is empty, so the algorithm skips the critical SSA call, records no critical event, and sets $\delta =h$. Thus $p=1\mathrm{/}2$. The guarded count $K_B$ is a sum of $n$ Bernoulli variables, each in $\{0,1\}$ by [F3], so $0\le K_B\le n$ and the updated population $n-K_B$ is nonnegative. Its mean is $np=n\mathrm{/}2$.

Independence of the Bernoulli trials gives $\mathrm{Var}(K_B)={\sum }_{i=1}^{n}p(1-p)=n\mathrm{/}4$. In contrast, the Poisson count has mean $a(n)h=n\mathrm{/}2$ and variance $n\mathrm{/}2$ by [F1, F4]. Thus the binomial guard retains the mean but halves the count variance in this finite model; it is not the same Poisson leap law.

The unguarded Poisson count has $\mathbb{P}(K_P>n)>0$, since every mass at $k>n$ is positive when its parameter $n\mathrm{/}2$ is positive. Therefore its leap can overdraw the available molecules, while the guarded count cannot.