Beyond the chemical master equation: Stochastic chemical kinetics coupled with auxiliary processes

The chemical master equation and its continuum approximations are indispensable tools in the modeling of chemical reaction networks. These are routinely used to capture complex nonlinear phenomena such as multimodality as well as transient events such as first-passage times, that accurately characterise a plethora of biological and chemical processes. However, some mechanisms, such as heterogeneous cellular growth or phenotypic selection at the population level, cannot be represented by the master equation and thus have been tackled separately. In this work, we propose a unifying framework that augments the chemical master equation to capture such auxiliary dynamics, and we develop and analyse a numerical solver that accurately simulates the system dynamics. We showcase these contributions by casting a diverse array of examples from the literature within this framework and applying the solver to both match and extend previous studies. Analytical calculations performed for each example validate our numerical results and benchmark the solver implementation.


Finite volume scheme
We employ the second-order, central, flux-limited finite-volume scheme introduced in Kurganov and Tadmor [1] applied for each reaction i ∈ I one-dimensionally in the e i direction. For the sake of concreteness, we demonstrate the scheme for a system comprising a single reaction in the direction e with rate r and no non-local effects or additional discrete states, which we write as ∂ ∂t p(x, t) = − ∂ ∂e [r(x, t)p(x, t)] + 1 2Ω The directional derivative highlights the one-dimensional nature of the differential terms. We write the scheme in a semi-discrete manner: we discretise the state operators while leaving the time derivative on the left-hand side to be dealt with by a separate timemarching scheme, namely an explicit Runge-Kutta method (or multi-level method) as discussed in Kurganov and Tadmor [1]. The scheme has the form ∂ ∂t p j (t) = − 1 ∆x r j (t) p j (t) + (∇ e p) j (t) ∆x 2 − r j−e (t) p j−e (t) + (∇ e p) j−e (t) ∆x 2 1/15 where the flux limiting modulates the gradient discretisation via (∇ e p) j = minmod θ p j − p j−e ∆x , p j+e − p j−e 2∆x , θ p j+e − p j ∆x , for minmod(a, b) = 1 2 (sgn(a) + sgn(b)) min(|a|, |b|), minmod(a, b, c) = minmod(minmod(a, b), c).
The flux limiting ensures no spurious oscillations are introduced into the numerical solution. The parameter θ ∈ [1,2] influences the nature of the flux limiter: interpolating between a less oscillatory scheme (θ = 1) and a less dissipative scheme (θ = 2). In all the examples presented here we fix θ = 1, but refer the reader to Kurganov and Tadmor [1] and references therein for further discussion. Zero normal flux at the boundaries (for each reaction) ensures conservation, and is achieved by neglecting flux terms whose support would extend beyond the grid. For the conservation equations studied in Kurganov and Tadmor [1], a sufficiently small time step ∆t n obeying the Courant-Friedrichs-Lewy (CFL) condition derived in Kurganov and Tadmor [1], the scheme is monotone and satisfies a maximum principle whereby the ∞ -norm of the discrete values is non-increasing (see Theorem 5.1, Corollaries 5.1 and 5.2 in Kurganov and Tadmor [1]).
The setting of conservation equations is more restrictive than our model (5), since we incorporate diffusion, non-local effects, discrete transitions, as well as state-dependent flux functions. Therefore, the associated CFL condition is not directly applicable. Moreover, we cannot expect to obtain a maximum principle, since the continuum dynamics do not respect a maximum principle: consider, for example, any reaction network whose stationary distribution is not uniform, and take initial conditions more spread than the stationary distribution. Nonetheless, the proof of Theorem 5.1 in Kurganov and Tadmor [1] is adapted in Section 1.3 to account for these additional contributions to provide a monotonicity condition: for a small enough time step ∆t n artificial oscillations will not be introduced into the numerical solution.

