of 20
Assessing Markovian and delay models for snRNA-seq
S1 DERIVATIONS AND SIMULATIONS
S1.1 Fully Markovian master equations
Section 2.1 is based on the theoretical foundations derived in full detail in Gorin et al. (2023). For completeness, we summarize
these ideas here for the special case
=
0
, i.e., with no continuous degrees of freedom.
First, we outline the scope of allowed reactions: gene state interconversion, degradation, isomerization, catalysis, and
production. Gene state interconversion includes the following probability flux contributions:
푑푃
(
,
x
,
)
푑푡
=
∑︁
=
1
푖푠
(
)
(
,
x
,
)
,
or
P
(
x
,
)
푑푡
=
(
)
T
P
,
(S1)
where the second line is a vector shorthand for the first, with entries of
P
indexed by
. As usual,
is a stochastic matrix with
푠푠
:
=
Í
푠푖
. By substituting this derivative into Equation 3 and exploiting the linearity of matrix-vector multiplication,
we obtain the corresponding spectral partial differential equation (PDE) term:
G
(
g
,
)
휕푡
=
(
)
T
G
.
(S2)
Degradation of species
X
contributes first-order terms that scale as
:
푑푃
(
,
x
,
)
푑푡
=
0
[
(
+
1
)
(
,
+
1,
)−
(
,
x
,
)
]
,
(S3)
where
with
do not factor into the computation and are elided. Conversion of
X
to
X
contributes similar terms:
푑푃
(
,
x
,
)
푑푡
=
푖푗

(
+
1
)
(
,
+
1,
1,
)−
(
,
x
,
)

.
(S4)
As above,
with
,
are elided. By differentiating
(
,
g
,
)
and
(
,
g
,
)
with respect to
and
and exploiting the
linearity of the master equation, we obtain the corresponding PDE terms:
휕퐺
(
,
g
,
)
휕푡
=
0
(
1
)
휕퐺
(
,
g
,
)
휕푔
for degradation and
휕퐺
(
,
g
,
)
휕푡
=
푖푗
(
)
휕퐺
(
,
g
,
)
휕푔
for conversion.
(S5)
These terms can be condensed using the Jacobian:
푠푖
:
=
휕퐺
(
,
g
,
)
휕푔
G
(
g
,
)
휕푡
=
̃
(
g
1
)
,
(S6)
where
̃
푖푗
=
푖푗
whenever
and
̃
푖푖
=
0
Í
푖푗
.
Autocatalysis of species
X
contributes the following first-order terms:
푑푃
(
,
x
,
)
푑푡
=
푖푖
[
(
1
)
(
,
1,
)−
(
,
x
,
)
]
.
(S7)
Catalytic production of species
X
from species
X
contributes analogous terms:
푑푃
(
,
x
,
)
푑푡
=
푖푗

푥푃
(
,
1,
)−
(
,
x
,
)

.
(S8)
By differentiating
,
,
, and
2
with respect to
and
and exploiting the linearity of the master equation, we
obtain the corresponding PDE terms:
휕퐺
(
,
g
,
)
휕푡
=
푖푖
(
1
)
휕퐺
(
,
g
,
)
휕푔
for autocatalysis and
휕퐺
(
,
g
,
)
휕푡
=
푖푗
(
1
)
휕퐺
(
,
g
,
)
휕푔
for catalysis.
(S9)
1
Gorin, Yoshida, and Pachter
These terms can also be condensed using the Jacobian:
G
(
g
,
)
휕푡
=

diag
g
̃
(
g
1
)

,
(S10)
where
̃
푖푗
=
푖푗
for all
,
.
Finally, each state
may have a fairly general set of bursty transcription processes, indexed by
. Each process has a
Poisson process arrival rate
. Each arrival event generates molecules; the number of molecules produced is drawn from an
-dimensional distribution
(
z
)
. The burst frequency
and burst size distribution
can arbitrarily vary in time.
To account for bursting, we combine probability flux contributions from all possible burst sizes
z
:
푑푃
(
,
x
,
)
푑푡
=
(
)
"
∑︁
z
(
z
,
)
(
,
x
z
,
)−
(
,
x
,
)
#
.
(S11)
By applying the convolution theorem, we find that each synthesis pathway contributes the following PDE term:
휕퐺
(
,
g
,
)
휕푡
=
(
1
)
(
,
g
,
)
,
(S12)
where
is the generating function of
(
z
,
)
, computed through Equation 3. For compactness, it is typically easier to aggregate
these reaction pathways into a vector function
A(
g
,
)
. Each entry of
A
, indexed by
, contains the following sum:
∑︁
(
)
(
(
g
,
)−
1
)
.
(S13)
In other words, transcriptional bursting contributes the following PDE terms:
G
(
g
,
)
휕푡
=
G
⊙A(
g
,
)
,
(S14)
where
denotes the Hadamard product. Aggregating Equations S2, S6, S10, and S14, we obtain the following system of partial
differential equations:
G
(
g
,
)
휕푡
=
(
)
T
G
+
G
⊙A(
g
,
)+
̃
(
g
1
)+

