Tau Leaping and Reaction Diffusion Model Error

If $A\sim \mathrm{Pois}(\alpha )$ and $B\sim \mathrm{Pois}(\beta )$ are independent, then for $\alpha ,\beta >0$ and $n\in {\mathbb{N}}_0$ their convolution is $\mathbb{P}(A+B=n)=e^{-(\alpha +\beta )}{\sum }_{k=0}^n\frac{{\alpha }^k{\beta }^{n-k}}{k!(n-k)!}=e^{-(\alpha +\beta )}\frac{(\alpha +\beta {)}^n}{n!},$ by the binomial theorem. If either parameter is zero, the corresponding count is identically zero by [F1], and the same conclusion follows. Thus an independent sum has Poisson law with the sum of its parameters.

Apply step 1.1 to $C+R_1$ and $C+R_2$. Their parameters are respectively $h\min ({\lambda }_1,{\lambda }_2)+h({\lambda }_1-{\lambda }_2{)}^{+}=h{\lambda }_1$ and $h\min ({\lambda }_1,{\lambda }_2)+h({\lambda }_2-{\lambda }_1{)}^{+}=h{\lambda }_2$. This proves the marginal laws.

The shared count cancels in $K_1-K_2=R_1-R_2$. Independence and [F2] give its mean as $h({\lambda }_1-{\lambda }_2{)}^{+}-h({\lambda }_2-{\lambda }_1{)}^{+}=h({\lambda }_1-{\lambda }_2)$ and its variance as $h({\lambda }_1-{\lambda }_2{)}^{+}+h({\lambda }_2-{\lambda }_1{)}^{+}=h\mathrm{\mid }{\lambda }_1-{\lambda }_2\mathrm{\mid }$. Also, since $R_1,R_2$ are independent of $C$ and of each other, $\mathrm{Cov}(K_1,K_2)=\mathrm{Var}(C)=h\min ({\lambda }_1,{\lambda }_2)$ by [F2]. This proves all identities.

Uniformly over $x\in D\cup F$, the no-event probability in a time interval of length $\delta$ is $e^{-\lambda (x)\delta }$, where $\lambda (x)=0$ on $F$ and $\lambda (x)\le A$ on $D$. Its difference from $1-\lambda (x)\delta$ is at most $A^2{\delta }^2\mathrm{/}2$. For a channel $j$ from $x\in D$, the probability that it fires once at time $r\in (0,\delta )$ and no later event occurs is ${\int }_0^{\delta }{a}_j(x)e^{-\lambda (x)r}{e}^{-\lambda (x+{\nu }_j)(\delta -r)}dr\mathrm{.}$ Both rates in the exponent are at most $A$, so its difference from $a_j(x)\delta$ is at most $a_j(x)A{\delta }^2$. Summing over channels gives at most $\lambda (x)A{\delta }^2\le A^2{\delta }^2$. After the first event, the conditional probability of another before $\delta$ is at most $A(\delta -r)$; integrating against the first-event rate, at most $A$, shows that the probability of at least two events is at most $A^2{\delta }^2\mathrm{/}2$. These estimates are uniform on the finite state set.

Since $D\cup F$ is finite and ${\mathrm{\partial }}_{t}v$ is continuous, uniformly in $s$ and $x$, $v(s+\delta ,x)=v(s,x)+\delta {\mathrm{\partial }}_{s}v(s,x)+o(\delta )\mathrm{.}$ The zero-event probability is $1-\lambda (x)\delta +O({\delta }^2)$; the single-channel probabilities are $a_j(x)\delta +O({\delta }^2)$ by step 1.1; and the multiple-event probability is $O({\delta }^2)$. In the zero-event case, the contribution is $v(s,x)+\delta ({\mathrm{\partial }}_{s}v(s,x)-\lambda (x)v(s,x))+o(\delta )$. The one-event channel contributions are $\delta {\sum }_j{a}_j(x)v(s,x+{\nu }_j)+o(\delta )$, because the time change in $v(s+\delta ,x+{\nu }_j)$ is multiplied by an $O(\delta )$ event probability. The multiple-event contribution is $O({\delta }^2)$. Adding the cases, the order-$\delta$ terms are $\delta {\mathrm{\partial }}_{s}v(s,x)+\delta {\sum }_j{a}_j(x)(v(s,x+{\nu }_j)-v(s,x)),$ with a uniform $o(\delta )$ remainder. This is $\delta ({\mathrm{\partial }}_{s}v(s,x)+Qv(s,\cdot )(x))+o(\delta )$; for $x\in F$ the chain is fixed and $Q=0$. Multiple-event payoffs contribute $O({\delta }^2)$ because $v$ is bounded on the finite set.