Discretising the burst operator
The non-local contributions on the second line of (5) are made up of probability mass lost to non-local effects, and the integral term describing mass gained due to these effects emanating from other states and coming into x. In this section, we describe the discretisation for the class of non-local kernels describing production in bursts. A discussion on the discretisation of fragmentation kernels appears in Section 4, which builds on the present analysis.
To simplify the notation, we consider a single bursting process, and thus remove the subscripts j and k, taking Burst kernels take the form where b is the mean burst size, and the burst kernel has an unbounded support. On the truncated state space, we respect the finite boundary by retaining only those bursts that 2/15 transfer probability mass between two states within the finite domain. This amounts to changing the domain of integration to x − ze ∈ [0, N 1 ∆x] × · · · × [0, N d ∆x]. To ensure conservation, we must account for the corresponding reduction in mass loss. This is achieved in our scheme by focusing on the bursts emanating from each x (as opposed to the integral representation in (S3), which describes the bursts into the state x). For each j∆x on the grid, we find the largest burst that remains within the domain: The scheme ensures that mass is lost and gained equally. The losses are balanced by the gains: where we have denoted the probability of a jump of me grid points by P m .
To determine the jump probability P m , we recall that in the finite-volume framework, the value at a point represents the probability mass of a cell around at that point. Retaining a one-dimensional projection, The probability to burst into a cell me units away, is given by where we assume a uniform distribution in the cell of origin. For jump kernels of the form (S4), this probability may be calculated analytically, namely This discretisation preserves the mean burst size (on an infinite domain). To see this, note that the mean burst size may be calculated via Since q < 1, we may evaluate the series via (S12) From (S11) and (S12) it follows that (S13)

Monotonicity condition
In this appendix, we outline the proof in Kurganov and Tadmor [1] that the numerical scheme satisfies a maximum principle, and show how the argument may be adapted to system (5). The equation under consideration is where x ∈ R is a scalar. We denote the values of the numerical scheme at time step n by p n j , where j indexes the finite volume. The CFL condition is where ∆t n denotes the size of the nth time step, and ∆x denotes the step size in state space. From Theorem 5.1 in Kurganov and Tadmor [1] it holds that, if the CFL condition (S15) is satisfied, then the ∞ -norm of the scheme is non-increasing with a first-order Euler time step, that is, Corollaries 5.1 and 5.2 in Kurganov and Tadmor [1] extend this result to higher-order time-stepping methods.
Our first observation is that, in preserving a uniform grid, the advection and diffusion operators may be projected onto one dimension, while the CFL condition (S15) was derived for two-dimensional advection. Restriction to one dimension allows us to relax the CFL condition (S15) by a factor of 2. To understand this, we highlight that the proof expresses the explicit Euler time step as where C i,k are functions satisfying a local maximum principle. The discretisation admits a decomposition of the p n j contribution of the form which allows the right-hand side of (S17) to be expressed as a linear combination of C i,k terms. The goal is to prove that the right-hand side may be expressed as a convex combination of C i,k terms, that is, where the coefficient of each C i,k term is positive and the sum of all coefficients is one. It may be seen from the formulae for the coefficients shown in Kurganov and Tadmor [1] that they sum to one, and we focus on demonstrating their positivity. For each dimension i, it may be shown that α i,1 and α i,2 are non-negative by construction of the scheme, while the other two terms are bounded by With no control over the sign of these α i,k terms for k = 3, 4, we may guarantee that the coefficients remain positive by choosing ∆t n small enough to ensure that the p n j contribution dominates the (possibly) negative contribution via From (S20) we see that the original CFL condition (S15) from Kurganov and Tadmor [1] may be relaxed by a factor of 2, even when system (5) is multidimensional, since the differential operators act on a single dimension, that is, d is effectively one in our discretisation. We may extend the equation under consideration to incorporate the time-and statevarying fluxes, discrete states, and non-local effects described in system (5). In this general case, we cannot hope to obtain a maximum principle since the flux function variations in state can act to concentrate probability mass, while discrete transition rates may be negative resulting in exponential growth (as in the phenotypic selection case study). Therefore, the continuous dynamics do not respect a maximum principle. Instead, we focus on monotonicity. A scheme is called monotone if it may be expressed in the form Artificial oscillations driving numerical instability may result from insufficiently small time steps which violate scheme monotonicity. In analogy with the previous argument, our aim is to express the explicit Euler scheme in a form akin to (S17), and find a time step ∆t n small enough to guarantee that the coefficients of C i,k remain positive. This suffices to ensure monotonicity since the local maximum principle satisfied by the C i,k terms guarantees that these terms may be expressed as convex combinations of local p n j terms.
Neglecting terms that are inhomogeneous in the discretised system (such as the integral terms whose non-local contributions come from states of indexĵ = j and sources from discrete states of indexk = k), we rewrite (5) as where we write the differential operators using directional derivatives to emphasise their one-dimensional nature, and the rate h k encapsulates the homogeneous contributions of discrete transitions and non-local effects, namely Focusing on a single value of k, we henceforth drop the subscript (but recycle the index k for separate use). The explicit Euler scheme may be written as The coefficients α i,1 and α i,2 are non-negative, while the remaining coefficients are bounded (see (S19)) via Collecting terms of p n j , we see that for non-negative constants β j ≥ 0. Decomposing p n j in analogy to (S18), via we may write Guaranteeing that the coefficients of the C i,k terms are non-negative requires a sufficiently small ∆t n . In light of the bound (S25), it suffices for ∆t n to satisfy, for example, Condition (S29) guarantees that, for each discrete state, the scheme is monotone with respect to the homogeneous terms. The terms in the square brackets in (S29) have direct interpretations: the first term represents the zeroth order exponential contributions of h from (S23) and the second term represents the contributions from the higher order differential terms. This latter term comprises two contributions in the round brackets, the first comes from the second-order term and the second from the first-order terms.