diag
g
̃
(
g
1
)

.
(S15)
It is typically more practical to use the variable
u
:
=
g
1
:
G
(
u
,
)
휕푡
=
(
)
T
G
+
G
⊙A(
u
,
)+
̃
u
+

diag
(
u
+
1
)
̃
u

=
(
)
T
G
+
G
⊙A(
u
,
)+
̃
u
+

diag
u
̃
u
+
̃
u

=
(
)
T
G
+
G
⊙A(
u
,
)+
[
u
+
diag
u
u
]
,
(S16)
where
:
=
̃
+
̃
and
:
=
̃
.
It remains to solve the partial differential equation system using the method of characteristics (John, 1978, Singh and Bokes,
2012). First, we define the
H
operator:
H(
u
,
)
G
:
=
(
)
T
G
G
⊙A(
u
,
)
,
i.e.,
G
(
u
,
)
휕푡
=
−H(
u
,
)
G
+
[
u
+
diag
u
u
]
,
H(
u
,
)
G
=
G
(
u
,
)
휕푡
+
[
u
+
diag
u
u
]
.
(S17)
Next, we take a total derivative of each component of
G
with respect to
s
, the characteristic variable. For each state, indexed by
, we yield
푑퐺
s
=
휕퐺
휕푇
푑푇
s
+
∑︁
=
1
휕퐺
휕푈
푑푈
s
.
(S18)
2
Assessing Markovian and delay models for snRNA-seq
This equation encodes the behavior of
characteristic curves
emanating from some point
(
,
u
)
. The derivative terms must match
those obtained from Equation S17:
(
H(
u
,
)
G
)
=
휕퐺
(
,
u
,
)
휕푡
+
∑︁
=
1
휕퐺
(
,
u
,
)
휕푢
(
u
+
diag
u
u
)
.
(S19)
The
curve, which propagates the time, satisfies
푑푇
s
=
1
such that
(
s
=
0
)
=
,
(S20)
yielding
(
s
)
=
s
. The
curve, which propagates the molecule amounts, satisfies
푑푈
s
=
(
U
(
s
)+
diag
U
(
s
)
U
(
s
)
)
such that
(
s
=
0
)
=
.
(S21)
These curves are generally coupled in a nontrivial way, and the system in Equation S21 needs to be solved simultaneously by
eigendecomposition. By matching the rest of the terms, we obtain the left-hand side
G
s
=
H(
U
(
u
,
s
)
,
(
s
))
G
.
(S22)
In other words, to solve the PDE system for
G
, we can first obtain the functional form of
U
, then use Equation 2 to propagate it
under a particular upstream driver. This equation is typically not analytically tractable, but straightforward to evaluate using
standard numerical schema.
The term-matching procedure relies on the Jacobian
being
left
-multiplied in Equation S17. This construction allows
us to exploit the fact that the
- and
-dependent terms in Equation S19 are identical for all
, and produce a common set
of characteristics
U
for all gene states. On the other hand, feedback regulation — for example, a molecule catalyzing one or
another state switching transition — is conceptually analogous to
right
-multiplication, which prevents such simplification
(Gorin et al., 2023). A comparably principled and numerically facile way of evaluating regulated systems does not exist as of
yet. Similarly, this procedure does not appear to generalize to multi-molecular reactions.
S1.2 Catalysis
S1.2.1 Derivation
The derivation of the characteristics for the system in Equation 5 requires instantiating the simpler system
X
1
=
⇒∅
,
X
1
−→X
1
+X
2
,
X
2
−→∅
.
(S23)
Evidently, if there is a single molecule of
X
2
at
=
0
, it vanishes after an exponentially distributed delay with rate
, i.e.,
2
=
2
s
.
If there is a single molecule of
X
1
at
=
0
, this molecule remains in the system with probability 1 until
=
. In fact, the
dynamics of
X
2
behave identically to the following non-catalytic system with a time-dependent transcription rate:
(
)
−−−→X
2
−→∅
,
(S24)
where
(
)
=
I
(
푡 < 휏
)
. From standard identities (Gorin and Pachter, 2022, Gorin et al., 2023, Jahnke and Huisinga, 2006), we
find:
log
(
2
,
)
=
0
I
(
s
< 휏
)
2
s
s
=
푞푢
2
0
I
(
s
< 휏
)
(
s
)
s
=
푞푢
2
0
s
s
=
푞푢
2
(
1
훾푡
)
whenever
푡 < 휏
and
=
푞푢
2
s
s
=
푞푢
2
(
(
)
훾푡
)
otherwise.
(S25)
3
Gorin, Yoshida, and Pachter
The second equality follows from the symmetry of convolutions. This function can be summarized as
2
2
(
)
, with
2
(
)
defined as
2
(
)
=