Sum the conditional identity in step 2.1 over a deterministic partition of $[0,t]$ and take expectations. The uniform $o(\delta )$ remainders sum to $o(1)$. By [F3] the chain has finitely many jumps almost surely on $[0,t]$; therefore $s\mapsto {\mathrm{\partial }}_{s}v(s,X_s)+Qv(s,\cdot )(X_s)$ is bounded and piecewise continuous, so its Riemann sums converge pathwise. Boundedness and dominated convergence let expectations pass to the integral. The left side telescopes to $\mathbb{E}[v(t,X_t)]-\mathbb{E}[v(0,X_0)]$, proving the identity.

Each parameter $a_j(x)h$ is finite, so [F1] defines the corresponding Poisson law and [F2] gives $\mathbb{E}[K_j\mid x]=a_j(x)h$. By [F3], conditional expectation is linear, hence $\mathbb{E}[\mathrm{\Delta }Y\mid x]={\sum }_j{\nu }_j{a}_j(x)h$.

Conditional independence in [F3] makes cross-channel covariances zero. The covariance of ${\nu }_j{K}_j$ is ${\nu }_j{\nu }_j^{\mathsf{T}}\mathrm{Var}(K_j\mid x)=h{a}_j(x){\nu }_j{\nu }_j^{\mathsf{T}}$ by [F2]. Summing proves the covariance formula.

For a square-integrable vector, $\mathbb{E}\mathrm{\parallel }\mathrm{\Delta }Y{\mathrm{\parallel }}^2=\mathrm{tr}\mathrm{Cov}(\mathrm{\Delta }Y)+\mathrm{\parallel }\mathbb{E}\mathrm{\Delta }Y{\mathrm{\parallel }}^2$. Substitution of steps 1.1 and 1.2 gives the asserted identity.

For the Poisson leap, the probability of zero total counts is $e^{-\lambda (x)h}$. The probability that exactly channel $j$ fires once is $a_j(x)h{e}^{-\lambda (x)h}$, so its difference from $a_j(x)h$ is at most $a_j(x)\lambda (x)h^2$. Summing these errors over $j$ gives at most $\lambda (x{)}^2{h}^2\le A^2{h}^2$. If $N={\sum }_j{K}_j$, independence and [F5] give $\mathbb{E}[N(N-1)\mid x]={\sum }_j(a_j(x)h{)}^2+2{\sum }_{i<j}{a}_i(x)a_j(x)h^2=(\lambda (x)h{)}^2,$ hence $\mathbb{P}(N\ge 2\mid x)\le (\lambda (x)h{)}^2\mathrm{/}2\le A^2{h}^2\mathrm{/}2$. The zero-count probability differs from $1-\lambda (x)h$ by at most $A^2{h}^2\mathrm{/}2$. A single allowed channel event remains in $E$ by [F1], so only the $N\ge 2$ part can trigger the invalid-state branch. The zero-count, total single-count, and multiple-count contributions are bounded respectively by $G{A}^2{h}^2\mathrm{/}2$, $G{A}^2{h}^2$, and $G{A}^2{h}^2\mathrm{/}2$. Thus $\mathrm{\mid }(K_{h}g)(x)-g(x)-h(Qg)(x)\mathrm{\mid }\le 2G{A}^2{h}^2$.

For the exact chain, the no-event probability is $e^{-\lambda (x)h}$, which differs from $1-\lambda (x)h$ by at most $A^2{h}^2\mathrm{/}2$. The probability of one channel-$j$ event and no later event is ${\int }_0^h{a}_j(x)e^{-\lambda (x)s}{e}^{-\lambda (x+{\nu }_j)(h-s)}ds\mathrm{.}$ Its difference from $a_j(x)h$ is at most $a_j(x)A{h}^2$, since all total rates on $E$ are at most $A$; summing over $j$ gives at most $\lambda (x)A{h}^2\le A^2{h}^2$. The probability of two or more events is at most $A^2{h}^2\mathrm{/}2$: conditional on a first event at time $s$, a later event has probability at most $A(h-s)$, and integration against the first-event rate bounded by $A$ gives the bound. The no-event, total single-event, and multiple-event contributions are bounded respectively by $G{A}^2{h}^2\mathrm{/}2$, $G{A}^2{h}^2$, and $G{A}^2{h}^2\mathrm{/}2$. Therefore $\mathrm{\mid }(Q_{h}g)(x)-g(x)-h(Qg)(x)\mathrm{\mid }\le 2G{A}^2{h}^2$.