Code example
In this appendix, we illustrate the Python implementation of the reaction networks (6) and (7), which comprise reactions, non-local bursting, and discrete states. Each network is reproduced alongside the code required to simulate the evolution of the distribution in Table A. The implementation of each network in the Flips software is a straightforward encoding of the network diagram: mRNA molecules are represented by the label 'm', protein molecules by the label 'p'. A solver object is created with the network structure, reactions are separated into those governing continuum species and those governing discrete states (for which a finite truncation is specified). Each reaction includes three pieces of information: the reactants, the products, and the reaction rate. Additionally, since we follow models in the limit as Ω → ∞, the diffusion is eliminated by setting diffusion=0. Finally, an initial distribution is set, the distribution is evolved until some time, and the result plotted.
The interface is a direct translation of the network diagram, and therefore accessible to the non-technical user. Simultaneously, options for the advanced user allow fine tuning of the state and time discretisation (see [2]).

6/15
(i) Reaction network (6):  (7):  Table A. Python code example for the Flips solver. Code to simulate the reaction networks (6) and (7) from the self-regulated gene expression case study.

Hermite equation
In this appendix, we seek separable solutions of equation (17), which we reproduce here: Writing q = T (t)X(x), and substituting into (S30) it follows that The left-hand side is a function only of t, while the right-hand side is a function only of x, therefore these must both be constant, taking the value of the growth rate, which we will denote r. This is sometimes called the Malthus parameter [3]. In other words, Equation (S32) is of Sturm-Liouville form, and does not have, in general, a closed-form solution. We are not aware of a closed-form solution for birth and death functions (20) even in the linear-growth case of G(x) = gx. Nonetheless, we may seek an asymptotic approximation in the limit as Ω → ∞. As explored in Lunz [4], the leading-order outer solution concentrates the probability mass at the critical point x c . There was no growth term in Lunz [4], however, while this influences the growth along the characteristics, the characteristics themselves are not affected and thus the same singular evolution unfolds in our case. This singular perturbation is regularised within a boundary layer in the vicinity of the critical point, by considering the scaling The leading-order inner form of (S32) is then given by where we assume that x c is a simple root of λ − µ, primes denoting differentiation with respect to the argument, and the constant is given by . whereby We now demonstrate that admissible solutions of the Hermite equation (S37) may be found if is a non-negative integer. We seek a formal series solution of the form Upon substituting (S38) into (S37) and matching coefficients of the powers of ξ, we obtain the recurrence relation The recurrence relation (S39) expresses a dependence of coefficients on those indexed of the same parity. Therefore, each coefficient c i may be expressed as a function of either c 0 or c 1 , for i even or odd, respectively. The values for c 0 and c 1 may be determined by conditions at, say, the origin: We argue, on physical grounds, that we expect the leading-order inner solution to be symmetric about the origin, whereby Y (0) = 0 = c 1 , and thus all odd coefficients vanish. For a non-negative integer ≥ 0, we see from (S39) that c = 0, from which it follows that c +2j = 0 for all j ≥ 0. Thus the series terminates at th order to give a polynomial. These polynomials are called Hermite polynomials. Further properties of the Hermite polynomials and connections to other special functions are discussed in Arfken et al. [5]. For many applications, this family of polynomial solutions is the only admissible solution, as other values of lead to solutions that are not polynomially bounded. However, in our case the justification must be slightly stronger, since from (S36) we see that Y can grow in the far field and X may still vanish. Since the outer solution concentrates probability mass in the vicinity of the stable equilibrium, we expect it to be zero to all orders at 8/15 asymptotically large time. On this basis, we now demonstrate that solutions for values of other than the non-negative integers are inadmissible.
Consider the series solution in Y for 0 < ∈ Z, and denote the smallest odd integer larger than by 2L − 1. Since only even indices c i are nonzero, it follows that We bound d i by noting from (S39) and (S41) that, for i > L, where C L = (2L − 1)!!|c 0 |/2 L and is independent of i. We thus deduce that for constants D i . Substituting the bound (S43) into (S36), we find that for other constantsD i . The second term on the right-hand side of (S44) is not beyond all orders for any L, therefore, such a solution X cannot match the outer solution and we neglect it. For < 0, Y diverges more rapidly in the far field, since Y (ξ; ) is monotonically increasing with respect to ∈ {0 ≤ n ∈ Z}. In fact, for = −1 it turns out that Y (ξ) = e ξ 2 , and X is constant. In summary, only solutions for non-negative integer are admissible. The growth rates r corresponding to mode are given by The sequence r is monotonically decreasing (since at a stable equilibrium λ − µ < 0), therefore, for large times, the = 0 mode will be dominant, dictating both the functional form and the growth rate. We thus deduce that the large-time behaviour will be described by constant Y (ξ), and thus where C is a normalisation constant, given by C = 1/ √ aΩπ for X to be a probability density. The profile is a Gaussian centred at x = x c of width O(1/ √ aΩ). The growth rate is approximated by r 0 , namely The last equality is for the specific forms G(x) = gx and the functions in (20). Defining the shorthand and integrating (18), we find that that is, the growth rate is given by the G-moment of the q density. For the separable solution q = T X = e r0t X(x), we see that from which it follows that the G-moment of p 0 is precisely r 0 .