I
(
푡 < 휏
)+
I
(
푡 > 휏
)
(
)
훾푡

.
(S26)
The function
2
is the mean of a Poisson distribution. It remains to note that, by definition, the generating function of this
system is
(
1
,
2
,
)
=
1
∑︁
1
=
0
∑︁
2
=
0
(
1
+
1
)
1
(
2
+
1
)
2
(
1
,
2
,
)
=
"
(
1
+
1
)
0
∑︁
2
=
0
(
2
+
1
)
2
(
1
=
0,
2
,
)+(
1
+
1
)
1
∑︁
2
=
0
(
2
+
1
)
2
(
1
=
1,
2
,
)
#
=
"
∑︁
2
=
0
(
2
+
1
)
2
(
1
=
0,
2
,
)+(
1
+
1
)
∑︁
2
=
0
(
2
+
1
)
2
(
1
=
1,
2
,
)
#
=
h
I
(
)
2
2
+(
1
+
1
)
I
(
푡 < 휏
)
2
2
i
=
[
I
(
)+(
1
+
1
)
I
(
푡 < 휏
)
]
2
2
.
(S27)
Subtracting unity per Equation 4 and using
s
as the integration variable, we yield the form of
1
given in Equation 6.
S1.2.2 Simulation and evaluation
As the system in Equation 5 is fairly simple, we did not implement a Gillespie-like stochastic simulation algorithm to generate
realizations. Instead, we used standard properties of Poisson point processes.
We generated the number of
X
1
arrival events on
[
0,
]
by sampling from a Poisson distribution with rate
훼푇
, then generated
the
X
1
arrival times from a uniform distribution on
[
0,
]
. Next, we iterated over the molecules of
X
1
, and generated the
X
2
arrival event counts and times by the same procedure, over
[
,
+
]
, where
is the arrival time for the molecule of
X
1
under
consideration. Finally, we generated
X
2
lifetimes by drawing random variates from an exponential distribution with rate
. We
evaluated the number of extant molecules at a grid of 100 uniformly spaced time points to obtain the time-dependent molecule
distribution. We performed 10,000 independent simulations of this system, using
=
2
,
=
1.2
,
=
1
,
=
3
, and
=
1.4
.
To evaluate the theoretical probabilities, we constructed a grid of spectral function values:
=
2
푖휋 푗
1
for
=
0,
...
,
1.
(S28)
To calculate the
X
1
marginal at time
, we calculated a vector
,
=
(
1
=
,
2
=
0,
)
by directly integrating Equation
S27 using order-60 Gaussian quadrature, implemented through the
SciPy
function
integrate.fixed_quad
(Virtanen et al.,
2020). To define the state space size
픰 = 픰
1
, we used the maximum
X
1
count observed at
, then added 30. To obtain the
probability distribution, we used Fourier inversion, implemented through the
SciPy
function
fft.ifftn
(Virtanen et al., 2020).
To calculate the
X
2
marginal, we used an identical procedure, but using the appropriate
픰 = 픰
2
and
,
=
(
1
=
0,
2
=
,
)
.
This procedure directly follows Bokes et al. (2012), Singh and Bokes (2012). The results are shown in Figure S1a.
S1.3 Competing Markovian reactions
In this section, we apply the derivation in Equation 12 to a species
X
that can be degraded or converted to multiple species
X
,
all in a Markovian fashion. As expected, the exit time distribution is exponential:
:
=
∑︁
=
0
푖푗
(
e
,
s
)
=
s
(S29)
Next, we need to construct an expression for the probability of having a single molecule of
X
at time
:
(
e
,
;
e
, 0
)
=
0
f
(
)
∑︁
=
1
(
e
,
;
e
,
)
푑푡
.
(S30)
4
Assessing Markovian and delay models for snRNA-seq
This equality follows from adding the probabilities of disjoint events. Here, a subtlety emerges:
f
is the probability density for
any
reaction firing, with rate
;
is the probability (
푖푗
/
) of
X
being converted to
X
. Interchanging the operators, we
obtain
(
e
,
;
e
, 0
)
=
∑︁
=
1
푖푗
0
(
e
,
)
(
s
)
푑푡
=
∑︁
=
1
0
(
e
,
)
푖푗
(
s
)
푑푡
,
(
e
,
s
)
=
∑︁
=
1
s
0
(
e
,
)
푖푗
(
s
)
푑푡
.
(S31)
An application of the Leibniz integral rule yields:
푑푈
(
e
,
s
)
s
=
s
∑︁
=
1
s
0
(
e
,
)
푖푗
(
s
)
푑푡
=
∑︁
=
1
s
s
0
(
e
,
)
푖푗
(
s
)
푑푡
=
∑︁
=
1
푖푗
(
e
,
s
)−
∑︁
=
1
s
0
(
e
,
)
푖푗
(
s
)
푑푡
=
∑︁
=
1
푖푗
(
e
,
s
)−
(
e
,
s
)
.
(S32)
This approach does not appear to be directly applicable to non-Markovian reactions, as it relies on Gillespie’s formulation
wherein the reaction index and timing are independent (Gillespie, 1992), which allows us to define an aggregated law
f
.
S1.4 Competing non-Markovian reaction simulation
S1.4.1 Derivation
To illustrate the case of competing non-Markovian reactions outlined in Section 2.2.4, we obtained the characteristics (Equation
21). This system is parametrized by degradation pathway probabilities
=
/
Í
:
1
for deterministic removal after time
1
,
2
for deterministic removal after time
2
, and
3
for stochastically timed removal after a standard half-Cauchy-distributed
waiting time, with the law
f
(
)
=
2
1
1
+
2
.
(S33)
To solve this system, we used the usual machinery for
=
1
and constitutive transcription:
log
(
,
)
=
0
훼푈
(
,
s
)
s
=
푢훼휋
1
(
I
(
푡 < 휏
1
)+
1
I
(
1
)
)
+
푢훼휋
2
(
I
(
푡 < 휏
2
)+
2
I
(
2
)
)
+
푢훼휋
3
휋푡
2
tan
1
+
log
(
1
+
2
)
,
(S34)
where the second equality arises from directly integrating
in Equation 21. As this log-generating function is proportional to
,
it defines a time-dependent Poisson distribution.
S1.4.2 Simulation and evaluation
As in Section S1.2.2, we used a direct (non-Gillespie) simulation owing to the simplicity of the problem.
5
Gorin, Yoshida, and Pachter
Figure S1
The distributions obtained by direct simulation of delayed systems agree with analytical results.
a.
The time-dependent distributions of a catalytic system with two species agree with numerically computed marginals
(histograms: simulated data; dashed red lines: distributions computed through generating function; rows: time points; columns:
X
1
and
X
2
marginals).
b.
The time-dependent distribution of a system with three non-Markovian decay pathways is Poisson, and agrees with the
analytical solution (gray line: simulation mean; red line: simulation variance; dashed blue line: analytical solution for the
process mean).
From the standard properties of Poisson point processes, we generated the number of arrival events on
[
0,
]
by sampling
from a Poisson distribution with rate
훼푇
, then generated the arrival times from a uniform distribution on
[
0,
]
. We generated
the pathway index from a categorical distribution over
and generated the molecule lifetimes accordingly. Finally, we evaluated
the number of extant molecules at a grid of 100 uniformly spaced time points to obtain the time-dependent molecule distribution.
We performed 10,000 independent simulations of this system, with
=
2
,
1
=
0.5
,
2
=
1
,
=
1.2
, and
=
[
0.5, 0.4, 0.1
]
.
As the distribution shape (Poisson) is so simple, we did not plot histograms. Instead, we plotted the mean and variance
of the simulations to verify that they are approximately equal, as well as to demonstrate that the process does not appear to
converge to a steady state. In addition, we plotted the analytical expression for the Poisson mean and variance, i.e., Equation
S34 evaluated at
=
1
. The results are shown in Figure S1b.
S1.5 Solving an arbitrary multi-state system
S1.5.1 Derivation
With the generic strategy described in Section 2.3, we sought to demonstrate its applicability to fairly complex systems with
multiple steps of processing, several non-Markovian reactions, and a multi-state promoter. We constructed the system shown in
Figure S2a, with four promoter states and four species. The promoter is characterized by the transition matrix:
=
00
01
02
0
10
11
0
13
20
0
22
0
0
31
0
33
,
(S35)
6
Assessing Markovian and delay models for snRNA-seq
where
푠푠
:
=
Í
푠푗
. This is a standard stochastic matrix. The transcription process is defined by the following operator:
A(
u
)
=
0,1
1
1,2
2
0
0
.
(S36)
In other words, promoter state 0 produces species
X
1
and promoter state 1 produces species
X
2
, both according to a Poisson
counting process. The operators
and
A
are time-independent.
Next, we need to construct the system characteristics. As outlined in Figure S2a, the conversion of
X
1
requires a deterministic
delay
1
and the degradation of
X
4
requires a deterministic delay
4
. All other interconversion and degradation, shown using
thin black lines, are Markovian. Their rates are
푖푗
for source
and destination
, where
=
0
indicates degradation. Evidently,
the terminal species’ characteristics take a simple, survival function-like form:
4
=
4
I
(
s
< 휏
4
)
,
3
=
3
30
s
,
(S37)
where the first equality applies Equation 12 and the second equality applies Equation 13 for terminal species. Next, we obtain
the characteristic for
X
2
by appropriately aggregating the above characteristics:
2
=
2
2
s
+
3
23
2
30

