Ginwidth=totalheight=keepaspectratio

# Inferring stochastic rates from heterogeneous snapshots of particle
positions

Christopher E. Miles Department of Mathematics, University of California,
Irvine<EMAIL_ADDRESS>0000-0001-5494-403X Scott A. McKinley Department
of Mathematics, Tulane University<EMAIL_ADDRESS>0000-0001-9434-9163 Fangyuan Ding Department of Biomedical Engineering,
University of California, Irvine.  0000-0003-0118-5441 Richard B. Lehoucq
Discrete Math and Optimization, Sandia National Laboratories

(November 3, 2023)

###### Abstract

Many imaging techniques for biological systems – like fixation of cells
coupled with fluorescence microscopy – provide sharp spatial resolution in
reporting locations of individuals at a single moment in time but also destroy
the dynamics they intend to capture. These _snapshot observations_ contain no
information about individual trajectories, but still encode information about
movement and demographic dynamics, especially when combined with a well-
motivated biophysical model. The relationship between spatially evolving
populations and single-moment representations of their collective locations is
well-established with partial differential equations (PDEs) and their inverse
problems. However, experimental data is commonly a set of locations whose
number is insufficient to approximate a continuous-in-space PDE solution.
Here, motivated by popular subcellular imaging data of gene expression, we
embrace the stochastic nature of the data and investigate the mathematical
foundations of parametrically inferring demographic rates from snapshots of
particles undergoing birth, diffusion, and death in a nuclear or cellular
domain. Toward inference, we rigorously derive a connection between individual
particle paths and their presentation as a Poisson spatial process. Using this
framework, we investigate the properties of the resulting inverse problem and
study factors that affect quality of inference. One pervasive feature of this
experimental regime is the presence of cell-to-cell heterogeneity. Rather than
being a hindrance, we show that cell-to-cell geometric heterogeneity can
increase the quality of inference on dynamics for certain parameter regimes.
Altogether, the results serve as a basis for more detailed investigations of
subcellular spatial patterns of RNA molecules and other stochastically
evolving populations that can only be observed for single instants in their
time evolution.

###### Contents

  1. 1 Introduction
    1. 1.1 Background and Motivation
    2. 1.2 Overview of work
  2. 2 Mathematical and Statistical Methods
    1. 2.1 Notation
    2. 2.2 Stochastic model for snapshot observations within a single cell
    3. 2.3 Framework for inference
  3. 3 Inferring dynamics of the 1D RNA model
    1. 3.1 Likelihood and information from analytical solution
    2. 3.2 Bayesian inference on synthetic data
    3. 3.3 Cell-to-cell heterogeneities and their impact on inference
  4. 4 Support for the inference framework.
    1. 4.1 Proof of Theorem 2.1
    2. 4.2 Nondimensionalization
    3. 4.3 Derivation of the likelihood function
    4. 4.4 Spatial binning encodes less statistical information
  5. 5 Discussion
  6. 6 Acknowledgements
  7. A Numerical implementation details
    1. A.1 Stochastic simulation algorithm for synthetic data
    2. A.2 Metropolis-in-Gibbs sampling for Bayesian inference

##  1 Introduction

###  1.1 Background and Motivation

Advances in microscopy and automated tracking have matured to an age of
observing of the spatial evolution of individuals at scales ranging from
tissue-scale collectives of cells [CLL14] down to single molecules [Moe07,
McL+20]. Such spatial tracking has led to discoveries in the understanding of
viral transmission [Che+14], intracellular transport [Osu+19, Ray+23], and
many other areas in cell and molecular biology [MG15]. In these
investigations, the dynamics of a population are probed through the tracked
motion of individual trajectories. The theory and practice of analyzing these
single-particle trajectories is rich and well-developed [QSE91, Mel+16, RW23].
In contrast, many imaging approaches require the fixation of cells (which
requires killing the tissue and cryo-preservation to stabilize agent
locations) and therefore provide only a snapshot of the population at one
moment of time [Bis+21]. Without direct access to the evolution of creation,
destruction, and motion of this population prior to imaging, what can be
inferred about these underlying dynamics?

One very active area of inferring dynamics from snapshots is in the
quantification of gene expression from spatial transcriptomics imaging of
individual RNA molecules [Wei+18, Mar21]. Although time-lapse imaging and
inference of these systems is possible in certain circumstances [Bow+22],
inferring the dynamics underlying populations from RNA snapshots is a far
wider and mainstream interest. Such spatial transcriptomics have served to be
invaluable in recent identifications of cell types [MP22] and disease
mechanisms [Wil+22]. However, the theory to analyze this data has lagged
behind the sophisticated techniques used to harvest it. Spatial imaging
techniques (such as smFISH, single-molecule FISH [Sha+13, Ome+19, DE19]) are
capable of resolving individual RNA molecules, multiplexed over several genes
[May+20] and many cells[CFX06, Eng+19]. This unprecedented single-molecule
resolution is largely neglected, with most analysis approaches focusing on RNA
counts, binned either per-cell or split into nuclear and cytoplasmic counts
[La ̵+18, SP18, Gor+22, Fu+22]. In doing so, the subcellular spatial factors
that crucially control gene expression [Hei+17, CE22], e.g., geometry-
dependent nuclear export [KH07, RH11], are largely understudied and
unincorporated into transcriptomics analyses.

Beyond the motivation for quantification of subcellular patterns of individual
RNA molecules, we also emphasize the goals and value of mechanistic modeling.
A zoo of sophisticated techniques has arisen for the analysis of (spatial)
transcriptomics data, but these largely phenomenological or statistical
studies [Imb+22, Mah+22] suffer a lack of methodological reproducibility and
interpretability [Gor+22a, CP23]. The proposed solution to these shortcomings
is the use of mechanistic, stochastic models [Lam+20, Gor+22a] as the basis of
inferring gene expression dynamics. Such methods have grown in popularity and
sophistication in recent years [MFN15, Góm+17, Her+17, GP22, Luo+22, Bow+22,
Gup+22, Kil+23a, Vo+23], but remain primarily focused on non-spatial models.
Thus, we have outlined the motivation for subceullar, mechanistic modeling of
RNA molecules as a natural next step in expanding the frontier of spatial
transcriptomics quantifications. We anticipate such techniques will be useful
in advancing understanding of subcellular, spatial factors that modulate gene
expression, including post-transcriptional splicing [DE19, Cot+23, Din+22] and
nuclear export [KH07, RH11].

The basis of this work is a mathematical connection between spatial point
processes and snapshot observations of stochastic birth-death-diffusion
processes. Most crucially, this connection provides statistical information of
stochastic particle positions and counts, while circumventing the need for
costly spatial stochastic simulation [LDS05]. The connection was suggested by
Gardiner in the 1970s [GC77, CG78, Cha+77] and recently reintroduced non-
rigorously in [SGS16] in broad generality of stochastic reaction-diffusion
processes. Here, we provide rigorous particle-perspective justification for
this connection in the setup of birth, death, and diffusion in a domain. The
machinery developed to prove this reveals new insight into generalizations and
limitations of the theory. Inference on spatial point processes is a well-
established topic in the statistics literature [Bad01], but focuses largely on
phenomenological or purely statistical modeling of the underlying intensities,
including in a recent application to the same data considered as our work
[WSV23]. In contrast, the resulting setup here is an intensity derived from
the solution of a birth-diffusion-death PDE. The model problem is similar in
structure to classical PDE-based "inverse problems" [ABT18, Fle+20] but the
crucial difference is that the data presents as collections of stochastic
particle locations rather than the PDE solution itself plus observational
noise. Here we embrace the stochastic particle-by-particle nature of the data
and study the system using spatial point process inference techniques, which
are naturally informed by the solutions of appropriate PDEs.

###  1.2 Overview of work

Figure 1: Motivation and setup of the work. A: Experimental FISH data from
[DE19] demonstrating the basis for modeling within this work. Individual mRNAs
are resolved in both the nucleus and cytoplasm. B: Schematic of the mRNA model
in nuclear domain $\Omega$: birth at rate $\kappa_{\mathrm{+}}$ at gene site
$z$, decay at rate $\kappa_{\mathrm{-}}$, and diffusion with diffusivity $D$,
and instantaneous nuclear export at the boundary $\partial\Omega$ (an
absorbing boundary). C: 1D simulations of the model demonstrating the
heterogeneity in nuclear (or cellular) domain size $L$ and gene site $z$
(black plus).

In this work, we investigate a "toy" model motivated by smFISH imaging data of
subcellular RNA [DE19] that distills salient features of the intrinsically
stochastic [Raj+06, RV08], spatial [Hei+17, CE22], and heterogeneous nature of
the data. In such imaging, such as the example shown in Fig. 1A, individual
mRNA can be resolved in both cytoplasm and the nucleus of cells. Panel B shows
the model considered throughout the remainder of the text for nuclear mRNA in
a domain $\Omega$, with stochastic production at rate $\kappa_{\mathrm{+}}$ at
a gene site $z$, decay at rate $\kappa_{\mathrm{-}}$, and diffusion with
diffusivity $D$. The nuclear boundary is assumed to be instantaneous export at
$\partial\Omega$. To faithfully reflect the cell-to-cell heterogeneity that
crucially affects the observed data, we incorporate different nuclear (or
cellular) domain sizes and gene site positions, shown through a simplified 1D
setup in C. The remainder of this work explores the mathematical nature of the
inverse problem on this collection of heterogeneous snapshots.

In Section 2 we outline the model setup of snapshot data as an observation of
a distribution of particles at time zero in a system in which the initial
condition has been pulled back in time to negative infinity. We show that the
snapshot data can be described as a nonhomogeneous spatial Poisson point
process with an intensity measure that is generally recognized as the
_occupation measure_ of diffusing particles before exiting the domain or being
removed by state-switching [Dur96]. The full proof is provided in Section 4.1
with a nondimensionalization argument found in Section 4.2. This connection is
then translated into a likelihood function that can then be used for inference
of parameters.

In Section 3, we put the theory to practice in the context of our toy model.
In this setting, the likelihood function has a fully explicit form, which we
report and visualize as a function of the number of cell snapshots that are
included in a given experiment. This raises questions concerning the impact on
inference from geometric heterogeneity among the cells and, in Section 3.1, we
calculate the (Fisher) information for experiments run under different
heterogeneity assumptions. In Section 3.2, we explore how these theoretical
results translate into uncertainty in the inference of model parameters. We
report on numerical experiments in which we generated synthetic snapshot data
from simulations of underlying particle processes, and employed Bayesian
methods to quantify joint uncertainty in the inference of particle production
and degradation rates. In Section 3.3 we report on similarly constructed
numerical experiments based on several different cell-heterogeneity
assumptions. We find that the statistical information analysis conducted in
Section 3.2 does faithfully predict Bayesian uncertainty quantifications that
result from simulated data. This encourages continued use of these methods for
both experimental design and practical statistical inference. Finally, in
Section 4.4 we use the information-theoretic framework to demonstrate the
benefit of using explicit particle locations for inference when compared to
spatial binning. Additional numerical implementation details regarding
stochastic process simulations and MCMC sampling for Bayesian inference are
outlined in Appendix A.

##  2 Mathematical and Statistical Methods

###  2.1 Notation

Throughout this study, we have it in mind that biological cells can be
heterogeneous in their geometry, but have similar internal reaction rates and
transport dynamics. Let $M\in\mathbb{N}$ be the number of cells in a study,
and for each $m\in\\{1,2,\ldots,M\\}$, let the compact set
$\Omega_{m}\subset\mathbb{R}^{d}$ denote the domain of the $m$th cell. Within
each cell there are diffusive particles that emerge at rate
$\kappa_{\mathrm{+}}$ at locations that are distributed according to a source
probability distribution $\varphi_{m}:\Omega_{m}\to\mathbb{R}_{+}$. They then
move according to the laws of Brownian motion with uniform diffusivity $D$,
possibly with drift $\alpha:\Omega\to\mathbb{R}^{d}$. If a particle in the
$m$th cell exits the domain $\Omega_{m}$, we assume that it cannot return.
Furthermore, at rate $\kappa_{\mathrm{-}}\geq 0$, particles undergo a change
of state (either they degrade, or engage in a chemical reaction, for example).
We will call a particle _active_ up until the time it either switches state or
exits the domain, and _inactive_ thereafter.

In each cell, at time $t=0$ we observe the locations of the active particles,
denoted
$\boldsymbol{x}_{m}=\\{x_{1m},x_{2m}\ldots,x_{n(\boldsymbol{x}_{m}),m}\\}$,
where $x_{im}\in\Omega_{m}$. We write $n(\boldsymbol{x})$ for the number of
particles in the collection $\boldsymbol{x}$, and emphasize that the number of
particles in each of the $M$ distinct cells $\Omega_{m}$ will be independent
of each other and random. In the foregoing discussion, we emphasize that the
boundary condition, the diffusivity constant, and the emergence and switch
rates are universal among all cells, but the initial location distribution
(often realistically in the cell nucleus) can be different in each different
cell.

As we will see later, the three parameters $D$, $\kappa_{\mathrm{+}}$ and
$\kappa_{\mathrm{-}}$ are not mutually identifiable. We therefore introduce
the ratios

$\lambda\coloneqq\frac{\kappa_{\mathrm{+}}L_{0}^{2}}{D},\text{ and
}\mu\coloneqq\frac{\kappa_{\mathrm{-}}L_{0}^{2}}{D}.$ (1)