Subtract the two first-order expansions from steps 1.1 and 1.2 and use the triangle inequality. This gives $\mathrm{\mid }(K_{h}g)(x)-(Q_{h}g)(x)\mathrm{\mid }\le 4G{A}^2{h}^2$ for every $x\in E$. Outside $E$ both kernels hold the state fixed, so the same bound holds there. Taking the supremum proves the final inequality.

Fix an interior lattice state $x$ and a time $s$. For channel $j$, apply one-variable Taylor's theorem to $r\mapsto u(s,x+r{\nu }_j\mathrm{/}V)$ on $[0,1]$. Its third-order remainder has absolute value at most $M\mathrm{\parallel }{\nu }_j{\mathrm{\parallel }}^3\mathrm{/}(6V^3)$. Multiplying by $V{a}_j(x)$ and summing over channels gives ${\mathcal{L}}_{V}u(s,x)={\sum }_j{a}_j(x){\nu }_j\cdot \mathrm{\nabla }u(s,x)+\frac{1}{2V}{\sum }_j{a}_j(x){\nu }_j^{\mathsf{T}}{D}^{2}u(s,x){\nu }_j+R_V(s,x),$ where $\mathrm{\mid }{R}_V(s,x)\mathrm{\mid }\le \frac{M}{6V^2}{\sum }_j{A}_j\mathrm{\parallel }{\nu }_j{\mathrm{\parallel }}^3\mathrm{.}$ The first two terms are exactly ${\mathcal{G}}_{V}u(s,x)$ by [F1].

Set ${\mathrm{\Psi }}_V(x)=V(x-x_0)$, $D_V={\mathrm{\Psi }}_V(D)$, and let $F_V$ be the images under ${\mathrm{\Psi }}_V$ of the one-jump exit destinations from $D$. These are disjoint finite subsets of ${\mathbb{Z}}^d$. The transported stopped chain $\tilde{X}={\mathrm{\Psi }}_V(\bar{X})$ moves from $y\in D_V$ by channel increment ${\nu }_j$ at rate $V{a}_j(x_0+y\mathrm{/}V)$ and is fixed on $F_V$. For $\tilde{u}(s,y)=u(s,x_0+y\mathrm{/}V)$, its time derivative is unchanged and its reaction generator on $D_V$ equals ${\mathcal{L}}_{V}u(s,x_0+y\mathrm{/}V)$. Apply [F2] to this integer-lattice chain and $\tilde{u}$. Before ${\tau }_X$, the integrand is ${\mathrm{\partial }}_{s}u+{\mathcal{L}}_{V}u=({\mathcal{L}}_V-{\mathcal{G}}_V)u$ by the backward equation [A3]. After exit, $\bar{X}$ is fixed at an exit destination where $u(s,{\bar{X}}_s)=g({\bar{X}}_s)$ is independent of $s$, so both its time derivative and stopped generator vanish. At time $T$, $u(T,{\bar{X}}_T)=g({\bar{X}}_T)$, and the common initial condition gives ${\bar{X}}_0={\bar{Z}}_0=x_0$ and $u(0,x_0)=\mathbb{E}g({\bar{Z}}_T)$. Therefore $\mathbb{E}g({\bar{X}}_T)-\mathbb{E}g({\bar{Z}}_T)=\mathbb{E}{\int }_0^{T\wedge {\tau }_X}({\mathcal{L}}_V-{\mathcal{G}}_V)u(s,X_s^V)ds\mathrm{.}$

Bound the integrand in step 1.2 by step 1.1 and use $T\wedge {\tau }_X\le T$. This proves the stopped-process estimate. Under [A4], the stopped and original terminal payoffs agree unless their respective process exits by time $T$; on that event their difference is at most $2G$. The triangle inequality for the two expectations adds at most $2G\mathbb{P}({\tau }_X\le T)+2G\mathbb{P}({\tau }_Z\le T)$, which proves the unstopped bound.