30
s
2
s

+
4
24
2
h
I
(
s
< 휏
4
)+
I
(
s
> 휏
4
)
2
(
s
4
)
2
s
i
,
(S38)
where
2
:
=
23
+
24
is the total efflux rate for
X
2
. This equality arises from directly applying and solving Equation 18. Finally,
we obtain the characteristic for
X
2
:
1
=
1
I
(
s
< 휏
1
)
+
2
I
(
s
> 휏
1
)
2
(
s
1
)
+
3
23
2
30
I
(
s
> 휏
1
)
h
s
(
30
1
)
2
(
s
1
)
i
+
4
24
2
h
I
(
s
∈ (
1
,
4
))−
I
(
s
> 휏
1
)
2
(
s
1
)
+
I
(
s
> 휏
1
+
4
)
2
(
s
1
4
)
i
.
(S39)
This equality follows immediately from applying Equation 12 for an intermediate species.
Although all of these characteristics are relatively straightforward to obtain manually, we used the MATLAB Symbolic
Toolbox (MathWorks, 2022a,b) to automate their computation, and rewrote the results (output using Heaviside functions) in
terms of indicator functions.
Given the vector
U
(
u
,
s
)
, it is straightforward to propagate an arbitrary initial condition
P
0
. Specifically, this initial condition
can be converted to a vector generating function
G
0
(
u
)
, where each entry
0
is the generating function of the molecule
distribution at
=
0
conditional on being in state
, scaled so
0
(
u
=
0
)
=
is the probability of starting in state
. This
identity follows immediately from the definition of the generating function (Equation 3). Thus, to encode a particular initial
condition, we compute the joint distribution:
0
(
u
)
=
Ö
=
1
0
,
(
)
,
(S40)
where
0
,
encodes either a deterministic initial condition with
,
molecules (
=
(
+
1
)
,
) or a Poisson initial condition with
,
molecules on average (
=
,
). To propagate the initial condition, we compute the vector
G
0
(
U
(
u
,
))
.
The full generating function at
is given by the solution to the ordinary differential equation system in Equation 2, using
G
0
(
U
(
u
,
))
as the initial condition at
s
=
and integrating backwards to
s
=
0
.
7