That is, we are nondimensionalizing with respect to the particle diffusivity
and the typical cell size $L_{0}$ over the population of observations.

###  2.2 Stochastic model for snapshot observations within a single cell

We now formally introduce the stochastic model that provides the basis for our
inference protocol. For the remainder of this section, we will suppress
dependence on the cell index $m$ and describe the dynamics within a single
focal cell.

In this work, we assume that the emergence (or birth) times of particles are
given by a stationary Poisson point process on the negative half of the real
line. Let $0>T_{1}>T_{2}>\dots$ be an enumeration of points in the negative
half-line, starting with the point closest to zero and decreasing from there.
The initial locations of particles will be given by the sequence of
independent and identically distributed (iid) random variables
$\xi_{i}\overset{iid}{\sim}\varphi$. Let $\\{\mathring{B}_{i}(t)\\}_{t\geq 0}$
be a sequence of iid standard $d$-dimensional Brownian motions (with
$\mathring{B}_{i}(0)=0$), and let $B_{i}(t)\coloneqq\mathring{B}_{i}(t-T_{i})$
be the Brownian motions shifted so that their initial time is the associated
particle emergence time $T_{i}$.

We use the term “time” in the preceding description, but as we show in Section
4.2,we can assume that the dynamics have already been rescaled to
nondimensional coordinates that result in a process that has diffusivity one,
switch rate $\mu$ and birth rate $\lambda$. In this nondimensional framework,
let a drift vector $\alpha\colon\mathbb{R}^{d}\to\mathbb{R}^{d}$ be given, and
define $\\{X_{i}(t)\\}_{t>T_{i}}$ to be the solution of the stochastic
differential equation (SDE)

$\displaystyle\mathrm{d}X_{i}(t)$
$\displaystyle=\alpha(X_{i}(t))\mathrm{d}t+\sqrt{2}\,\mathrm{d}B(t),\quad
t>T_{i};$ (2) $\displaystyle X_{i}(T_{i})$ $\displaystyle=\xi_{i}.$

The partial differential operator associated with this diffusion is defined in
terms of its action on twice-differentiable functions
$f\colon\mathbb{R}^{d}\to\mathbb{R}$,

$\mathcal{L}_{X}\coloneqq\alpha(x)\cdot\nabla f(x)+\Delta f(x)$ (3)

where $\Delta$ is the Laplacian.

Associated with each path is a continuous-time Markov chain (CTMC)
$\\{J_{i}(t)\\}_{t>T_{i}}$ that records the state of the particle at time $t$.
We assume that all particles start in the _active_ state $J_{i}(0)=1$. At rate
$\mu$, the particles switch to the inactive state 0. We label the switch times
$\tau_{\mathrm{switch},i}$ and, due to the Markov assumption, the switch times
are exponentially distributed with densities $p_{i}(t)=\mu
e^{-\mu(t-T_{i})}1_{[T_{i},\infty)}(t)$. Because particles are assumed to not
be able to return to the interior of a cell once they exit, we define the
stopping times $\tau_{\mathrm{exit},i}\coloneqq\inf\\{t\geq
T_{i}\,:\,X_{i}(t)\notin\Omega\\}$.

Figure 2: Visualization of the stochastic model. Particle populations observed
at the snapshot instant ($t=0$) have long histories that cannot be observed.
Here, particles are “born” at a rate $\lambda$ (yellow circles) in the source
location $z$ (yellow line). The particles diffuse until they are either
degraded (green circles) or exit the domain (pink circles). The faded
trajectories are for particles that left the system before time zero, while
the purple trajectories are those that remain active and are observed in the
$t=0$ snapshot. Nondimensional $x,t$ shown with rates $\lambda=10$, $\mu=1.5$
corresponding original simulation parameters
$\kappa_{\mathrm{+}}=5,\kappa_{\mathrm{-}}=0.75,D=0.5,L=3,$ $L_{0}=1,z=1$.

Our model for the “snapshot” data is to define a spatial point process
$\boldsymbol{N}\colon\mathcal{B}(\Omega)\to\mathbb{Z}$ that counts the number
of _active_ particles in subsets of the domain at time 0. To be precise, for
every Borel subset $A\in\mathcal{B}(\Omega)$, we define