For lattice states $x,z$, write $e=x-z$ and ${\delta }_j=a_j(x)-a_j(z)$. Bounded rates and the integer lattice's minimum nonzero distance one imply $\mathrm{\mid }{\delta }_j\mathrm{\mid }\le L\mathrm{\parallel }e\mathrm{\parallel }$ for a constant $L$; this also holds when one state is outside $E$, where the extended rate is zero. In one half interval, shared counts cancel from the increment difference $\eta$. The independent residual Poisson counts and [F5] give $m:=\mathbb{E}[\eta \mid x,z]=h{\sum }_j{\nu }_j{\delta }_j,\mathbb{E}[\mathrm{\parallel }\eta -m{\mathrm{\parallel }}^2\mid x,z]=h{\sum }_j\mathrm{\mid }{\delta }_j\mathrm{\mid }\mathrm{\parallel }{\nu }_j{\mathrm{\parallel }}^2\mathrm{.}$ Thus $\mathrm{\parallel }m\mathrm{\parallel }\le Ch\mathrm{\parallel }e\mathrm{\parallel }$ and the centered second moment is at most $Ch\mathrm{\parallel }e\mathrm{\parallel }$. If $e=0$, then $\eta =0$. Otherwise $\mathrm{\parallel }e\mathrm{\parallel }\ge 1$, so $h\mathrm{\parallel }e\mathrm{\parallel }\le h\mathrm{\parallel }e{\mathrm{\parallel }}^2$. For $h\le 1$ this proves $\mathbb{E}[\mathrm{\parallel }\eta {\mathrm{\parallel }}^2\mid x,z]\le Ch\mathrm{\parallel }e{\mathrm{\parallel }}^2\mathrm{.}$

If $x\mathrm{\notin }E$, the invalid-state convention freezes the fine path, so its increment $R$ is zero and the estimates below are immediate. Suppose $x\in E$, and let $R$ be its fine-path increment on the first half interval. By [F3], [F4], and the triangle inequality, $\mathbb{E}\mathrm{\parallel }R\mathrm{\parallel }\le {\sum }_j\mathrm{\parallel }{\nu }_j\mathrm{\parallel }h{a}_j(x)\le Ch\mathrm{.}$ By [F2] and bounded propensities, $\mathbb{E}\mathrm{\parallel }R{\mathrm{\parallel }}^2\le Ch$ for $h\le 1$. If $x^{\mathrm{\prime }}=x+R$, then $\mathbb{E}\mathrm{\parallel }{x}^{\mathrm{\prime }}-z\mathrm{\parallel }\le \mathrm{\parallel }e\mathrm{\parallel }+Ch,\mathbb{E}\mathrm{\parallel }{x}^{\mathrm{\prime }}-z{\mathrm{\parallel }}^2\le C(\mathrm{\parallel }e{\mathrm{\parallel }}^2+h)\mathrm{.}$ For the second inequality, expand $\mathrm{\parallel }e+R{\mathrm{\parallel }}^2$ and use $\mathrm{\parallel }\mathbb{E}R\mathrm{\parallel }\le Ch$ and $\mathbb{E}\mathrm{\parallel }R{\mathrm{\parallel }}^2\le Ch$; the cross term is bounded by $Ch\mathrm{\parallel }e\mathrm{\parallel }\le Ch\mathrm{\parallel }e{\mathrm{\parallel }}^2$ when $e\mathrm{\ne }0$, and is zero when $e=0$.

Write the two half-interval differences as ${\eta }_i=m_i+{\xi }_i$, where $m_i$ is the conditional mean and ${\xi }_i$ is centered given the starting states of that half interval. Step 1.1 gives $\mathrm{\parallel }{m}_1\mathrm{\parallel }\le Ch\mathrm{\parallel }e\mathrm{\parallel },\mathbb{E}\mathrm{\parallel }{\xi }_1{\mathrm{\parallel }}^2\le Ch\mathrm{\parallel }e{\mathrm{\parallel }}^2\mathrm{.}$ The second half starts at $x^{\mathrm{\prime }}$ and $z$, so $\mathrm{\parallel }{m}_2\mathrm{\parallel }\le Ch\mathrm{\parallel }{x}^{\mathrm{\prime }}-z\mathrm{\parallel }$. Steps 1.1 and 1.2 therefore give $\mathbb{E}\mathrm{\parallel }{m}_2{\mathrm{\parallel }}^2\le C{h}^2(\mathrm{\parallel }e{\mathrm{\parallel }}^2+h),\mathbb{E}\mathrm{\parallel }{\xi }_2{\mathrm{\parallel }}^2\le Ch(\mathrm{\parallel }e\mathrm{\parallel }+h)\le Ch\mathrm{\parallel }e{\mathrm{\parallel }}^2+C{h}^2\mathrm{.}$ For the latter inequality, use $\mathbb{E}\mathrm{\parallel }{x}^{\mathrm{\prime }}-z\mathrm{\parallel }\le \mathrm{\parallel }e\mathrm{\parallel }+Ch$ and $\mathrm{\parallel }e\mathrm{\parallel }\le \mathrm{\parallel }e{\mathrm{\parallel }}^2$ when $e\mathrm{\ne }0$. Let ${\mathcal{H}}_1$ be the complete first-half history. The fresh second-half copies are independent of that history, so $\mathbb{E}[{\xi }_2\mid {\mathcal{H}}_1]=0$ and $\mathbb{E}[{\xi }_1\cdot {\xi }_2]=0$.