Fragmentation models
In this appendix, we discuss analytical and numerical aspects of the growth-fragmentation equation (24), which we rewrite here: First, we demonstrate that ensuring the fragmentation splits one cell into two while conserving volume is achieved by imposing the constraints [3] 2B(x, t) = respectively. Consistency further requires that the rate (per unit size) of cell division for cells of size y into cells of size x and y − x is equal, that is, From (S51), we find that the rate of change of the number of cells in the system is governed by 10/15 The first term on the right-hand side of (S54) vanishes due to the no-flux and vanishing farfield conditions. By changing the order of integration and using the first condition (S52), we find that The rate of change of volume in the system is governed by By changing the order of integration in the final term on the right-hand side of (S56) and using the second condition (S52), we find that the second and third terms cancel, thus where we again used the no-flux and vanishing far-field conditions. In fact, here we need The calculations confirm the physical features of the model: the rate of change of the number of cells (S55) is given by the average rate of fragmentation B, but not explicitly dependent of the growth g, while the rate of change of volume (S57) is given by the average growth rate g, but not explicitly dependent of fragmentation.
Given a kernel b, we may derive B via one of the conditions in (S52), however, the other condition in (S52) as well as condition (S53) still constrain b in a highly non-trivial way. With the aim of simplifying the fragmentation description, as well as bringing it in line with the framework based on system (5), we introduce a new representation of the fragmentation function b (from which we will derive B) that will decouple the constraints without loss of generality. We define b(y, x, t) = 2r(y, t)f y (x, t), where r is the fragmentation rate of cells of size y, and f y is the probability density that the fragmentation from size y produces, for a given daughter cell, a cell of size x. We define the corresponding cumulative probability function F y via and note that F y (y, t) = 1.
We satisfy the symmetry condition (S53) by requiring that We may satisfy conditions (S52) by requiring that where the last equality comes from integrating by parts and using (S60). Upon substituting (S62) into (S63), the constraint takes the form We proceed to show how constraint (S64) is satisfied by the representation (S58). From the probability property (S60) and the symmetry condition (S61) we see that To recap, by introducing the representation (S58) for b based on a decomposition into a rate and a symmetric probability density, we trivially satisfy both conditions in (S52) and condition (S53). The fragmentation rate B is equivalent to the rate r in this representation (and to the rate f in the ACME formulation (5)). Importantly, this representation expresses the kernel b in terms of identifiable quantities. In particular, the practitioner needs to specify the rate r at which fragmentation occurs for each cell size, and the distribution of fragments f y . The kernel B in the ACME formulation (5) is given by 2f y , and is thus twice a probability density.
There are additional advantages to this representation. This formulation allows for a straightforward description of self-similar fragmentation kernels, that is, the dependence of f y on y is only a rescaling, namely where f (z, t) is a density on z ∈ (0, 1). Symmetric and asymmetric cell division, for example, is commonly modeled with size distributions in proportion with the original cell [3], in line with this description. Moreover, it is easy to incorporate Dirac delta functions within f y since the finite-volume scheme is built via integrating over finite regions of the state space. Therefore, the cumulative representation of f y via F y suffices, and this is a well-defined function. To see this, we detail the discretisation of the fragmentation operator.
As with bursts, for each volume centred at i∆x we subtract the propensity of fragmenting to within the volume centred at j∆x for 0 ≤ j < i and add that propensity to the j volume, thereby achieving conservation of volume. Note that in this paragraph we have used the word "volume" to refer to what is normally called a "cell", so as not to confuse the use of the word "cell" elsewhere in this paper referring to the biological entity. The probability of fragmenting from j to i is given by The probability (S68) is expressed via its cumulative density. Further simplification is possible, for example, in the self-similar case, If the distribution is only Dirac masses, for example, time-invariant asymmetric cell division at size fractions c and 1 − c for c ∈ (0, 1), then f (x, t) = [δ(x − c) + δ(x − 1 + c)]/2, in which case (S69) further simplifies to where 1 x>0 is the Heaviside function and Q is the integral Q(c) = 1 0 1 x<c dx = min(1, max(0, c)).