$N(A)\coloneqq\sum_{i=1}^{\infty}\delta_{A}\big{(}(X_{i}(0\wedge\tau_{\mathrm{exit},i})\big{)}1_{\\{J(0\wedge\tau_{\mathrm{exit},i})=1\\}},$
(4)

where $\delta_{A}(x)=\int_{A}\delta(x-y)\mathrm{d}y$ and $t\wedge
s\coloneqq\min(t,s)$. In other words, $N(A)$ counts all particles that have
neither exited the system through the boundary $\partial\Omega$, nor switched
from active to inactive, and are located in the set $A$ at time $t=0$. A
schematic of this can observation model can be seen in Fig. 2.

We can now state the main mathematical result, which identifies the
relationship between collections of particle locations in steady state and the
solution of a well-known boundary value problem (BVP).

######  Theorem 2.1.

Let ${(X_{i}(t),J(t))}_{i\in\mathbb{N}}$ be the particle/state pairs defined
by (2) and the foregoing discussion. Let $\Omega\subset\mathbb{R}^{d}$ be an
open domain with compact closure that has a boundary with all regular points
with respect to the SDE (2). Then the associated spatial point process
$\boldsymbol{N}$, defined by (4), giving the locations of _active_ particles
at time zero, is a Poisson spatial process with intensity measure $u$
satisfying the BVP

$\displaystyle\mathcal{L}^{*}u(x)-\mu u(x)$
$\displaystyle=-\lambda\varphi(x),$ $\displaystyle x\in\Omega,$ (5)
$\displaystyle u(x)$ $\displaystyle=0,$ $\displaystyle x\in\partial\Omega.$

where $\mathcal{L}^{*}u=-\nabla\cdot(\alpha(x)u(x))+\Delta u(x)$ is the
adjoint of the operator (3).

######  Remark 2.2.

This is to say, $\boldsymbol{N}$ satisfies the condition that, given any
disjoint sets $A_{1},\ldots A_{K}\subset\mathcal{B}(\Omega)$, the random
variables $\\{N(A_{k})\\}_{k=1}^{K}$ are independent with respective
distributions

$N(A_{k})\sim\mathrm{Pois}\Big{(}\int_{A_{k}}\\!\\!u(x)\mathrm{d}x\Big{)}.$
(6)

The proof of this theorem is the content of Section 4.1.

###  2.3 Framework for inference

The Poisson spatial point process representation of the particle locations
organizes the data into a form that readily allows for statistical inference.
Suppose that we partition the domain $\Omega$ into a fixed collection of sets.
We then “bin the data,” meaning we record the particle locations in terms of
the number of particles in each bin. By Theorem 2.1, the number of particles
in each bin is Poisson distributed and independent of number of particles in
other bins. This yields a likelihood function that can be used for parameter
inference, but depends explicitly on the choice of bins. For example, suppose
that the sets $(A_{1},A_{2},\ldots,A_{K})$ form a partition of cell domain,
meaning that they are mutually disjoint and their union is $\Omega$. For a
given set of particle locations $\boldsymbol{x}$, define
$n_{k}(\boldsymbol{x})$ to be the number of particles observed in the set
$A_{k}$. Then the likelihood function for the parameter vector
$\theta=(\lambda,\mu)$ given the binned data is

$\displaystyle P_{\theta}\big{(}N(A_{1})=n_{1}(\boldsymbol{x}),$
$\displaystyle\ldots,N(A_{K})=n_{K}(\boldsymbol{x}))$ (7)
$\displaystyle=\prod_{k=1}^{K}e^{-\int_{A_{k}}\\!\\!u(x;\,\theta)\mathrm{d}x}\Big{(}\int_{A_{k}}\\!\\!u(x;\,\theta)\mathrm{d}x\Big{)}^{n_{k}(\boldsymbol{x})}\frac{1}{n_{k}(\boldsymbol{x})!}.$

We introduce the notation $u(x;\,\theta)$ to emphasize the dependence of the
intensity measure on $\theta$ through the BVP (5).

Now, if we change the partition of the domain, the associated likelihood
function for the new binning of the data will be different, leading to
different estimates and uncertainty regions for the parameters. Naturally
among all choices for which partition to use, we would prefer to use one that
yields the “best” statistical inference in some rigorous sense. One standard
approach is to consider the likelihood that emerges in a limit in which the
partition mesh size (often defined to be the size of the largest bin) is taken
to zero. The form of the limiting likelihood is well known for spatial Poisson
processes [Bad01, DV+03] and can be expressed directly in terms of the
particle locations, rather than a summary from binning.

######  Definition 2.3 (Likelihood function).

Let $\Omega\subset\mathbb{R}^{d}$ satisfy the conditions of Theorem 2.1. For
each $\boldsymbol{\theta}=(\lambda,\mu)\in\Theta$, let
$u(\cdot;\,\boldsymbol{\theta})$ satisfy the BVP (5) for the given parameter
vector. Then for any set of particle locations
$\boldsymbol{x}=\big{\\{}x_{1},x_{2},\ldots,x_{n(\boldsymbol{x})}\big{\\}}\subset\Omega$,
we define the likelihood function
$L(\cdot;\,\boldsymbol{x})\colon\Theta\to\mathbb{R}_{+}$ as follows:

$L(\boldsymbol{\theta};\,\boldsymbol{x})\coloneqq\left(\prod_{i=1}^{n(\boldsymbol{x})}u(x_{i};\,\boldsymbol{\theta})\right)e^{-\int_{\Omega}u(x;\,\boldsymbol{\theta})\mathrm{d}x}.$
(8)

While the use of this form of likelihood function is standard, in Section 4.3
we include a derivation using the Lebesgue Differentiation Theorem in the
context of Bayesian methods that we have not explicitly seen in the
literature. Moreover, in Section 4.4, we use the notion of statistical
information to quantify the claim that the likelihood function defined by (8)
is “better” than any likelihood function of the form (7) that arises from
binning the particle locations.

We note that the above discussion contains a mix of Bayesian and frequentist
ideas and we will continue to draw on tools from each school throughout this
work. Generally speaking, when we are considering a specific data set and wish
to establish point estimates for the parameters and to quantify uncertainty,
we will use a Bayesian method. This is because the notion of a highest-density
posterior region (HDPR) allows for clean visualization of 2d “joint
uncertainty” for the parameters $\lambda$ and $\mu$, which can then be
summarized with a single quantity: its area. Meanwhile, when we wish to make
comparisons that are not specific to a given data set, we will use the
determinant of the information matrix $\det(I)$. Optimization of this
determinant, or D-optimality [ADT07], is a popular choice in experimental
design (including in other studies of gene expression [FM19]) due to its
ability to capture correlations between inferred parameters in the (inverse)
area of confidence ellipsoids. However, confidence is a frequentist notion,
and we use $\det(I)$ as an analytically tractable proxy for the inverse of the
area of the joint uncertainty region in Bayesian inference. The relationship
between maximizing $\det(I)$ and minimizing expected HDPR is not precise for
finite data, in part because the role of the prior must be considered in the
Bayesian setting [GB15, Sha+22]. But, as we show in Section 3.3, $\det(I)$
does make useful predictions in the present setting.

##  3 Inferring dynamics of the 1D RNA model

Using the framework established in previous sections, we now explore their
consequences and implementation in the context of the 1D model approximating
the key features seen in subcellular spatial RNA transcriptomics data. The
principal result was the connection between the PDE (5) as the intensity for
the corresponding point process. For a single cell, this PDE becomes

$0=\Delta u(x)+\lambda\delta(x-z)-\mu u(x)\,\text{ for }x\in[0,L],\quad
u(0)=0,\,u(L)=0,$ (9)

where $(\lambda,\mu)$ are nondimensional (1) dynamic quantities sought to be
inferred, and $L$ is relative to the typical length seen in the population of
images.

The solution to (9) can be explicitly computed:

$u(x)=\frac{\lambda}{\sqrt{\mu}}\operatorname{csch}\left(\sqrt{\mu}L\right)\sinh\left(\sqrt{\mu}\min\\{x,z\\}\right)\sinh\left(\sqrt{\mu}(L-\max\\{x,z\\})\right).$
(10)

The expected number of particles is then, from (6),

$\mathbb{E}[N]=\int_{0}^{L}u(x)\,\mathrm{d}x=\frac{\lambda}{\mu}\left(1-\text{sech}\left(\frac{\sqrt{\mu}L}{2}\right)\cosh\left(\frac{1}{2}\sqrt{\mu}(L-2z)\right)\right)$
(11)

The expectation in (11) shows an intuitive structure of the product of the
pure birth-death mean $\lambda/\mu$, and a geometry-dependent term (always
strictly less than 1) that scales the former.

###  3.1 Likelihood and information from analytical solution

Equipped with the statistical characterization from Theorem 2.1 and the
analytical expression of the intensity for a single snapshot in (10), we can
now examine the behavior and emergent lessons arising from this inference
setup.

Figure 3: Behavior of the likelihood. A: $\lambda-\mu$ evaluations of the
likelihood $L(\boldsymbol{\theta};\,\boldsymbol{x})$ for increasing numbers of
snapshots $M$ with fixed $\lambda=500,\mu=10,L=1,z=0.5$. At low $M$, the
likelihood peaks around the line of constant $\mathbb{E}[N]$, and as $M$
increases, sharply peaks around the true parameter values. B/C: Slices of the
likelihood shown in panel A evaluated at $\mu_{\text{true}}$ and
$\lambda_{\text{true}}$ respectively.

We first explore the behavior of the likelihood (8) as a function of the
system parameters and quantity of data. To do so, we perform stochastic
simulations of the underlying birth-death-diffusion trajectories and then
terminate them at an experimental time $T$. The particle positions at this
observation time are the data used for inference. Further details on the
stochastic simulation can be found in Appendix A. Beginning with the simplest
setup, we fix $\lambda=500,\mu=10$, $L=1$ and $z=1/2$ and explore how
increasing the number of snapshots $M$ shapes the likelihood, the results of
which can be seen in Fig. 3A. As $M$ increases, the likelihood first
concentrates around a line and then within this line, begins to peak at the
true parameter values. To understand this behavior, we plot the $\lambda-\mu$
line of constant $\mathbb{E}[N]$ from (11) superimposed on the likelihood.
From this, we see that the likelihood concentrates around the line defining
the true number of particles, and then further resolves within this line from
their spatial information. In panels B and C of the same figure, the same
information is shown as slices of the likelihood with one parameter at a time
held to be the true value. From these panels, we see again that sufficiently
large $M$ (here, $M\approx 1000$) yields likelihoods that are tightly
concentrated around the true values used to generate the synthetic data. In
the next section, we will explore more concretely how the spread of the
likelihood translates into uncertainty in the inference itself.

Before performing inference on the synthetic data, other useful quantities can
be computed directly from the likelihood. For a likelihood
$L(x;\boldsymbol{\theta})$, the corresponding information matrix [VHO20] is
defined by the elements

$I_{ij}\coloneqq\mathbb{E}_{\boldsymbol{\theta}}\left[(\partial_{\theta_{i}}\log
L(\boldsymbol{\theta};\,\boldsymbol{x}))(\partial_{\theta_{j}}\log
L(\boldsymbol{\theta};\,\boldsymbol{x}))\right]=-\mathbb{E}_{\theta}\left[\partial_{\theta_{i}}\partial_{\theta_{j}}L(\boldsymbol{\theta};\,\boldsymbol{x})\right],$
(12)

where the latter equality holds under regularity conditions straightforward to
verify for Poisson point processes [Cla22]. Our parameters of interest are
$\boldsymbol{\theta}=(\lambda,\mu)$ so the information matrix takes the form

$I=-\mathbb{E}_{\theta}\begin{bmatrix}\partial_{\lambda\lambda}L(x;\lambda,\mu)&\partial_{\lambda\mu}L(x;\lambda,\mu)\\\
\partial_{\lambda\mu}L(x;\lambda,\mu)&\partial_{\mu\mu}L(x;\lambda,\mu)\end{bmatrix}.$
(13)

Importantly, the expectation is over both stochastic particle positions and
number. To do so, we employ machinery from the theory of spatial point
processes, specifically the Campbell-Hardy formula stated in Theorem 3.1.

######  Theorem 3.1 (Campbell-Hardy, [Bad01]).

Let $N$ be governed by a general point process on domain $x\in\Omega$ with
intensity $u(x)$ so that $\mathbb{E}[N(B)]=\int_{B}u(x)\mathrm{d}x$.Let $f(x)$
be any measurable function. Then, the expectation of the sum over observations
of the point process is

$\mathbb{E}\left[\sum_{i=1}^{N}f(x_{i})\right]=\int_{\Omega}f(x)u(x)\mathrm{d}x.$
(14)

Through Theorem 3.1, the elements of our information matrix in (13) can be
related to the PDE solution $u$ by (10). To do so, the log-likelihood
$\ell(\boldsymbol{\theta};\boldsymbol{x})$ from (8) is

$\ell(\boldsymbol{\theta};\boldsymbol{x})\coloneqq\log
L(\boldsymbol{\theta};\,\boldsymbol{x})=-\int_{\Omega}u(x;\boldsymbol{\theta})\mathrm{d}x+\sum_{i=1}^{N}\log
u(x_{i};\boldsymbol{\theta})$

and
$\partial_{\theta_{i}\theta_{j}}\ell=-\int_{\Omega}\partial_{\theta_{i}\theta_{j}}u(x;\boldsymbol{\theta})\mathrm{d}x+\sum_{i=1}^{N}\partial_{\theta_{i}\theta_{j}}\log
u(x_{i};\boldsymbol{\theta})$.

$\displaystyle
I_{ij}(\boldsymbol{\theta})=-\mathbb{E}[\partial_{\theta_{i}\theta_{j}}\ell(\boldsymbol{\theta};\boldsymbol{x})]$
$\displaystyle=-\mathbb{E}\left[-\int_{\Omega}\partial_{\theta_{i}\theta_{j}}u\mathrm{d}x\right]-\mathbb{E}\left[\sum_{i=1}^{N}\partial_{\theta_{i}\theta_{j}}\log
u(x_{i};\boldsymbol{\theta})\right]$ (15)
$\displaystyle=\int_{\Omega}\partial_{\theta_{i}\theta_{j}}u(x;\boldsymbol{\theta})\mathrm{d}x-\mathbb{E}\left[\sum_{i=1}^{N}\partial_{\theta_{i}\theta_{j}}\log
u(x_{i};\boldsymbol{\theta})\right]$
$\displaystyle=\int_{\Omega}\partial_{\theta_{i}\theta_{j}}u(x;\boldsymbol{\theta})\mathrm{d}x-\int_{\Omega}(\partial_{\theta_{i}\theta_{j}}\log
u(x;\boldsymbol{\theta})u(x;\boldsymbol{\theta})\,\mathrm{d}x.$

The first term is constant, and therefore is its own expected value. The
second, stochastic sum, takes the form of (14) with
$f(x)=\partial_{\theta_{i}\theta_{j}}\log u(x)$ and is evaluated with Theorem
3.1. The integrands of (15) can be manipulated to a more tractable form by
expanding derivatives

$\displaystyle I_{ij}(\boldsymbol{\theta})$
$\displaystyle=-\int_{\Omega}(\partial_{\theta_{i}\theta_{j}}\log
u(x;\boldsymbol{\theta}))u(x;\boldsymbol{\theta})\,\mathrm{d}x+\int_{\Omega}\partial_{\theta_{i}\theta_{j}}u(x;\boldsymbol{\theta})\,\mathrm{d}x$
(16) $\displaystyle=\int_{\Omega}(\partial_{\theta_{i}}\log
u)(\partial_{\theta_{j}}\log
u(x;\boldsymbol{\theta}))u(x;\boldsymbol{\theta})\,\mathrm{d}x.$
$\displaystyle=\int_{\Omega}\frac{(\partial_{\theta_{i}}u(x;\boldsymbol{\theta}))(\partial_{\theta_{j}}u(x;\boldsymbol{\theta}))}{u(x;\boldsymbol{\theta})}\,\mathrm{d}x.$

If $\lambda>0$, then $u>0$ and therefore $I_{ij}>0$ so identifiability
[Wie+21, Bro+20] of both parameters can be concluded. Moreover, the analytical
solution to the PDE allows us to compute the information explicitly. For
instance,

$I_{11}(\boldsymbol{\theta})=-\mathbb{E}_{\theta}\left[\partial_{\lambda\lambda}\ell(\boldsymbol{\theta};\boldsymbol{x})\right]=\left(\lambda\mu\right)^{-1}\left[1-\operatorname{sech}\left(\sqrt{\mu}L/2\right)\cosh\left(\sqrt{\mu}(L-2z)/2\right)\right].$
(17)

The remaining entries can be computed straightforwardly using Mathematica but
have unwieldy expressions not shown here.

Figure 4: Behavior of the information. A: Determinant of the information $\det
I$ in the $z-\mu$ plane with $\lambda=500$, $L=1$ fixed. B: same information
as panel B, with $\det I$ as a function of $z$ shown for various values of
$\mu$, with the maximizing value of $z$ labeled. B inset: the-information
maximizing source location $z$ as a function of $\mu$. C: Determinant of the
information $\det I$ in the $L-\mu$ plane with $\lambda=500$, $z=L/2$ fixed.
D: same information as panel C, with $\det I$ as a function of $L$ shown for
various values of $\mu$.

In Fig. 4, we employ the analytical expressions for the information matrix to
explore the behavior of $\det I$ as a function of the underlying parameters.
We interpret this quantity’s importance through the classical lens of the
Cramér-Rao bound [Cla22]: any unbiased estimator of the parameter’s variance
is bounded below by the inverse of information. Informally, we interpret
larger information setups as "easier" inverse problems. We later verify this
interpretation through posterior widths on actual inference setups.

In Fig. 4A and B, we fix $\lambda=500,L=1$ and vary both source location $z$
and degradation rate $\mu$. For large values of $\mu$, the behavior is
intuitive: the information is maximized when the source location is at the
center, as this minimizes the overall number of particles exiting the boundary
and corresponding information loss. For smaller values of $\mu$, a bifurcation
in the behavior occurs and the optimal source location becomes non-trivial:
occurring at two symmetric locations somewhere between the cell center and
boundary. This non-trivial behavior arises only in the simultaneous inference
of parameters, suggesting its origin. We interpret this non-trivial behavior
as arising from the trade-off between seeing more particles (sources closer to
the center) or seeing a wider range of the spatial gradient (for off-center
source locations).

Panels C and D investigate the cell size’s impact on the information by fixing
$z/L$ and varying $L$ and $\mu$. Intuitively, larger cells (relative to the
typical size) yield more information due to less information loss through the
absorbing boundary. Interestingly, this effect saturates for sufficiently
large values of $\mu$, when the typical lengthscale of a trajectory is far
shorter than the domain size, so no additional information is gained through
enlarging the domain.

###  3.2 Bayesian inference on synthetic data

For its noted strengths in robustness and natural uncertainty quantification
that we defer to arguments elsewhere [YC19], we will deploy Bayesian
inference. We must therefore specify prior distributions [Gel+95] for
$\lambda$ and $\mu$. For the production rate $\lambda$, we take

$\lambda\sim\mathrm{Gamma}(\alpha,\beta),\qquad
p(\lambda)=\frac{\lambda^{\alpha-1}e^{-\beta\lambda}\beta^{\alpha}}{\Gamma(\alpha)}.$
(18)

Assuming the total prior is independent for each parameter,
$p(\boldsymbol{\theta})=p(\lambda)p(\mu)$, then for a single observation, we
have

$\displaystyle p(\boldsymbol{\theta}|x)$ $\displaystyle\propto
p(\boldsymbol{\theta})L(\boldsymbol{\theta};\boldsymbol{x})=p(\mu)p(\lambda)e^{-\lambda\int_{\Omega}v(x;\mu)\mathrm{d}x}\prod_{i=1}^{N}\lambda
v(x_{i};\mu)$ (19)
$\displaystyle=p(\mu)\prod_{i=1}^{N}v(x_{i};\mu)\frac{\lambda^{\alpha-1}e^{-\beta\lambda}\beta^{\alpha}}{\Gamma(\alpha)}\lambda^{n}e^{-\lambda\int_{\Omega}v(x;\mu)\mathrm{d}x}$
$\displaystyle=p(\mu)\prod_{i=1}^{N}v(x_{i};\mu)\frac{\lambda^{(\alpha+n)-1}e^{-\lambda\left(\beta+\int_{\Omega}v(x;\mu)\mathrm{d}x\right)}\beta^{\alpha}}{\Gamma(\alpha)}.$

Here, we have introduced the definition of $v$ as $u(x;\mu,\lambda)=\lambda
v(x;\mu)$, noting that (10) is linear in $\lambda$ and following standard
practice of identifying a "scale" parameter in the point process literature
[Bad01]. This yields the convenient conditional posterior for $\lambda$
arising from conjugacy

$p(\lambda|\boldsymbol{x},\mu)=\mathrm{Gamma}\Big{(}\alpha+n(\boldsymbol{x}),\beta+\int_{\Omega}v(x;\mu)\,\mathrm{d}x\Big{)}.$
(20)

For $M$ observations, each with a different domain, this generalizes to

$p(\lambda|\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{m},\mu)=\mathrm{Gamma}\left(\alpha+\sum_{m=1}^{M}n(\boldsymbol{x}_{m}),\beta+\sum_{m=1}^{M}\int_{\Omega_{m}}v(x;\mu)\,\mathrm{d}x\right).$
(21)

For $\mu$, we choose a uniform prior
$\mu\sim\mathrm{unif}(0,\mu_{\mathrm{max}})$ with large $\mu_{\mathrm{max}}$
for uninformedness. With these choices, of priors, it is straightforward to
implement a Gibbs-in-Metropolis sampler for the posterior distribution, with a
Gibbs step using the conjugacy for $\lambda$ and Metropolis step for $\mu$.
Further details on the MCMC implementation (including verification of
convergence) can be found in Appendix A.

Figure 5: Bayesian inference on synthetic data. A: Shaded regions show 89%
highest-density posterior region (HDPR) for $M=10,50,500$ snapshots and fixed
parameters $\lambda=500,\mu=10,z=1/2,L=1$. As $M$ increases, the posterior
regions concentrate around the true parameters. B,C: Marginal posterior
samples for $\lambda$ and $\mu$ respectively, along with the initial prior
distributions for each (with prior means shown with dots on the curves). A
inset: area of the HDP region as a function of $M$, showing $\sim M^{-1}$
scaling.

The results of performing Bayesian inference on synthetic data can be seen in
Fig. 5. The data considered here is the same as that in Fig. 3, with fixed
parameters and varied number of snapshots. In panel A, we depict the 89% HDPR
(highest density posterior region) [Gel+95], computed by evaluating
$L(\boldsymbol{\theta};\,\boldsymbol{x})p(\boldsymbol{\theta})$ for each
posterior sample $\theta$, retaining the top 89%, and then computing the
convex hull of these samples. The resulting uncertainty qualitatively mimics
the lesson seen from the likelihood. For small values of $M$, the posterior
concentrates around the parameter region that produces an expected number of
particles equal to the observed count, and then resolves to be sharply peaked
around the true parameters with data from more snapshots. Qualitatively, the
same result can be seen in the marginal posterior samples seen in Fig. 5B and
C. To more explicitly quantify the dependence of uncertainty on data, we plot
the HDPR area as a function of the number of snapshots $M$ in the inset of
Fig. 5A, and find a convincing $\sim M^{-1}$ scaling. This scaling can be
interpreted as arising from the product (an area) of the uncertainty of each
parameter $\sim M^{-1/2}$ arising from the central limit theorem.

###  3.3 Cell-to-cell heterogeneities and their impact on inference

In the previous figure, we assumed all cells were identical in shape and
source location. However, in experimental setups, there is enormous cell-to-
cell heterogeneity. The likelihood results shown in Fig. 3 suggest that this
heterogeneity should affect inference quality, as both domain size and source
location both modulate the total information. However, those quantities were
generated through the computation of the information for a single observation.
How do these propagate to a heterogeneous population? To investigate this
heterogeneity, we extend the information matrix calculation to now assume that
the domain parameters $\boldsymbol{\vartheta}=(z,L)$ follow some probability
distribution. The expectation can be extended straightforwardly by the Law of
Total Expectation:

$I_{ij}(\boldsymbol{\theta})=-\mathbb{E}_{\boldsymbol{\theta},\boldsymbol{\vartheta}}[\partial_{\theta_{i}}\partial_{\theta_{j}}\log
L(x;\boldsymbol{\theta},\boldsymbol{\vartheta})]=-\mathbb{E}_{\boldsymbol{\vartheta}}\mathbb{E}_{\boldsymbol{\theta}}[\partial_{\theta_{i}}\partial_{\theta_{j}}\log
L(x;\boldsymbol{\theta}|\boldsymbol{\vartheta})].$ (22)

Ultimately, this just reduces to integrating the results of previous figures
against the probability distributions for $z$ and $L$. To parameterize the
noise, we take $z/L\,|\,L\sim\mathrm{Beta}(\sigma_{z}^{-1},\sigma_{z}^{-1})$
and $L\sim\Gamma(\sigma_{L}^{-2},\sigma_{L}^{2})$. The choice of distribution
for source location was chosen so that
$\mathbb{E}_{\boldsymbol{\vartheta}}[z\,|\,L]=L/2$ and increasing $\sigma_{z}$
spreads out the source symmetrically around this location, with $\sigma_{z}=1$
corresponding to spatially uniform locations. Similarly, $L$ is supported on
the non-negative numbers, and the distribution is chosen so that
$\mathbb{E}_{\boldsymbol{\vartheta}}[L]=1$ and increasing $\sigma_{L}$ spreads
out the distribution of cell sizes in a plausible way.

Figure 6: Impact of cell-to-cell heterogeneities. A: Information with $L=1$,
$\lambda=500$ fixed and source location varied by a spatially symmetric
parameterization with $\sigma_{z}\to 0$ is a fixed source and $\sigma_{z}=1$
is spatially uniform. For low values of $\mu$, spatial heterogeneity increases
information, whereas it decreases the information for large $\mu$. B:
Information with $z=L/2$ and varied $L$ parameterized by a distribution such
that $\mathbb{E}_{\boldsymbol{\vartheta}}[L]=1$ and increasing $\sigma_{L}$
spreads out the distribution. Again, lower $\mu$ scenarios have increased
information from spatial heterogeneity. C: Repeated Bayesian inference on
synthetic data with $\mu=1$, varied source location $\sigma_{z}=0.75$, cell
size $\sigma_{L}=0.25$, and both random. All heterogeneous setups have lower
uncertainty (smaller HDPR area) than static $z,L$.

The impact of cell-to-cell heterogeneity on the information can be seen in
Fig. 6. Panel A the information $\det I$ as a function of $\sigma_{z}$ for
various values of $\mu$ and fixed $\lambda=500,L=1$. For low values of $\mu$,
varying source locations increase the amount of information in observations,
whereas for large $\mu$, they lower the information. This is in alignment with
the results seen in Fig. 4 for fixed source locations, where low values of
$\mu$ corresponded to non-trivial optimal source locations, and large $\mu$
have optimal source locations in the center.

A similar effect can be seen for varying $L$ in panel B of Fig. 6: for low
values of $\mu$, varied cell sizes increase the information, whereas for large
$\mu$, they decrease it. This result is somewhat more surprising. For
increasing values of $\sigma_{L}$, the median $L$ decreases, meaning that more
mass of the distribution is shifted toward cell sizes smaller than the typical
size. The fixed-size results say that smaller cells encode less information.
However, taken together, the population-heterogeneous results suggest that the
information from large cells offsets the lack thereof from small, and creates
a net gain in information.

An outstanding curiosity remains of whether these information results are
meaningful for actual inference. To investigate this, we repeat the Bayesian
inference pipeline from the previous section for various setups of
heterogeneities: varied $z$, varied $L$, and both varied. For a fixed
experiment of $M=100$ images, synthetic data (with $\mu$ chosen in the regime
of increase information) is generated and inference is performed. The HDPR
area is computed as a measure of the overall uncertainty in the estimation.
This inference experiment is then repeated 1000 times, and the results can be
seen in Fig. 6 panel C. The results from estimation support the previous:
heterogeneity does indeed increase the amount of information in the snapshots,
and lowers uncertainty in estimation. All three scenarios with random source
location, cell size, or both random yield smaller HDPR areas than the static
scenario.

##  4 Support for the inference framework.

###  4.1 Proof of Theorem 2.1

In order to prove that the locations of particles is a Poisson spatial
process, we need only show that for any finite collection of disjoint sets
$A_{1},\ldots,A_{k}$, the number of particles in the sets are independent and
Poisson distributed. If we find a function $u$ such that
$\int_{A}u(x)\mathrm{d}x=P(A)$ for all $A$, then this the intensity measure.

Poisson distribution. To fix notation, let $\rho_{A}(t)$ be the probability
that a particle that is distributed according to $\varphi$ at time 0 and is in
the set $A$ at time $t>0$. Recall, however, that in the model the particles
emerge at a rate $\lambda$ according to a Poisson arrival process for some
interval of time $(-h,0)$ and we are counting the particles still in the
system at time zero. Since this means that the initial times are iid and
uniformly distributed over $(-h,0)$, we can define

$\rho^{(h)}_{A}=\frac{1}{h}\int_{-h}^{0}\rho_{A}(-s)\mathrm{d}s=\frac{1}{h}\int_{0}^{h}\rho_{A}(t)\mathrm{d}t.$
(23)

By Poisson thinning, if the total number of particles that emerge is
$\text{Pois}(\lambda h)$, then the distribution of particles in $A$ at time
zero is $\text{Pois}(\lambda h\rho_{A}^{(h)})$. Taking the limit as
$h\to\infty$, we have that the number of particles in $A$ has a Poisson
distribution with mean

$\lambda\rho_{A}=\lambda\int_{0}^{\infty}\rho_{A}(t)\mathrm{d}t.$ (24)

Intensity measure. The time-dependent PDE that describes the probability
density of particles that have not yet been removed by state-switching or by
being absorbed at the boundary of a domain is [CM77, p. 230]

$\displaystyle\partial_{t}\rho(x,t)$
$\displaystyle=L^{*}\rho(x,t)-\mu\rho(x,t),$ $\displaystyle x\in\Omega,t>0,$
(25) $\displaystyle\rho(x,t)$ $\displaystyle=0,$ $\displaystyle
x\in\partial\Omega,t>0,$ $\displaystyle\rho(x,0)$ $\displaystyle=\varphi(x),$
$\displaystyle x\in\Omega.$

Defining $\rho(x)=\int_{0}^{\infty}\rho(x,t)\mathrm{d}x$, [Gar04, p. 163-167]
we integrate (25) to find

$\displaystyle\big{(}\lim_{t\to\infty}\rho(x,t)\big{)}-\rho(x,0)$
$\displaystyle=L^{*}\rho(x)-\mu\rho(x),$ $\displaystyle x\in\Omega,$ (26)
$\displaystyle\rho(x)$ $\displaystyle=0,$ $\displaystyle x\in\partial\Omega.$

Since particles exit the domain in finite time with probability one,
$\rho(x,t)\to 0$ as $t\to\infty$. Imposing the initial condition yields

$\displaystyle L^{*}\rho(x)-\mu\rho(x)$ $\displaystyle=-\varphi(x),$
$\displaystyle x\in\Omega$ (27) $\displaystyle\rho(x)$ $\displaystyle=0,$
$\displaystyle x\in\partial\Omega.$

In order to achieve the desired mean for the Poisson distribution in (24), we
define $u(x)=\lambda\rho(x)$ and note that $u$ can be attained by multiplying
the initial distribution function by $\lambda$. This yields the BVP stated in
the theorem and the proof is complete.

$\square$

###  4.2 Nondimensionalization

In this subsection, we outline the nondimensionalization argument justifying
the SDE (2). Suppose that the original SDE satisfies the SDE

$\mathrm{d}X(t)=\alpha(X(t))\mathrm{d}t+\sqrt{2D}\mathrm{d}W(t)$ (28)

with initial condition $X(t_{0})=x_{0}$. The particle can meet one of two
fates: a state-switch due to degradation, at rate $\kappa_{\mathrm{-}}$, or
exit through the boundary of the domain. Define

$\tilde{t}\coloneqq\frac{D}{L_{0}^{2}}t,\,\,\tilde{x}\coloneqq\frac{1}{L_{0}}x,\text{
and }\tilde{\alpha}(\tilde{x}):=\frac{L_{0}}{D}\alpha(x).$

Then we claim that $\tilde{X}(\tilde{t})\coloneqq\frac{1}{L}X(t)$ satisfies

$\mathrm{d}\tilde{X}(\tilde{t})=\tilde{\alpha}(\tilde{X}(t))\mathrm{d}\tilde{t}+\sqrt{2}\mathrm{d}\widetilde{B}(\tilde{t})$
(29)

where $\tilde{B}(\tilde{t})$ is a standard Brownian motion. The argument
follows from the integral form of the SDE:

$\displaystyle\tilde{X}(\tilde{t})-\tilde{X}(\tilde{t}_{0})$
$\displaystyle=\frac{1}{L_{0}}\left(X(t)-X(t_{0})\right)$
$\displaystyle=\int_{t_{0}}^{t}\frac{1}{L_{0}}\alpha(X(s))\mathrm{d}s+\sqrt{\frac{2D}{L_{0}^{2}}}\left(B(t)-B(t_{0})\right).$

For the drift term, under the substitution $\tilde{s}=\frac{D}{L_{0}^{2}}s$ we
have

$\displaystyle\int_{t_{0}}^{t}\frac{1}{L_{0}}\alpha(X(s))\mathrm{d}s$
$\displaystyle=\int_{\frac{Dt_{0}}{L_{0}^{2}}}^{\frac{Dt}{L_{0}^{2}}}\frac{L_{0}}{D}\alpha\left(X\Big{(}\frac{L_{0}^{2}}{D}\tilde{s}\Big{)}\right)\mathrm{d}\tilde{s}$
$\displaystyle=\int_{\tilde{t_{0}}}^{\tilde{t}}\frac{L_{0}}{D}\alpha\left(L_{0}\tilde{X}(\tilde{s})\right)\mathrm{d}\tilde{s}$
$\displaystyle=\int_{\tilde{t}_{0}}^{\tilde{t}}\tilde{\alpha}\left(\tilde{X}(\tilde{s})\right)\mathrm{d}\tilde{s}$

Meanwhile, define $\tilde{B}(\tilde{t})=\sqrt{\frac{D}{L_{0}^{2}}}B(t)$. Then
$E\big{(}\tilde{B}(\tilde{t})\big{)}=0$ and

$\textrm{Var}\big{(}\tilde{B}(\tilde{t})\big{)}=\frac{D}{L_{0}^{2}}\textrm{Var}(B(t))=\frac{D}{L_{0}^{2}}t=\tilde{t}$
(30)

meaning that $\tilde{B}$ is a standard Brownian motion.

In terms of the essential stopping times, we define $\mu=L_{0}^{2}/D$ (as in
(1)), and so $\tilde{\tau}_{\mathrm{switch}}=\tau_{\text{switch}}D/L_{0}^{2}$
is the corresponding nondimensional switch time. Moreover, defining
$\tilde{\Omega}=\\{x\,:\,L_{0}x\in\Omega\\}$ and
$\tilde{\tau}_{\mathrm{exit}}:=\inf\\{\tilde{t}>\tilde{t}_{0}\,:\,\tilde{X}(\tilde{t})\notin\tilde{\Omega}\\}$.
We can readily show path-by-path that
$\tilde{\tau}_{\mathrm{exit}}=\tau_{\text{exit}}D/L_{0}^{2}$. Analogous to the
dimensional case, we define
$\tilde{\tau}\coloneqq\min(\tilde{\tau}_{\mathrm{switch}},\tilde{\tau}_{\mathrm{switch}})$

Altogether, it follows that, for any $A\in\Omega$, if we let
$\tilde{A}=\\{x\,:\,L_{0}x\in\Omega\\}$, then

$P(\tilde{X}(0\wedge\tilde{\tau})\in\tilde{A})=P(X(0\wedge\tau)\in A)).$ (31)

Hence, if we wish to solve for the intensity measure of a given dimensional
system, we can rescale space by $L_{0}$ and the rates as in (1) and use the
intensity measure of the associated nondimensionalized system for statistical
inference.

###  4.3 Derivation of the likelihood function

######  Proposition 4.1.

Let $\\{\mathcal{A}_{K}\\}_{K=1}^{\infty}$ (with
$\mathcal{A}_{K}=(A_{1,K},\ldots,A_{K,K})$) be a sequence of
$\Omega$-partitions that has the following properties:

  1. (i)

the partitions are nested in the sense that if $K_{1}<K_{2}$, then every set
in $\mathcal{A}_{K_{2}}$ is contained within some member of
$\mathcal{A}_{K_{1}}$; and

  2. (ii)

the partition mesh sizes,
$\big{(}\max_{k\in\\{1,\ldots,K\\}}\\{|A_{k,K}|\\}\big{)}$, where $|A|$
denotes the volume of a set $A$, decrease to zero with $K$.

Fix $\boldsymbol{x}=\big{(}x_{1},x_{2},\ldots,x_{n(\boldsymbol{x})}\big{)}$ to
be a single snapshot of particle locations arising from the dynamics defined
for Theorem 2.1 and let $\pi(\boldsymbol{\theta})$ be a prior distribution for
$\boldsymbol{\theta}=(\lambda,\mu)$. Define
$\pi_{K}(\boldsymbol{\theta},\boldsymbol{x})$ be the posterior distribution of
$\boldsymbol{\theta}$ under the prior $\pi$ when the locations are reported
through binning in the sets of the partition $A_{K}$. Then,

$\pi(\theta\,|\,\boldsymbol{x})\coloneqq\lim_{K\to\infty}\pi_{K}(\theta\,|\,\boldsymbol{x})=\frac{L(\theta;\,\boldsymbol{x})\pi(\boldsymbol{\theta})}{\int_{\Theta}L(\vartheta;\,\boldsymbol{x})\pi(\vartheta)\mathrm{d}\vartheta}.$
(32)

where $L(\boldsymbol{\theta};\,\boldsymbol{x})$ is the likelihood function
given in Definition 2.3.

######  Remark 4.2.

A construction of such a sequence of partitions can be achieved, for example,
by starting with any partition of $\Omega$ and creating each refinement by
randomly selecting one of the partition’s sets and splitting in half.

###### Proof.

We can represent the collection of particle locations as a measure
$\delta_{\boldsymbol{x}}(x)\coloneqq\sum_{i=1}^{n(\boldsymbol{x})}\delta_{x_{i}}(x)$
[DV+03] and for each member $A_{k,K}$ of a partition $\mathcal{A}_{K}$, we can
define

$n_{k,K}(\boldsymbol{x})\coloneqq\int_{A_{k,K}}\delta_{\boldsymbol{x}}(x)\mathrm{d}x.$
(33)

We then can write the posterior distribution $\pi_{K}(\cdot|\boldsymbol{x})$
associated with the partition $\mathcal{A}_{K}$ as follows:

$\pi_{K}(\theta\,|\,\boldsymbol{x})=\frac{P_{\theta}\big{(}N(A_{1,K})=n_{1,K}(\boldsymbol{x}),\ldots,N(A_{K,K})=n_{K,K}(\boldsymbol{x})\big{)}\pi(\theta)}{\int_{\Theta}P_{\vartheta}\big{(}N(A_{1,K})=n_{1,K}(\boldsymbol{x}),\ldots,N(A_{K,K})=n_{K,K}(\boldsymbol{x}))\pi(\vartheta)\mathrm{d}\vartheta}$
(34)

Here the subscript of the probability measure denotes that we are evaluating
the probability with the indicated parameter vector. The full parameter space
is denoted $\Theta$, which in our case is
$\mathbb{R}_{+}\times\mathbb{R}_{+}$. Since the probability in the above
formula is expressed in terms of an event comprised of independent Poisson
random variable results, we have the following explicit formula for the
likelihood (suppressing dependence on $K$), which follows from Remark 2.2:

$\displaystyle P_{\theta}\big{(}N(A_{1})=n_{1}(\boldsymbol{x}),$
$\displaystyle\ldots,N(A_{K})=n_{K}(\boldsymbol{x}))$ (35)
$\displaystyle=\prod_{k=1}^{K}e^{-\int_{A_{k}}\\!\\!u(x;\,\theta)\mathrm{d}x}\Big{(}\int_{A_{k}}\\!\\!u(x;\,\theta)\mathrm{d}x\Big{)}^{n_{k}(\boldsymbol{x})}\frac{1}{n_{k}(\boldsymbol{x})!}.$

Now, because the partitions are nested, there exists a $K_{*}$ such that for
all $K>K_{*}$, all partition sets contain at most one particle. If we write
$A_{k_{i}}$ to be the bin containing the particle location $x_{i}$, we have
that $n_{k_{i}}(\boldsymbol{x})=1$ for each $i$, and $n_{k}(\boldsymbol{x})=0$
otherwise. For the purpose of writing the likelihood function, we gather all
bins not containing a particle into one set

$A_{0}:=\Omega\setminus\Big{\\{}\bigcup_{i=1}^{n(\boldsymbol{x})}A_{k_{i}}\Big{\\}}.$

For such a partition, (35) has the form

$\displaystyle
P_{\theta}\big{(}N(A_{1})=n_{1}(\boldsymbol{x}),\ldots,N(A_{K})=n_{K}(\boldsymbol{x}))$
(36)
$\displaystyle\qquad\qquad=P_{\theta}\big{(}N(A_{0})=0,N(A_{k_{1}})=1,\ldots,N(A_{k_{n}(\boldsymbol{x}})=1\big{)}$
$\displaystyle\qquad\qquad=P_{\theta}\big{(}N(A_{0})=0\big{)}\prod_{i=1}^{n(\boldsymbol{x})}P_{\theta}\big{(}N(A_{k_{i}})=1\big{)}$
$\displaystyle\qquad\qquad=e^{-\int_{A_{0}}\\!\\!u(x;\,\theta)\mathrm{d}x}\prod_{i=1}^{n(\boldsymbol{x})}e^{-\int_{A_{k_{i}}}\\!\\!u(x;\,\theta)\mathrm{d}x}\Big{(}\int_{A_{k_{i}}}\\!\\!u(x;\,\theta)\mathrm{d}x\Big{)}$
$\displaystyle\qquad\qquad=e^{-\int_{\Omega}u(x;\,\theta)\mathrm{d}x}\prod_{i=1}^{n(\boldsymbol{x})}\int_{A_{k_{i}}}\\!\\!u(x;\,\theta)\mathrm{d}x.$

Each term in the final product is going to zero, but by the Lebesgue
Differentiation theorem (and observing that for the $i$th particle location,
the sequence of $x_{i}$-containing partition sets $A_{k_{i},K}$ is nested), we
have that

$\lim_{K\to\infty}\frac{1}{|A_{k_{i},K}|}\int_{A_{k_{i}}}\\!\\!u(x;\,\theta)\mathrm{d}x=u(x_{i}).$

It follows that

$\displaystyle\lim_{K\to\infty}\frac{P_{\theta}\big{(}N(A_{\epsilon,1})=n_{1}(\boldsymbol{x}),\ldots,N(A_{\epsilon,K})=n_{K}(\boldsymbol{x})\big{)}}{\prod_{i=1}^{n(\boldsymbol{x})}|A_{k_{i},K}|}=e^{-\int_{\Omega}u(x)\mathrm{d}x}\prod_{i=1}^{n(\boldsymbol{x})}u(x_{i}).$
(37)

The right-hand side is precisely $L(\theta;\,\boldsymbol{x})$, the function
appearing in the likelihood definition (8). Multiplying the numerator and
denominator of (34) by the rescaling factor
$\prod_{i=1}^{n(\boldsymbol{x})}|A_{k_{i},K}|^{-1}$, we obtain the limit of a
sequence of posterior distributions $\pi_{K}(\cdot\,|\,\boldsymbol{x})$, which
can be written

$\pi(\theta\,|\,\boldsymbol{x})\coloneqq\lim_{K\to\infty}\pi_{K}(\theta\,|\,\boldsymbol{x})=\frac{L(\theta;\,\boldsymbol{x})\pi(\boldsymbol{\theta})}{\int_{\Theta}L(\vartheta;\,\boldsymbol{x})\pi(\vartheta)\mathrm{d}\vartheta}.$

This justifies the choice of $L$ appearing in Definition 2.3. ∎

###  4.4 Spatial binning encodes less statistical information

The previous section provides a formal derivation of the point process
likelihood (8) by considering increasingly fine-scale partitions of the
spatial domain. In this subsection, we motivate this further by showing that
the resulting likelihood encodes higher statistical information than any
spatial binning of the process, stated as the following proposition.

######  Proposition 4.3.

Let $u(x;\theta)$ be the intensity of a Poisson point process on a domain
$x\in\Omega$ with underlying parameter $\theta$. Let $\\{A_{j}\\}_{j=1}^{j=K}$
be any non-overlapping partition of the domain and denote observations of the
binned process

$N_{j}=N(A_{j})\sim\mathrm{Poisson}\left(\int_{A_{j}}u(x)\mathrm{d}x\right),\qquad
j=1,\ldots,K$ (38)

with corresponding log-likelihood
$\ell_{\mathrm{bin}}(\theta;\,\boldsymbol{x})$. Let
$I_{\mathrm{bin}}(\theta)=-\mathbb{E}_{\theta}[\partial_{\theta\theta}\ell_{\mathrm{bin}}(\theta;\,\boldsymbol{x})]$
be the information for the binned process and
$I_{\mathrm{unbin}}(\theta)=-\mathbb{E}_{\theta}[\partial_{\theta\theta}\log
L(\theta;\,\boldsymbol{x})]$ be the point process information from likelihood
(8). Then,

$I_{\mathrm{bin}}(\theta)\leq I_{\mathrm{unbin}}(\theta).$

###### Proof.

The likelihood for the binned process is this process is (35), and the log-
likelihood is

$\displaystyle\ell_{\mathrm{bin}}(\theta;\,\boldsymbol{x})$
$\displaystyle\coloneqq\log L_{\mathrm{bin}}(\theta;\,\boldsymbol{x})$ (39)
$\displaystyle=-\int_{\Omega}u(x;\,\theta)\mathrm{d}x+\sum_{j}n_{j}(\boldsymbol{x})\log\int_{A_{j}}u(x;\,\theta)\mathrm{d}x+\sum_{j}\log
n_{j}(\boldsymbol{x})!.$ (40)

For a parameter $\theta$,

$\displaystyle I_{\mathrm{bin}}(\theta)$
$\displaystyle=-\mathbb{E}_{\theta}[\partial_{\theta\theta}\ell_{\mathrm{bin}}(\theta;\,\boldsymbol{x})]$
(41)
$\displaystyle=\int_{\Omega}u_{\theta\theta}(x;\,\theta)\mathrm{d}x-\sum_{j}\mathbb{E}[n_{j}]\partial_{\theta\theta}\log\left(\int_{A_{j}}u(x;\,\theta)\mathrm{d}x\right)$
$\displaystyle=\int_{\Omega}u_{\theta\theta}(x;\,\theta)\mathrm{d}x-\sum_{j}\left(\int_{A_{j}}u(x;\,\theta)\mathrm{d}x\right)\partial_{\theta\theta}\log\left(\int_{A_{j}}u(x;\,\theta)\mathrm{d}x\right),$
$\displaystyle=\sum_{j}\frac{\left(\int_{A_{j}}u_{\theta}(x;\,\theta)\mathrm{d}x\right)^{2}}{\int_{A_{j}}u(x;\,\theta)\mathrm{d}x}.$

where we have used $\partial_{\theta\theta}n_{j}=0$ and
$\mathbb{E}_{\theta}[n_{j}]=\int_{A_{j}}u(x;\,\theta)\mathrm{d}x$. We can
compare this directly to the unbinned information obtained via Campbell’s
formula (16), restated as

$I_{\mathrm{unbin}}(\theta)=\int_{\Omega}u_{\theta\theta}(x;\,\theta)\mathrm{d}x-\int_{\Omega}u\partial_{\theta\theta}(x;\,\theta)\log
u(x;\,\theta)\mathrm{d}x=\int_{\Omega}\frac{u_{\theta}(x;\,\theta)^{2}}{u(x;\,\theta)}\mathrm{d}x.$
(42)

To prove $I_{\mathrm{bin}}\leq I_{\mathrm{unbin}}$, it must be shown that

$\sum_{j}\frac{\left(\int_{A_{j}}u_{\theta}(x;\,\theta)\mathrm{d}x\right)^{2}}{\int_{A_{j}}u(x;\,\theta)\mathrm{d}x}\leq\int_{\Omega}\frac{u_{\theta}(x;\,\theta)^{2}}{u(x;\,\theta)}\mathrm{d}x=\sum_{j}\int_{A_{j}}\frac{u_{\theta}(x;\,\theta)^{2}}{u(x;\,\theta)}\mathrm{d}x,$
(43)

and noting this must hold for any $A_{j}$, this further reduces to showing

$\frac{\left(\int_{A_{j}}u_{\theta}(x;\,\theta)\mathrm{d}x\right)^{2}}{\int_{A_{j}}u(x;\,\theta)\mathrm{d}x}\leq\int_{A_{j}}\frac{u_{\theta}(x;\,\theta)^{2}}{u(x;\,\theta)}\mathrm{d}x.$
(44)

Denote the inner product $\langle f,g\rangle_{u}\coloneqq\int
f(x)g(x)u(x)\mathrm{d}x$. Then the assertion (44) can be written

$\frac{\langle\partial_{\theta}\log u,1\rangle_{u}^{2}}{\langle
1,1\rangle_{u}}\leq\langle\partial_{\theta}\log u,\partial_{\theta}\log
u\rangle_{u}.$ (45)

In this form (44) is seen to be an expression of the Cauchy-Schwarz
inequality, and the claim $I_{\mathrm{bin}}\leq I_{\mathrm{unbin}}$ is proved.
∎

Cauchy-Schwarz provides an interpretation of this result. Equality is achieved
when $\partial_{\theta}\log u\propto\mathrm{constant}$, or when
$\partial_{\theta}u\propto u$. Intuitively, this says that scale parameters
(such as $\lambda$ in the main text) are insensitive to binning, but
parameters that modulate the spatial variation of the intensity magnify the
information loss from binning.

##  5 Discussion

In summary, we have established new understanding in the theory and practice
of inferring dynamics of spatially stochastic birth-death-diffusion processes
in heterogeneous domains from particle snapshot data alone. This pursuit was
motivated by technical challenges arising in the study of spatial
transcriptomics imaging data, where observable data often consists of
positions of a stochastic population of RNA molecules, at a single time (due
to cell fixation), in a cell-specific geometry. Using an argument based on the
occupation measure of particles, the principal result is a rigorous
statistical connection between particle paths and their observation as a
Poisson point process with intensity governed by a birth-death-diffusion PDE
with appropriate boundaries. This explicit statistical description enables
inference that faithfully incorporates often-ignored cell-to-cell
heterogeneities. Surprisingly, we find that variation in source location and
domain size can increase the accuracy of inference on the underlying dynamics
in the regime where particles are long-lived relative to their typical
diffusive exit time from the domain. Our study has broader context and future
directions in the areas of inverse problems and mathematical modeling of gene
expression data.

As an inverse-problem, the most important aspect of the work is the formal
derivation and justification for the point process likelihood. This non-
trivial likelihood must describe a stochastic number of particles and their
position, in contrast to previous inference studies with single particle paths
[BF16, FC20] or population positions with fixed number [Che+21]. Our work
clarifies how point processes are the natural resolution of this technical
hurdle. Moreover, we show that the point process description encodes more
statistical information than spatial binning used in other inference of
stochastic reaction-diffusion systems [Dew+10]. We are far from the first to
consider mechanistic models in partnership with spatial point processes in
biological systems, e.g., [PL15, Ova+14]. The distinction of these studies in
the ecological literature is subtle but important: in our work, the point
process is an emergent statistical description from a birth-death-diffusion
process rather than an assumed statistical model for observations. Although
the PDE-constraint on the intensity bares resemblance to classical PDE inverse
problems [ABT18, Fle+20], we emphasize the key difference in the observation
noise. Rather than observations of the PDE solution (often boundary data) plus
additive noise, our inverse problem corresponds to the observations of
stochastic particle positions that serve as a proxy for the PDE solution. We
provide a preliminary investigation of the statistical properties of this
inverse problem, but there are several avenues of promising investigation in
discerning similarities and differences in this PDE inverse problem against
their classical counterpart.

On the computational side, this connection to other inverse problems suggests
there are are likely improvements toward more efficient inference. Such
improvements will likely be necessary for large-scale inference on now
available tissue-scale data [Eng+19]. The primary bottleneck of our approach
is each evaluation of the likelihood requires solving the PDE model over all
observed domains. The simplicity of the 1D model allowed for an analytical
solution to the PDE and straightforward MCMC sampling. However, these luxuries
break down quickly for more realistic investigations. Therefore, one may be
able to construct improved samplers seen in other gene model
inference[Kil+23], or leverage evaluations of the likelihood in approximate
techniques such as variational Bayesian methods [BKM17, Ace18]. However, these
approaches still suffer the bottleneck of each likelihood evaluation requiring
the solution to the model over all observed domains. An alternative avenue to
circumvent this may be the use of neural network or related approaches. Such
techniques have served useful in inverse problems on stochastic chemical
reactions [Jia+21], classical PDE inverse problems [RPK19], and particle
positions [Che+21], including with domain generalization [Cao+23]. Therefore,
this seems like a promising avenue of study for this problem at the
intersection of these topics. Lastly, a variety of investigations remain in
adapting the zoo of techniques in the toolbelt for inverse problems to this
setup, ranging from model selection [Nar+21] to optimal experimental design
[Rya+16].

The steady-state birth-death-diffusion model was chosen at a deliberate level
of complexity to distill the salient features of the data: snapshots of
fluctuating populations in a domain. However, there is considerable work that
must be done in incorporating known biological complexities before our
framework can be used on experimental smFISH data. We acknowledge the context
of a vast literature on forward modeling of various non-spatial aspects of
stochastic gene expression e.g., [Pau05, SS13, CG20, GP20, KM23]. The use of
these stochastic models in inference for transcriptomics data is also
increasingly prevalent and sophisticated [Her+17, Góm+17, Luo+22]. We hope
that our investigate serves as a segue for the field to transition into
spatially resolved models that can be linked to this already available data. A
non-exhaustive list of interesting complexities that may be included in our
model are as follows. The production of RNA for many genes is known to be
bursty [JE18, Don+19], with short stochastic intermittent windows of
production. A spatial bursting model may be tractable with generating function
techniques [CST12] that serve fruitful for analytical understanding of non-
spatial models [GP20, GP22]. On the geometric side, the model must be extended
to biological datasets of 2D and 3D FISH imaging [Ome+19]. The models can also
be extended to multiple compartments (the nucleus, cytoplasm) and serve as a
platform for investigating the nature of RNA export through different possible
boundary conditions, e.g., Robin boundary conditions arising through
stochastic gating of particles [LK15] that may approximate the complex
dynamics of RNA export [KH07, Azi+14]. Lastly, the movement model could be
extended to include subdiffusive motion seen [Cab+06], known to profoundly
dictate spatial patterns, especially near boundaries [Hol19].

In summary, our work serves as a basis for more detailed investigations of
inferring dynamics from spatial dynamics from snapshots of stochastic particle
populations, especially those arising from spatial transcriptomics data. The
pursuit of understanding this biological data could only be accomplished by
distilling it to its essential components and constructing that reflects these
complexities. The analytical tractability of the model allowed for us to
carefully examine the role of cell-to-cell heterogeneities which are surely
present in a wide array of biological datasets. This surprising insight of
robustness in inference arising from noise contributes a new dimension to the
broader theme in mathematical biology of how noise can sharpen signals
[Hän02]. In an era where phenomenological statistical analysis of biological
data is commonplace, we hope that our study will inspire others to consider
rigorous modeling of biological data to discover other surprises and insights.

##  6 Acknowledgements

FD acknowledges NIH support 1DP2GM149554. The work of RBL and SAM was
supported by the Laboratory Directed Research and Development program at
Sandia National Laboratories, a multimission laboratory managed and operated
by National Technology and Engineering Solutions of Sandia LLC, a wholly owned
subsidiary of Honeywell International Inc. for the U.S. Department of Energy’s
National Nuclear Security Administration under contract DE-NA0003525.

##  Appendix A Numerical implementation details

###  A.1 Stochastic simulation algorithm for synthetic data

Algorithm 1 describes the generation procedure for synthetic data used in the
inference setups of Figs. 5 and 6. Importantly, we simulate trajectories of
the particle process, and do inference using the point process theory from
Theorem 2.1. The simulation procedure is exact except for a small
$\Delta_{t}$-controlled approximation error in computing whether a Brownian
path has exited the domain [SG19]. In all simulations, we take
$\Delta_{t}=10^{-5}$, several orders of magnitude smaller than any other
timescale.

Algorithm 1  Stochastic simulation to generate particle positions for
synthetic data.

SimulateParticles$k_{\text{birth}},k_{\text{death}},D,z,L,T,\Delta_{t}$

1:$N\text{births}\sim\operatorname{Poisson}(\kappa_{\mathrm{+}}T)$
$n=1,\ldots,N{\text{births}}$

2:$T\text{birth}_{n}\sim\operatorname{unif}(0,T)$

3:$\tau\text{decay}_{n}\sim\operatorname{exp}(1/\kappa_{\mathrm{-}})$

4:$T\text{decay}_{n}\leftarrow T\text{birth}_{n}+\tau\text{decay}_{n}$
$T\text{decay}_{n}>T$$\triangleright$ only simulate if particle is possibly
alive $T$

5:$x\leftarrow z$, $t_{n}\leftarrow T\text{birth}_{n}$ $\triangleright$ birth
particle at source location

6:$\text{survived}_{n}\leftarrow 1$ $t_{n}\leq T$

7:$x\leftarrow x+\sqrt{2D\Delta_{t}}\cdot\texttt{randn}$ $\triangleright$
Brownian step

8:$t_{n}\leftarrow t_{n}+\Delta_{t}$ $x<0$ or $x>L$ $\triangleright$ if exits
domain

9:$\text{survived}_{n}\leftarrow 0$ break $\triangleright$ instantly absorbed,
stop simulating $\text{survived}_{n}$ $\triangleright$ if alive at observation
$T$

10:append $x$ $\triangleright$ store final position

11:return $x_{1},\ldots,x_{N}$

###  A.2 Metropolis-in-Gibbs sampling for Bayesian inference

See Section 3.2 for an explanation of prior choices for $\lambda,\mu$. The
choice of prior for $\lambda$ provides a (marginal) conjugate posterior that
can be computed analytically, so a Gibbs-in-Metropolis MCMC approach is used
and described in Algorithm 2.

Algorithm 2  MCMC Gibbs-in-Metropolis sampling of posterior distribution.

PosteriorSampMCMC$X,\sigma,\alpha_{\lambda},\beta_{\lambda},\mu_{\mathrm{max}}$

1:init
$\lambda_{0}\sim\operatorname{Gamma}(\alpha_{\lambda},\beta_{\lambda})$,
$\mu_{0}\sim\operatorname{unif}(0,\mu_{\mathrm{max}})$ $s=1\ldots,S$ MH step
for $\mu$ first

2:$\mu_{*}\leftarrow\mu+\sigma\cdot\texttt{randn}$ $\mu_{*}<0$ or
$\mu_{*}>\mu_{\text{max}}$ $\mu_{*}\leftarrow\mu+\sigma\cdot\texttt{randn}$
$\triangleright$ only accept proposals for $\mu$ in range of prior

3:$\pi_{*}\leftarrow\ell(\lambda_{s-1},\mu_{*},X)$,
$\pi_{s-1}\leftarrow\ell(\lambda_{s-1},\mu_{s-1},X)$ $\triangleright$ proposal
& current log-likelihood

4:$p\leftarrow\min\\{\exp(\pi_{*}-\pi_{s-1}),1\\}$ $\texttt{rand}<p$
$\mu_{s}\leftarrow\mu_{*}$ $\triangleright$ accept proposal
$\mu_{s}\leftarrow\mu_{s-1}$ $\triangleright$ keep current Gibbs sample
$\lambda$ using Gamma conjugacy

5:$\alpha,\beta$ $\leftarrow$ using prior $\alpha_{\lambda},\beta_{\lambda}$,
data $X$, and $\mu_{s}$ in (21) $\triangleright$ Gibbs step uses current value
of $\mu$

6:$\lambda_{s}\sim\operatorname{Gamma}(\alpha,\beta)$

7:return $\lambda_{0},\ldots,\lambda_{S},\quad\mu_{0},\ldots,\mu_{S}$

Figure 7: Rubin-Gelman $\hat{R}$-1 demonstrating convergence of MCMC (on the
setup in Fig 5 with $M=10$) for two chains. Convergence is considered good
when $\hat{R}-1<10^{-1}$, achieved in approximately $10^{3}$ MCMC steps.

As is standard practice [Gel+95], we do sequence of warm-up phases to target
the range of acceptance probabilities between $0.3$ and $0.5$, starting with
$\sigma=1$ and rounds of 500 steps. After each warm-up round, if the
acceptance frequency is below the target range, $\sigma\to\sigma/2$, and if
above $\sigma\to 2\sigma$. Once in the range, the ending values of
$\lambda,\mu$ are used as initial values and the full algorithm is run with
this step size $\sigma$ for much greater number of steps. The values from the
chains for this burn-in phase are discarded.

To choose the number of MCMC steps, we investigate the Gelman-Rubin diagnostic
$\hat{R}$ [Gel+95] for the presumed worst-case scenario of $M=10$ snapshots in
the setup of Fig. 5. For $M=10$, we initialize multiple chains and compute
$\hat{R}$ as a function of the number of MCMC steps, shown in Fig. 7.
Convergence is considered satisfactory for $\hat{R}<1.1$ [Gel+95], which is
achieved around $10^{3}$ MCMC steps. We take $10^{5}$ MCMC steps for all
setups in the main text and do not monitor convergence further.

## References

  * [Ace18] Luigi Acerbi “Variational Bayesian Monte Carlo” In _Proceedings of the 32nd International Conference on Neural Information Processing Systems_ , NeurIPS’18 Montréal, Canada: Curran Associates Inc., 2018, pp. 8223–8233
  * [ABT18] Richard C Aster, Brian Borchers and Clifford H Thurber “Parameter estimation and inverse problems” Elsevier, 2018 DOI: 10.1016/c2009-0-61134-x
  * [ADT07] Anthony Atkinson, Alexander Donev and Randall Tobias “Optimum experimental designs, with SAS” OUP Oxford, 2007 DOI: 10.1093/oso/9780199296590.001.0001
  * [Azi+14] Mohammad Azimi, Evgeny Bulat, Karsten Weis and Mohammad R.. Mofrad “An agent-based model for mRNA export through the nuclear pore complex” In _Molecular Biology of the Cell_ 25.22, 2014, pp. 3643–3653 DOI: 10.1091/mbc.e14-06-1065
  * [Bad01] A.. Baddeley “Likelihoods and pseudolikelihoods for Markov spatial processes” In _Institute of Mathematical Statistics Lecture Notes - Monograph Series_ Beachwood, OH: Institute of Mathematical Statistics, 2001, pp. 21–49 DOI: 10.1214/lnms/1215090061
  * [BF16] Jason Bernstein and John Fricks “Analysis of single particle diffusion with transient binding using particle filtering” In _Journal of theoretical biology_ 401 Elsevier, 2016, pp. 109–121 DOI: 10.1016/j.jtbi.2016.04.013
  * [Bis+21] Jeetayu Biswas, Weihan Li, Robert H Singer and Robert A Coleman “Imaging organization of RNA processing within the nucleus” In _Cold Spring Harbor Perspectives in Biology_ 13.12 Cold Spring Harbor Lab, 2021, pp. a039453 DOI: 10.1101/cshperspect.a039453.
  * [BKM17] David M. Blei, Alp Kucukelbir and Jon D. McAuliffe “Variational Inference: A Review for Statisticians” In _Journal of the American Statistical Association_ 112.518, 2017, pp. 859–877 DOI: 10.1080/01621459.2017.1285773
  * [Bow+22] Jonathan R Bowles, Caroline Hoppe, Hilary L Ashe and Magnus Rattray “Scalable Inference of Transcriptional Kinetic Parameters from MS2 Time Series Data” In _Bioinformatics_ 38.4, 2022, pp. 1030–1036 DOI: 10.1093/bioinformatics/btab765
  * [Bro+20] Alexander P Browning et al. “Identifiability analysis for stochastic differential equation models in systems biology” In _Journal of the Royal Society Interface_ 17.173 The Royal Society, 2020, pp. 20200652 DOI: doi.org/10.1098/rsif.2020.0652
  * [Cab+06] Ghislain G Cabal et al. “SAGA interacting factors confine sub-diffusion of transcribed genes to the nuclear envelope” In _Nature_ 441.7094 Nature Publishing Group UK London, 2006, pp. 770–773 DOI: 10.1038/nature04752
  * [CFX06] Long Cai, Nir Friedman and X. Xie “Stochastic protein expression in individual cells at the single molecule level” In _Nature_ 440.7082, 2006, pp. 358–362 DOI: 10.1038/nature04599
  * [Cao+23] Zhixing Cao et al. “Efficient and Scalable Prediction of Spatio-Temporal Stochastic Gene Expression in Cells and Tissues Using Graph Neural Networks”, 2023 DOI: 10.1101/2023.02.28.530379
  * [CG20] Zhixing Cao and Ramon Grima “Analytical distributions for detailed models of stochastic gene expression in eukaryotic cells” In _Proceedings of the National Academy of Sciences_ 117.9 National Acad Sciences, 2020, pp. 4682–4692 DOI: 10.1073/pnas.1910888117
  * [CE22] Lucia Cassella and Anne Ephrussi “Subcellular spatial transcriptomics identifies three mechanistically different classes of localizing RNAs” In _Nature Communications_ 13.1, 2022, pp. 6355 DOI: 10.1038/s41467-022-34004-2
  * [CP23] Tara Chari and Lior Pachter “The specious art of single-cell genomics” In _PLoS Computational Biology_ 19.8 Public Library of Science San Francisco, CA USA, 2023, pp. e1011288 DOI: 10.1371/journal.pcbi.1011288
  * [CG78] S Chaturvedi and CW Gardiner “The Poisson Representation. II Two-time correlation functions” In _Journal of Statistical Physics_ 18.5 Springer, 1978, pp. 501–522 DOI: 10.1007/bf01014520
  * [Cha+77] S Chaturvedi, CW Gardiner, IS Matheson and DF Walls “Stochastic analysis of a chemical reaction with spatial and temporal structures” In _Journal of Statistical Physics_ 17.6 Springer, 1977, pp. 469–489 DOI: 10.1007/bf01014350
  * [Che+14] Alex Chen et al. “Transient antibody-mucin interactions produce a dynamic molecular shield against viral invasion” In _Biophysical journal_ 106.9 Elsevier, 2014, pp. 2028–2036 DOI: 10.1016/j.bpj.2014.02.038
  * [Che+21] Xiaoli Chen, Liu Yang, Jinqiao Duan and George Em Karniadakis “Solving Inverse Stochastic Problems from Discrete Particle Observations Using the Fokker–Planck Equation and Physics-Informed Neural Networks” In _SIAM Journal on Scientific Computing_ 43.3, 2021, pp. B811–B830 DOI: 10.1137/20M1360153
  * [Cla22] Daniel E Clark “A Cramér Rao bound for point processes” In _IEEE Transactions on Information Theory_ 68.4 IEEE, 2022, pp. 2147–2155 DOI: 10.1109/tit.2022.3140374
  * [CLL14] Laurent Cognet, Cécile Leduc and Brahim Lounis “Advances in live-cell single-particle tracking and dynamic super-resolution imaging” In _Current opinion in chemical biology_ 20 Elsevier, 2014, pp. 78–85 DOI: 10.1016/j.cbpa.2014.04.015
  * [Cot+23] Allison Coté et al. “Post-transcriptional splicing can occur in a slow-moving zone around the gene” In _bioRxiv_ Cold Spring Harbor Laboratory, 2023 DOI: 10.1101/2020.04.06.028092
  * [CST12] David Cottrell, Peter S. Swain and Paul F. Tupper “Stochastic branching-diffusion models for gene expression” In _Proceedings of the National Academy of Sciences_ 109.25, 2012, pp. 9699–9704 DOI: 10.1073/pnas.1201103109
  * [CM77] David Roxbee Cox and Hilton David Miller “The theory of stochastic processes” CRC press, 1977 DOI: 10.1201/9780203719152
  * [DV+03] Daryl J Daley and David Vere-Jones “An introduction to the theory of point processes: volume I: elementary theory and methods” Springer, 2003
  * [Dew+10] Michael A Dewar, Visakan Kadirkamanathan, Manfred Opper and Guido Sanguinetti “Parameter estimation and inference for stochastic reaction-diffusion systems: application to morphogenesis in D. melanogaster” In _BMC Systems Biology_ 4.1 BioMed Central, 2010, pp. 1–9 DOI: 10.1186/1752-0509-4-21
  * [DE19] Fangyuan Ding and Michael B. Elowitz “Constitutive splicing and economies of scale in gene expression” In _Nature Structural and Molecular Biology_ 26.6 Springer US, 2019, pp. 424–432 DOI: 10.1038/s41594-019-0226-x
  * [Din+22] Fangyuan Ding et al. “Dynamics and functional roles of splicing factor autoregulation” In _Cell reports_ 39.12 Elsevier, 2022 DOI: 10.1016/j.celrep.2022.110985
  * [Don+19] Benjamin T Donovan et al. “Live-cell imaging reveals the interplay between transcription factors, nucleosomes, and bursting” In _The EMBO journal_ 38.12, 2019, pp. e100809 DOI: 10.15252/embj.2018100809
  * [Dur96] Richard Durrett “Stochastic calculus: a practical introduction” CRC press, 1996
  * [Eng+19] Chee-Huat Linus Eng et al. “Transcriptome-scale super-resolved imaging in tissues by RNA seqFISH+” In _Nature_ 568.7751, 2019, pp. 235–239 DOI: 10.1038/s41586-019-1049-y
  * [FC20] Rebeca Cardim Falcao and Daniel Coombs “Diffusion analysis of single particle trajectories in a Bayesian nonparametrics framework” In _Physical Biology_ 17.2 IOP Publishing, 2020, pp. 025001 DOI: 10.1088/1478-3975/ab64b3
  * [Fle+20] Mark B Flegg et al. “Parameter estimation for a point-source diffusion-decay morphogen model” In _Journal of Mathematical Biology_ 80.7 Springer, 2020, pp. 2227–2255 DOI: 10.1007/s00285-020-01494-x
  * [FM19] Zachary R Fox and Brian Munsky “The finite state projection based Fisher information matrix approach to estimate information and optimize single-cell experiments” In _PLoS computational biology_ 15.1 Public Library of Science San Francisco, CA USA, 2019, pp. e1006365 DOI: 10.1371/journal.pcbi.1006365
  * [Fu+22] Xiaoming Fu et al. “Quantifying how post-transcriptional noise and gene copy number variation bias transcriptional parameter inference from mRNA distributions” In _Elife_ 11 eLife Sciences Publications Limited, 2022, pp. e82493 DOI: 10.1101/2021.11.09.467882
  * [GB15] Attila Gábor and Julio R Banga “Robust and efficient parameter estimation in dynamic models of biological systems” In _BMC systems biology_ 9.1 BioMed Central, 2015, pp. 1–25 DOI: 10.1186/s12918-015-0219-2
  * [Gar04] Crispin W Gardiner “Handbook of Stochastic Methods” Springer Berlin, 2004
  * [GC77] Crispin W Gardiner and S Chaturvedi “The Poisson representation. I. A new technique for chemical master equations” In _Journal of Statistical Physics_ 17.6 Springer, 1977, pp. 429–468 DOI: 10.1007/bf01014349
  * [Gel+95] Andrew Gelman, John B Carlin, Hal S Stern and Donald B Rubin “Bayesian Data Analysis” ChapmanHall/CRC, 1995 DOI: 10.1201/9780429258480
  * [Góm+17] Mariana Gómez-Schiavon, Liang-Fu Chen, Anne E West and Nicolas E Buchler “BayFish: Bayesian inference of transcription dynamics from population snapshots of single-molecule RNA FISH in single cells” In _Genome Biology_ 18 Springer, 2017, pp. 1–12 DOI: 10.1186/s13059-017-1297-9
  * [Gor+22] Gennady Gorin, Meichen Fang, Tara Chari and Lior Pachter “RNA Velocity Unraveled” In _PLoS Computational Biology_ 18.9 Public Library of Science San Francisco, CA USA, 2022, pp. e1010492 DOI: 10.1101/2022.02.12.480214
  * [GP20] Gennady Gorin and Lior Pachter “Special function methods for bursty models of transcription” In _Phys. Rev. E_ 102 American Physical Society, 2020, pp. 022409 DOI: 10.1103/PhysRevE.102.022409
  * [GP22] Gennady Gorin and Lior Pachter “Modeling bursty transcription and splicing with the chemical master equation” In _Biophysical Journal_ 121.6 Elsevier, 2022, pp. 1056–1069 DOI: 10.1016/j.bpj.2022.02.004
  * [Gor+22a] Gennady Gorin, John J Vastola, Meichen Fang and Lior Pachter “Interpretable and tractable models of transcriptional noise for the rational design of single-molecule quantification experiments” In _Nature Communications_ 13.1 Nature Publishing Group UK London, 2022, pp. 7620 DOI: 10.1038/s41467-022-34857-7
  * [Gup+22] Anika Gupta et al. “Inferring Gene Regulation from Stochastic Transcriptional Variation across Single Cells at Steady State” In _Proceedings of the National Academy of Sciences_ 119.34, 2022, pp. e2207392119 DOI: 10.1073/pnas.2207392119
  * [Hän02] Peter Hänggi “Stochastic resonance in biology how noise can enhance detection of weak signals and help improve biological information processing” In _ChemPhysChem_ 3.3 Wiley Online Library, 2002, pp. 285–290 DOI: doi.org/10.1002/1439-7641(20020315)3:3<285::AID-CPHC285>3.0.CO;2-A
  * [Hei+17] Stephanie Heinrich et al. “Temporal and spatial regulation of mRNA export: Single particle RNA-imaging provides new tools and insights” In _BioEssays_ 39.2, 2017, pp. 1600124 DOI: 10.1002/bies.201600124
  * [Her+17] Ulysse Herbach, Arnaud Bonnaffoux, Thibault Espinasse and Olivier Gandrillon “Inferring gene regulatory networks from single-cell data: a mechanistic approach” In _BMC Systems Biology_ 11.1 BioMed Central, 2017, pp. 1–15 DOI: 10.1186/s12918-017-0487-0
  * [Hol19] William R Holmes “Subdiffusive dynamics lead to depleted particle densities near cellular borders” In _Biophysical Journal_ 116.8 Elsevier, 2019, pp. 1538–1546 DOI: 10.1016/j.bpj.2019.02.021
  * [Imb+22] Arthur Imbert et al. “FISH-quant v2: A Scalable and Modular Tool for smFISH Image Analysis” In _RNA_ 28.6, 2022, pp. 786–795 DOI: 10.1261/rna.079073.121
  * [Jia+21] Qingchao Jiang et al. “Neural network aided approximation and parameter inference of non-Markovian models of gene expression” In _Nature Communications_ 12.1, 2021, pp. 2618 DOI: 10.1038/s41467-021-22919-1
  * [JE18] Daniel Jones and Johan Elf “Bursting onto the scene? Exploring stochastic mRNA production in bacteria” In _Current Opinion in Microbiology_ 45 Elsevier, 2018, pp. 124–130 DOI: 10.1016/j.mib.2018.04.001
  * [KM23] Bhargav R. Karamched and Christopher E. Miles “Stochastic switching of delayed feedback suppresses oscillations in genetic regulatory systems” In _Journal of The Royal Society Interface_ 20.203, 2023, pp. 20230059 DOI: 10.1098/rsif.2023.0059
  * [Kil+23] Zeliha Kilic, Max Schweiger, Camille Moyer and Steve Pressé “Monte Carlo samplers for efficient network inference” In _PLOS Computational Biology_ 19.7 Public Library of Science San Francisco, CA USA, 2023, pp. e1011256 DOI: 10.1371/journal.pcbi.1011256
  * [Kil+23a] Zeliha Kilic et al. “Gene expression model inference from snapshot RNA data using Bayesian non-parametrics” In _Nature Computational Science_ 3.2 Nature Publishing Group US New York, 2023, pp. 174–183 DOI: 10.1038/s43588-022-00392-0
  * [KH07] Alwin Köhler and Ed Hurt “Exporting RNA from the nucleus to the cytoplasm” In _Nature Reviews Molecular Cell Biology_ 8.10 Nature Publishing Group, 2007, pp. 761–773 DOI: 10.1038/nrm2255
  * [La ̵+18] Gioele La Manno et al. “RNA velocity of single cells” In _Nature_ 560.7719 Nature Publishing Group, 2018, pp. 494–498 DOI: 10.1038/s41586-018-0414-6
  * [Lam+20] Nicholas C Lammers, Yang Joon Kim, Jiaxi Zhao and Hernan G Garcia “A matter of time: Using dynamics and theory to uncover mechanisms of transcriptional bursting” In _Current opinion in cell biology_ 67 Elsevier, 2020, pp. 147–157 DOI: 10.1016/j.ceb.2020.08.001
  * [LK15] Sean D Lawley and James P Keener “A new derivation of Robin boundary conditions through homogenization of a stochastically switching boundary” In _SIAM Journal on Applied Dynamical Systems_ 14.4 SIAM, 2015, pp. 1845–1867 DOI: 10.1137/15M1015182
  * [LDS05] Caroline Lemerle, Barbara Di Ventura and Luis Serrano “Space as the final frontier in stochastic simulations of biological systems” In _FEBS letters_ 579.8 Elsevier, 2005, pp. 1789–1794 DOI: 10.1016/j.febslet.2005.02.009
  * [Luo+22] Xizhi Luo, Fei Qin, Feifei Xiao and Guoshuai Cai “BISC: accurate inference of transcriptional bursting kinetics from single-cell transcriptomic data” In _Briefings in Bioinformatics_ 23.6 Oxford Academic, 2022 DOI: 10.1093/bib/bbac464.
  * [Mah+22] Clarence K Mah et al. “Bento: A toolkit for subcellular analysis of spatial transcriptomics data” In _BioRxiv_ Cold Spring Harbor Laboratory, 2022, pp. 2022–06 DOI: 10.1101/2022.06.10.495510
  * [MG15] Carlo Manzo and Maria F Garcia-Parajo “A review of progress in single particle tracking: from methods to biophysical insights” In _Reports on progress in physics_ 78.12 IOP Publishing, 2015, pp. 124601
  * [Mar21] Vivien Marx “Method of the Year: spatially resolved transcriptomics” In _Nature methods_ 18.1 Nature Publishing Group US New York, 2021, pp. 9–14 DOI: 10.1038/s41592-020-01033-y
  * [May+20] Kristen R Maynard et al. “Dotdotdot: An Automated Approach to Quantify Multiplex Single Molecule Fluorescent in Situ Hybridization (smFISH) Images in Complex Tissues” In _Nucleic Acids Research_ 48.11, 2020, pp. e66 DOI: 10.1093/nar/gkaa312
  * [McL+20] Grace A McLaughlin et al. “Spatial heterogeneity of the cytosol revealed by machine learning-based 3D particle tracking” In _Molecular Biology of the Cell_ 31.14 Am Soc Cell Biol, 2020, pp. 1498–1511 DOI: 10.1091/mbc.E20-03-0210
  * [Mel+16] John WR Mellnik et al. “Maximum likelihood estimation for single particle, passive microrheology data with drift” In _Journal of Rheology_ 60.3 The Society of Rheology, 2016, pp. 379–392 DOI: 10.1122/1.4943988
  * [Moe07] William E Moerner “New directions in single-molecule imaging and analysis” In _Proceedings of the National Academy of Sciences_ 104.31 National Acad Sciences, 2007, pp. 12596–12602 DOI: 10.1073/pnas.0610081104
  * [MP22] Lambda Moses and Lior Pachter “Museum of spatial transcriptomics” In _Nature Methods_ 19.5, 2022, pp. 534–546 DOI: 10.1038/s41592-022-01409-2
  * [MFN15] Brian Munsky, Zachary Fox and Gregor Neuert “Integrating single-molecule experiments and discrete stochastic models to understand heterogeneous gene transcription dynamics” In _Methods_ 85 Elsevier, 2015, pp. 12–21 DOI: 10.1016/j.ymeth.2015.06.009
  * [Nar+21] John T Nardini, Ruth E Baker, Matthew J Simpson and Kevin B Flores “Learning differential equation models from stochastic agent-based model simulations” In _Journal of the Royal Society Interface_ 18.176 The Royal Society, 2021, pp. 20200987 DOI: 10.1098/rsif.2020.0987
  * [Ome+19] Manja Omerzu et al. “Three-Dimensional Analysis of Single Molecule FISH in Human Colon Organoids” In _Biology Open_ 8.8, 2019, pp. bio042812 DOI: 10.1242/bio.042812
  * [Osu+19] Olaolu Osunbayo et al. “Complex nearly immotile behaviour of enzymatically driven cargos” In _Soft Matter_ 15.8, 2019, pp. 1847–1852 DOI: 10.1039/c8sm01893f
  * [Ova+14] Otso Ovaskainen et al. “A general mathematical framework for the analysis of spatiotemporal point processes” In _Theoretical Ecology_ 7 Springer, 2014, pp. 101–113 DOI: 10.1007/s12080-013-0202-8
  * [Pau05] Johan Paulsson “Models of stochastic gene expression” In _Physics of life reviews_ 2.2 Elsevier, 2005, pp. 157–175 DOI: 10.1016/j.plrev.2005.03.003
  * [PL15] Michael J Plank and Richard Law “Spatial point processes and moment dynamics in the life sciences: a parsimonious derivation and some extensions” In _Bulletin of mathematical biology_ 77 Springer, 2015, pp. 586–613 DOI: 10.1007/s11538-014-0018-8
  * [QSE91] Hong Qian, Michael P Sheetz and Elliot L Elson “Single particle tracking. Analysis of diffusion and flow in two-dimensional systems” In _Biophysical Journal_ 60.4 Elsevier, 1991, pp. 910–921 DOI: 10.1016/s0006-3495(91)82125-7
  * [RPK19] Maziar Raissi, Paris Perdikaris and George E Karniadakis “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations” In _Journal of Computational physics_ 378 Elsevier, 2019, pp. 686–707 DOI: 10.1016/j.jcp.2018.10.045
  * [Raj+06] Arjun Raj et al. “Stochastic mRNA synthesis in mammalian cells” In _PLoS biology_ 4.10 Public Library of Science San Francisco, USA, 2006, pp. e309 DOI: 10.1371/journal.pbio.0040309
  * [RV08] Arjun Raj and Alexander Van Oudenaarden “Nature, nurture, or chance: stochastic gene expression and its consequences” In _Cell_ 135.2 Elsevier, 2008, pp. 216–226 DOI: 10.1016/j.cell.2008.09.050
  * [Ray+23] Nathan T Rayens, Keisha J Cook, Scott A McKinley and Christine K Payne “Palmitate-mediated disruption of the endoplasmic reticulum decreases intracellular vesicle motility” In _Biophysical Journal_ 122.7 Elsevier, 2023, pp. 1355–1363 DOI: 10.1016/j.bpj.2023.03.001
  * [RW23] Florian Rehfeldt and Matthias Weiss “The random walker’s toolbox for analyzing single-particle tracking data” In _Soft Matter_ Royal Society of Chemistry, 2023 DOI: 10.1039/D3SM00557G
  * [RH11] Susana Rodríguez-Navarro and Ed Hurt “Linking Gene Regulation to mRNA Production and Export” In _Current opinion in Cell Biology_ 23.3 Elsevier, 2011, pp. 302–309 DOI: 10.1016/j.ceb.2010.12.002
  * [Rya+16] Elizabeth G Ryan, Christopher C Drovandi, James M McGree and Anthony N Pettitt “A review of modern computational algorithms for Bayesian optimal design” In _International Statistical Review_ 84.1 Wiley Online Library, 2016, pp. 128–154 DOI: 10.1088/0034-4885/78/12/124601
  * [SGS16] David Schnoerr, Ramon Grima and Guido Sanguinetti “Cox process representation and inference for stochastic reaction–diffusion processes” In _Nature Communications_ 7.1 Nature Publishing Group, 2016, pp. 1–11 DOI: 10.1038/ncomms11729
  * [Sha+13] Sydney M Shaffer, Min-Tzu Wu, Marshall J Levesque and Arjun Raj “Turbo FISH: a method for rapid single molecule RNA FISH” In _PloS One_ 8.9 Public Library of Science San Francisco, USA, 2013, pp. e75120 DOI: 10.1371/journal.pone.0075120
  * [Sha+22] Jesse A Sharp, Alexander P Browning, Kevin Burrage and Matthew J Simpson “Parameter estimation and uncertainty quantification using information geometry” In _Journal of the Royal Society Interface_ 19.189 The Royal Society, 2022, pp. 20210940
  * [SS13] Abhyudai Singh and Mohammad Soltani “Quantifying intrinsic and extrinsic variability in stochastic gene expression models” In _Plos one_ 8.12 Public Library of Science San Francisco, USA, 2013, pp. e84301 DOI: 10.1371/journal.pone.0084301
  * [SG19] Stephen Smith and Ramon Grima “Spatial stochastic intracellular kinetics: A review of modelling approaches” In _Bulletin of Mathematical Biology_ 81.8 Springer, 2019, pp. 2960–3009 DOI: 10.1007/s11538-018-0443-1
  * [SP18] Valentine Svensson and Lior Pachter “RNA velocity: molecular kinetics from single-cell RNA-Seq” In _Molecular Cell_ 72.1 Elsevier, 2018, pp. 7–9 DOI: 10.1016/j.molcel.2018.09.026
  * [VHO20] Milad R Vahid, Bernard Hanzon and Raimund J Ober “Fisher information matrix for single molecules with stochastic trajectories” In _SIAM Journal on Imaging Sciences_ 13.1 SIAM, 2020, pp. 234–264 DOI: 10.1137/19M1242562
  * [Vo+23] Huy D. Vo, Linda S. Forero-Quintero, Luis U. Aguilera and Brian Munsky “Analysis and design of single-cell experiments to harvest fluctuation information while rejecting measurement noise” In _Frontiers in Cell and Developmental Biology_ 11, 2023, pp. 1133994 DOI: 10.3389/fcell.2023.1133994
  * [WSV23] Florin C Walter, Oliver Stegle and Britta Velten “FISHFactor: a probabilistic factor model for spatial transcriptomics data with subcellular resolution” In _Bioinformatics_ 39.5 Oxford University Press, 2023, pp. btad183
  * [Wei+18] Caleb Weinreb et al. “Fundamental limits on dynamic inference from single-cell snapshots” In _Proceedings of the National Academy of Sciences_ 115.10 National Acad Sciences, 2018, pp. E2467–E2476
  * [Wie+21] Franz-Georg Wieland et al. “On structural and practical identifiability” In _Current Opinion in Systems Biology_ 25, 2021, pp. 60–69 DOI: doi.org/10.1016/j.coisb.2021.03.005
  * [Wil+22] Cameron G. Williams et al. “An introduction to spatial transcriptomics for biomedical research” In _Genome Medicine_ 14.1, 2022, pp. 68 DOI: 10.1186/s13073-022-01075-1
  * [YC19] Christopher Yau and Kieran Campbell “Bayesian statistical learning for big data biology” In _Biophysical reviews_ 11.1 Springer, 2019, pp. 95–102 DOI: 10.1007/s12551-019-00499-1