Fix the starting pair $(x,z)$ during this expansion; the estimates are uniform in the pair and depend on it only through $e=x-z$. The endpoint error is $e^{\mathrm{\prime }}=e+m_1+{\xi }_1+m_2+{\xi }_2$. First, $\mathbb{E}\mathrm{\parallel }e+{\eta }_1{\mathrm{\parallel }}^2=\mathrm{\parallel }e+m_1{\mathrm{\parallel }}^2+\mathbb{E}\mathrm{\parallel }{\xi }_1{\mathrm{\parallel }}^2\le (1+Ch)\mathrm{\parallel }e{\mathrm{\parallel }}^2\mathrm{.}$ The vector $e+m_1+{\xi }_1+m_2$ is ${\mathcal{H}}_1$-measurable, so $\mathbb{E}[{\xi }_2\mid {\mathcal{H}}_1]=0$ makes its expected inner product with ${\xi }_2$ zero. The remaining cross terms satisfy $2\mathrm{\mid }\mathbb{E}({\xi }_1\cdot m_2)\mathrm{\mid }\le 2(\mathbb{E}\mathrm{\parallel }{\xi }_1{\mathrm{\parallel }}^2{)}^{1\mathrm{/}2}(\mathbb{E}\mathrm{\parallel }{m}_2{\mathrm{\parallel }}^2{)}^{1\mathrm{/}2}\le C{h}^{3\mathrm{/}2}\mathrm{\parallel }e\mathrm{\parallel }\sqrt{\mathrm{\parallel }e{\mathrm{\parallel }}^2+h}\le Ch\mathrm{\parallel }e{\mathrm{\parallel }}^2,$ with value zero if $e=0$, and $2\mathrm{\mid }\mathbb{E}((e+m_1)\cdot m_2)\mathrm{\mid }\le h\mathrm{\parallel }e+m_1{\mathrm{\parallel }}^2+h^{-1}\mathbb{E}\mathrm{\parallel }{m}_2{\mathrm{\parallel }}^2\le Ch\mathrm{\parallel }e{\mathrm{\parallel }}^2+C{h}^2\mathrm{.}$ The additional terms satisfy $\mathbb{E}\mathrm{\parallel }{m}_2{\mathrm{\parallel }}^2\le Ch\mathrm{\parallel }e{\mathrm{\parallel }}^2+C{h}^2$ and $\mathbb{E}\mathrm{\parallel }{\xi }_2{\mathrm{\parallel }}^2\le Ch\mathrm{\parallel }e{\mathrm{\parallel }}^2+C{h}^2$ by step 2.1. Combining these expansion terms proves the recurrence conditional on $(x,z)$; averaging over starting pairs with difference $e$ gives $\mathbb{E}[\mathrm{\parallel }{e}^{\mathrm{\prime }}{\mathrm{\parallel }}^2\mid e]\le (1+Ch)\mathrm{\parallel }e{\mathrm{\parallel }}^2+C{h}^2\mathrm{.}$

Let $u_n=\mathbb{E}\mathrm{\parallel }{e}_n{\mathrm{\parallel }}^2$ at coarse endpoints. Since $u_0=0$, step 3.1 gives $u_{n+1}\le (1+Ch)u_n+C{h}^2$. There are $N=T\mathrm{/}(2h)$ coarse intervals, so iteration yields $u_N\le C{h}^2{\sum }_{r=0}^{N-1}(1+Ch{)}^r\le C{h}^{2}N{e}^{ChN}\le C_{T}h\mathrm{.}$