(S71)
The formula (S70) accelerates construction of the scheme since no integrals need to be computed, and demonstrates how Dirac delta functions are accurately captured in the finite volume discretisation. The fragmentation probability is multiplied by the fragmentation rate r to give the local fragmentation propensity, and probability mass is transferred from volume j in proportion to the propensity, to volume i with twice this proportion to account for the cell creation (this is the factor of 2 in (S58)).

Explicit growth-fragmentation solutions
In this appendix, we compile analytical solutions of growth-fragmentation equations for benchmark use. We begin by considering the growth-free problem considered in Rooney et al. [6], namely ∂ ∂t ρ(x, t) = ∞ x 2ay k−1 ρ(y, t) dy − ax k ρ(x, t).

13/15
We refer the reader to Rooney et al. [6] for a step-by-step solution of (S72). In brief, we seek a similarity solution of the form ρ(x, t) = x α h(η), η = xt β .
It can be shown that the scalings must satisfy α = −2 and β = 1/k. The equation (S72) may then be reduced to an ODE for h, admitting the solution where c is an arbitrary constant. We now employ a change of variables presented in Cáceres et al. [7] which transforms the solution of the general growth-free equation such as (S72), to the solution of a growth-fragmentation equation with non-uniform growth and uniform degradation. We assume that the fragmentation kernel is timeinvariant homogeneous of degree γ − 1 for γ > 0, that is, b(sy, sx) = s γ−1 b(x, y), for any s > 0. For equation (S72) the kernel is homogeneous with γ = k. We now defineρ (x, t) = e −2t ρ(e −t x, e γt − 1), and look for an equation governing the time evolution ofρ by taking partial derivatives. First with respect to x, we find that ∂ ∂xρ (x, t) = e −3t ∂ρ ∂x (e −t x, e kt − 1).
Then, using the dynamics (S72), and (S78), we see that Upon substituting (S77), changing variables, and using homogeneity assumption (S75), the square brackets on the right-hand side of (S79) may be expressed as whereby we see thatρ satisfies the growth-fragmentation equation 14/15