Incorporating age and delay into models for biophysical systems

It has been reported that some biological processes do not take place instantaneously. In other words, there is a time lag between the initiation and the completion of the process. Time delays are observed inherently in many biological systems, such as gene transcription [9–11] and translation [12], cell cycle in cancer treatment [13], intracellular viral dynamics [14, 15], control of infectious diseases [16], population growth [17, 18], RNA and protein folding [19, 20], and enzyme catalyzed reactions [21, 22]. Sometimes time delays are introduced purposefully as a useful means to reduce model complexity and to compensate for the lack of experimental observation in both deterministic and stochastic models of biological processes.

Intermediate, ancillary processes or unobserved reactions can be replaced by time delays. For example, production of hes1 mRNA from hes1 gene has been modeled using delay differential equations where detailed mRNA synthesis and processing steps are replaced by a time delayed reaction [23]. While modeling the mammalian circadian clock, intermediate protein dynamics can be simplified as transcriptional feedback loops with time delayed variables in delay differential equations [24]. In enzyme catalyzed reactions with multiple intermediates, the production of the final product can be expressed as a distributed delay equation, which is a useful tool when measurements on multiple intermediates in the experiment are not available [25].

Introduction of time delays as a model reduction technique has also been applied in discrete stochastic models for CRNs. For instance, model complexity of unimolecular reaction networks is reduced by generating delay distributions with key model features that are derived by computing first passage times of target species [26]. In [27], the production of yellow fluorescent protein has been described using a time-delayed birth and death process where a randomly distributed time delay was generated to simplify a sequence of steps in gene activation.

In most previous works in this area, the focus was on investigating stochastic models for CRNs with constant or random time delays. In those models, the probability that a reaction occurs within the next short interval of time is commonly described by a propensity (also known as intensity) function of the reaction. The waiting time for non-delayed reactions is exponentially distributed [28]. In practice, the occurrence of some reactions is not only determined by the molecular counts of the reactants but also affected by the age distributions or lifetimes of the reactant molecules. For example, mRNA decay rates vary depending on the age of each mRNA. Moreover, the age of the mRNAs determines polysome size distributions and protein synthesis rates in translation ([29, 30], chapters 3 and 5 in [31]). It was also reported that an mRNA tail length distribution depends on the average age of mRNA population and that the tail-length distribution plays a significant role in deadenylation and decay dynamics of mRNA populations [32, 33].

When time delays are used to aggregate out ancillary or unobserved processes and reduce model complexity, it makes more sense that the length of time delay depends on the age of each reactant molecule (e.g., mRNA, protein, and enzymes). Therefore, it is worthwhile to consider an individual-based age-structured stochastic model for CRNs.

In this work, we develop a way to describe CRNs with random time delays and non-delayed reaction rates incorporating the age of each reactant and making use of hazard functions in survival analysis [34, 35]. See appendix A for some preliminaries on relevant mathematical and statistical concepts. In our approach, the hazard functions are set as constant, time-dependent, or age-dependent functions generalizing the notion of reaction rate constants in propensity functions. Our model keeps track of the age of each reactant molecule and provides a new way to express time delays in non-Markovian models. Moreover, the method also allows us to describe discrete stochastic CRNs with constant or random time delays without age dependence, as considered in previous works. We study the large-volume limit of the proposed non-Markov CRN and provide a mean-field PDE limit for the age densities by virtue of the law of large numbers (LLN), as opposed to an ODE limit in the classical theory. The PDE limit is based on existing results in the literature [36, 37] and follows from the standard limit theory for measure-valued Markov processes. However, novel usage of the PDE limit can provide further approximations and pave the way for efficient simulation algorithms. For the sake of illustration, we show how the PDE limit can be used to approximate mean first passage times (MFPTs) in the context of CRNs. As another by-product of the LLN, we show how an efficient (fast) hybrid simulation algorithm can be devised when a subset of the CRN is abundantly available, giving a flavor of multiscale approximation. Finally, as simple applications of our approach, we briefly discuss a prokaryotic auto-regulation and the quasi-steady state approximation (QSSA) in the context of the Michaelis–Menten enzyme kinetic reactions. Numerical examples have been provided wherever deemed necessary. For the sake of ready usage of our methods, the Julia scripts used in the numerical examples have been made available via a GitHub repository [38].

The following notational conventions are adhered to throughout the paper. We use 1 _{A}(x) to denote the indicator (or characteristic) function of a set A, i.e., 1 _{A}(x) = 1 if and only if x ∈ A. Given a suitable space E, let D([0, ∞), E) (or D([0, T], E)) denote the space of E-valued càdlàg functions defined on [0, ∞) (or [0, T], for some T > 0). The set of Borel subsets of a set A will be denoted by $\mathcal{B}\left(A\right)$. The set of natural numbers are denoted by $\mathbb{N}$. The set of real numbers is denoted by $\mathbb{R}$. Other notations will be introduced as and when needed.

The standard way to model the CRN in equation (2.1) is to use a CTMC to describe the counts of molecules of the species A and B over time. In such a model, whenever each reaction fires, the consumption and the production of molecules are instantaneous. Let ${\tilde {X}}_{A},\;{\tilde {X}}_{B}$ denote the stochastic processes counting the copy numbers of the species A and B respectively. Here, the quantities b, τ, and d are reaction rate constants. The propensity functions corresponding to the three chemical reactions are defined as

$$\begin{align*}\hfill & {\lambda }_{b}\left(t\right)=b,\hspace{25.0pt}{\lambda }_{\tau }\left(t\right)=\tau {\times}{x}_{A}\left(t\right),\hspace{25.0pt}\hfill \\ \hfill & {\lambda }_{d}\left(t\right)=d{\times}{x}_{A}\left(t\right),\hfill \end{align*}$$

where x_i (t) denotes the number of molecules of the chemical species i at time t, for i = A, B. Define T_k to be the waiting time until the next reaction of type birth (k = b), conversion (k = τ), and death (k = d). Then, T_k is exponentially distributed with rate λ_k (t) for k = b, τ, d. The probability of each reaction's occurrence is expressed in terms of the corresponding propensity function as follow:

$$\begin{align*}\hfill & \mathsf{P}\left(t{\leqslant}{T}_{k}{< }t+{\Delta}t\vert {\tilde {X}}_{A}\left(t\right)={x}_{A},{\tilde {X}}_{B}\left(t\right)={x}_{B}\right)\hfill \\ \hfill & \hspace{25.0pt}\approx {\lambda }_{k}\left(t\right){\Delta}t+o\left({\Delta}t\right),\hfill \end{align*}$$

for k = b, τ, d when Δt is small enough. Then, the trajectory equations can be written in a straightforward fashion following the random time changed representation of Poisson processes as

$$\begin{equation*}\begin{aligned}\hfill {\tilde {X}}_{A}\left(t\right)& ={\tilde {X}}_{A}\left(0\right)+{R}_{1}\left(bt\right)-{R}_{2}\left({\int }_{0}^{t}\tau {\tilde {X}}_{A}\left(s\right)\enspace \mathrm{d}s\right)\hfill \\ \hfill & \quad -{R}_{3}\left({\int }_{0}^{t}d\hspace{2.0pt}{\tilde {X}}_{A}\left(s\right)\enspace \mathrm{d}s\right),\hfill \\ \hfill {\tilde {X}}_{B}\left(t\right)& ={R}_{2}\left({\int }_{0}^{t}\tau {\tilde {X}}_{A}\left(s\right)\enspace \mathrm{d}s\right),\hfill \end{aligned}\end{equation*}$$

where R₁, R₂, and R₃ are unit rate Poisson processes [2]. We assume we do not have any B molecules in the system initially, i.e., ${\tilde {X}}_{b}\left(0\right)=0$. Now, if we scale the stochastic processes by a scaling parameter n, e.g., volume of the system, it follows directly from the LLN for Poisson processes [41, 42] that the scaled process $\left({n}^{-1}{\tilde {X}}_{A},{n}^{-1}{\tilde {X}}_{B}\right)$ can be approximated by the solution to the following system of ODEs:

$$\begin{align*}\hfill \begin{aligned}\hfill \frac{\enspace \mathrm{d}}{\enspace \mathrm{d}t}{x}_{A}\left(t\right)& =-\left(\tau +d\right){x}_{A}\left(t\right),\hfill \\ \hfill \frac{\enspace \mathrm{d}}{\enspace \mathrm{d}t}{x}_{B}\left(t\right)& =\tau {x}_{A}\left(t\right).\hfill \end{aligned}\end{align*}$$

Notice that the birth rate b vanishes in the limit because we did not assume any scaling of b with respect to n. In general, one would assume that the overall birth rate scales linearly with n so that it is retained in the limit.

Now, let us introduce age and delay into the CRN described by equation (2.1). We assume that the production rate of B and the degradation rate of A depend on the age of the reactant molecule of species A. We use 'age' as an umbrella term to refer to the amount of elapsed time since a specific event. Thus, 'age' could mean different things depending on the application area. The most straightforward way is the biological or the physical age, which we take as the time duration since the molecule was born or created. In systems where a certain reaction can fire only when a gene is activated, one could define age as the time duration since activation of the gene. In some cases, it may be desirable to define delays in terms of time duration since the initiation of a reaction. The notion of age is sufficiently general to account for those cases as well. For example, a reaction A → B in which the delay is defined purely in terms of time difference between initiation and completion of the reaction, can be replaced by the reaction system A → F → B where F is a fictitious species. The physical age of this fictitious species F is precisely the time since the initiation of the reaction A → B. Now, putting an appropriate hazard function on the reaction F → B, we can introduce a random or a deterministic delay in the reaction A → B. Therefore, for the CRN in equation (2.1), it seems sufficient to define the age to be the physical age of the molecules of A.

When we have an age-structured model, the counts (copy numbers of the species A and B) are inherently non-Markovian unless the holding times are exponentially distributed. However, if we keep track of the ages of the molecules in addition to the counts, we can get a Markov system, albeit on a more abstract state space. A neat way to do so is to use measure-valued processes that keep track of the age distribution of the molecules over time. Moreover, the measure-valued processes are also Markovian, which allows us to make use of the already existing limit theory for Banach space-valued Markov processes. This approach to age-structured modeling in biology is not new. Our work builds on the existing literature [36, 37, 43, 44]. In the next section, we describe how the measure-valued processes can be utilized in the context of the CRN in equation (2.1).

Let us denote by N_A (t) and N_B (t) the numbers of molecules of the chemical species A and B at time t. Then, individual molecules of A are labeled 1, 2, ..., N_A (t). We denote the age of the ith molecule of the species A by a_i (t) for i = 1, 2, ..., N_A (t). Similarly, we denote by b_j (t) the age of the jth molecule of the species B at time t. Now, we define a measure-valued process ${X}_{t}=\left({X}_{t}^{A},{X}_{t}^{B}\right)$ where ${X}_{t}^{A}$ and ${X}_{t}^{B}$ describe the age distributions of chemical species A and B at time t. To be more precise, we define

$$\begin{equation}{X}_{t}^{A}{:=}\sum _{i=1}^{{N}_{A}\left(t\right)}{\delta }_{{a}_{i}\left(t\right)},\hspace{25.0pt}{X}_{t}^{B}{:=}\sum _{i=1}^{{N}_{B}\left(t\right)}{\delta }_{{b}_{i}\left(t\right)},\end{equation} \tag{ 2.3 }$$

where δ_x is the Dirac measure, a function that takes value 1 if the argument to the function (a measurable set) contains x and zero otherwise. The components ${X}_{t}^{A}$ and ${X}_{t}^{B}$ of X_t are finite point measures with atoms placed on the individual ages of the molecules. For example, ${X}_{t}^{A}\left(\left(0.5,11.25\right]\right)={\sum }_{i=1}^{{N}_{A}\left(t\right)}{\delta }_{{a}_{i}\left(t\right)}\left(\left(0.5,11.25\right]\right)$ gives us the count of species A molecules with ages in the set (0.5, 11.25] at time t. In general, ${X}_{t}^{A}\left(F\right)$ gives us the count of species A molecules whose ages lie in the set F at time t.

For any point measure $\mu ={\sum }_{i=1}^{n}{\delta }_{{x}_{i}}$ and a measurable function f, we denote the integration of the function f with respect to the measure μ by

$$\begin{equation*}\langle \mu ,f\rangle {:=}\int f\enspace d\mu =\sum _{i=1}^{n}f\left({x}_{i}\right).\end{equation*}$$

If μ := (μ₁, μ₂, ..., μ_L ), for some positive integer L, is a vector of point measures and f is a measurable function, we use the notation ⟨⟨μ, f⟩⟩ to denote

$$\begin{equation*}\langle \langle \mu ,f\rangle \rangle {:=}\sum _{i=1}^{L}\langle {\mu }_{i},f\rangle .\end{equation*}$$

Therefore, we have ${N}_{A}\left(t\right)=\langle {X}_{t}^{A},1\rangle ={X}_{t}^{A}\left({\mathbb{R}}_{+}\right)$ and ${N}_{B}\left(t\right)=\langle {X}_{t}^{B},1\rangle ={X}_{t}^{B}\left({\mathbb{R}}_{+}\right)$ where 1 stands for the identity function. The set of non-negative real numbers is denoted by ${\mathbb{R}}_{+}$. The total population size is given by

$$\begin{equation*}N\left(t\right){:=}\langle \langle {X}_{t},1\rangle \rangle ={N}_{A}\left(t\right)+{N}_{B}\left(t\right).\end{equation*}$$

The process X_t is a Markov process on the space $D\left(\left[0,T\right],{\mathcal{M}}_{P}\left({\mathbb{R}}_{+}\right){\times}{\mathcal{M}}_{P}\left({\mathbb{R}}_{+}\right)\right)$ where T > 0 is a finite time horizon and ${\mathcal{M}}_{P}\left({\mathbb{R}}_{+}\right)$ is the space of finite, point measures on ${\mathbb{R}}_{+}$.

In order to simplify notations, we introduce maps ${\sigma }_{i}:{\mathcal{M}}_{P}\left({\mathbb{R}}_{+}\right)\to {\mathbb{R}}_{+}$, for i = 1, 2, 3, ..., the purpose of which is to extract the ith atom (the age of the ith molecule) from a point measure following some partial order (e.g., 'greater or equal to' relation). Therefore, ${\sigma }_{i}\left({X}_{t}^{A}\right)$ gives us the age of the ith molecule of the species A at time t. We can now write down the trajectory equations:

$$\begin{align}\hfill {X}_{t}^{A}& =\sum _{k=1}^{{N}_{A}\left(0\right)}{\delta }_{t+{\sigma }_{k}\left({X}_{0}^{A}\right)}+{\int }_{0}^{t}{\int }_{0}^{\infty }{\delta }_{t-s}\enspace {\mathsf{1}}_{\left\{\theta {\leqslant}b\right\}}\enspace \hfill \\ \hfill & \quad {\times}{Q}_{1}\left(\mathrm{d}s,\mathrm{d}\theta \right)-{\int }_{0}^{t}{\int }_{\mathbb{N}}{\int }_{0}^{\infty }{\delta }_{t-s+{\sigma }_{i}\left({X}_{s-}^{A}\right)}\hfill \\ \hfill & \quad {\times}\enspace {\mathsf{1}}_{\left\{i{\leqslant}{N}_{A}\left(s-\right)\right\}}\enspace {\mathsf{1}}_{\left\{\theta {\leqslant}\tau \left({\sigma }_{i}\left({X}_{s-}^{A}\right)\right)\right\}}\enspace {Q}_{2}\left(\mathrm{d}s,di,\mathrm{d}\theta \right)\hfill \\ \hfill & \quad -{\int }_{0}^{t}{\int }_{\mathbb{N}}{\int }_{0}^{\infty }{\delta }_{t-s+{\sigma }_{i}\left({X}_{s-}^{A}\right)}\enspace {\mathsf{1}}_{\left\{i{\leqslant}{N}_{A}\left(s-\right)\right\}}\enspace \hfill \\ \hfill & \quad {\times}{\mathsf{1}}_{\left\{\theta {\leqslant}d\left({\sigma }_{i}\left({X}_{s-}^{A}\right)\right)\right\}}\enspace {Q}_{3}\left(\mathrm{d}s,di,\mathrm{d}\theta \right),\hfill \end{align} \tag{ 2.4 }$$

$$\begin{align}\hfill {X}_{t}^{B}& ={\int }_{0}^{t}{\int }_{\mathbb{N}}{\int }_{0}^{\infty }{\delta }_{t-s}\enspace {\mathsf{1}}_{\left\{i{\leqslant}{N}_{A}\left(s-\right)\right\}}\enspace {\mathsf{1}}_{\left\{\theta {\leqslant}\tau \left({\sigma }_{i}\left({X}_{s-}^{A}\right)\right)\right\}}\enspace \hfill \\ \hfill & \quad {\times}{Q}_{2}\left(\mathrm{d}s,di,\mathrm{d}\theta \right),\hfill \end{align} \tag{ 2.5 }$$

where Q₁, Q₂, Q₃ are independent Poisson point measures (PPMs) with intensity measures ds × dθ, ds × di × dθ, and ds × di × dθ respectively, where di is a counting measure on $\mathbb{N}$, and ds and dθ are Lebesgue measures on ${\mathbb{R}}_{+}$. Provided the global jump rates are upper bounded by a finite quantity and the initial population size does not explode ${\mathrm{sup}}_{n}\mathsf{E}\left[{n}^{-1}{N}_{A}\left(0\right)\right]{< }\infty$, the trajectory equations admit a unique pathwise solution $\left({X}_{t}^{A},{X}_{t}^{B}\right)$ (see [37, theorem 2.5] for a similar derivation).

Under some assumptions on the hazard functions and the initial age distribution of the A molecules, the scaled process n⁻¹ X_t converges to a deterministic, continuous function ${x}_{t}{:=}\left({x}_{t}^{A},{x}_{t}^{B}\right)$ whose components ${x}_{t}^{A}$ and ${x}_{t}^{B}$ are themselves measure-valued functions satisfying

$$\begin{align*}\hfill \langle {x}_{t}^{A},{f}_{t}\rangle & =\langle {x}_{0}^{A},{f}_{0}\rangle +{\int }_{0}^{t}{\int }_{0}^{\infty }\left(\frac{\partial }{\partial a}{f}_{s}\left(a\right)+\frac{\partial }{\partial s}{f}_{s}\left(a\right)\right.\hfill \\ \hfill & \quad \left.-{f}_{s}\left(a\right)\left(\tau \left(a\right)+d\left(a\right)\right)\right){x}_{s}^{A}\left(\enspace \mathrm{d}a\right)\enspace \mathrm{d}s\hfill \\ \hfill \langle {x}_{t}^{B},{f}_{t}\rangle & ={\int }_{0}^{t}{\int }_{0}^{\infty }\left(\frac{\partial }{\partial a}{f}_{s}\left(a\right)+\frac{\partial }{\partial s}{f}_{s}\left(a\right)\right.\hfill \\ \hfill & \quad \left.+{f}_{s}\left(0\right)\tau \left(a\right)\right){x}_{s}^{A}\left(\enspace \mathrm{d}a\right)\enspace \mathrm{d}s,\hfill \end{align*}$$

for a sufficiently large class of test functions f : (a, s) → f_s (a). The convergence of the scaled stochastic process n⁻¹ X_t to the deterministic function x_t can be proved using techniques similar to those in [36, 37, 43–45]. However, for the sake of completeness, a brief, intuitive argument is presented in appendix B.

Since the measure-valued function ${x}_{t}^{B}$ is determined entirely by ${x}_{t}^{A}$, it suffices to study ${x}_{t}^{A}$. The densities y_A (t, a) of the measure ${x}_{t}^{A}$, when they exist, are an important quantity describing the distribution of age of the species A molecules in the large-volume mean-field limit. The density function y_A should satisfy

$$\begin{equation}\left({\partial }_{t}+{\partial }_{s}\right){y}_{A}\left(t,s\right)=-\left(\tau \left(s\right)+d\left(s\right)\right){y}_{A}\left(t,s\right),\end{equation} \tag{ 2.6 }$$

with the initial and the boundary conditions

$$\begin{equation*}{y}_{A}\left(0,s\right)={f}_{A}\left(s\right),\hspace{25.0pt}{y}_{A}\left(t,0\right)=0,\end{equation*}$$

where f_A (s) specifies the age distribution of A molecules at time t = 0. To be more precise, it is the density of the limiting measure ${x}_{0}^{A}$, which we assume exists, with respect to the Lebesgue measure. Notice that the birth rate b vanishes in the limit, as in case of CTMC model, because we did not assume any scaling of the birth rate with respect to n.

Let y_B denote the limiting proportion of B molecules in the system. Then, y_B can be described entirely in terms of the density y_A as a solution to the ODE:

$$\begin{equation}\frac{\enspace \mathrm{d}}{\enspace \mathrm{d}t}{y}_{B}\left(t\right)={\int }_{0}^{\infty }\tau \left(s\right){y}_{A}\left(t,s\right)\enspace \mathrm{d}s,\end{equation} \tag{ 2.7 }$$

with the initial condition y_B (0) = 0. Luckily, the limiting system equation (2.6) can be solved explicitly using standard analysis techniques:

$$\begin{equation*}{y}_{A}\left(t,s\right)={f}_{A}\left(s-t\right){S}_{\tau }\left(s\right){S}_{d}\left(s\right)/\left({S}_{\tau }\left(s-t\right){S}_{d}\left(s-t\right)\right),\end{equation*}$$

where S_τ and S_d are the survival functions of the probability distributions characterized by the hazard functions τ and d respectively (see appendix A for the definition of a survival function). Therefore, the limiting concentration of B molecules can be described by

$$\begin{equation*}{y}_{B}\left(t\right)={\int }_{0}^{t}{\int }_{0}^{\infty }\tau \left(v\right){y}_{A}\left(u,v\right)\enspace \mathrm{d}v\enspace \mathrm{d}u.\end{equation*}$$

In figure 1, we numerically show the agreement between the theoretical limits in equations (2.6) and (2.7) and the stochastic simulation. More specifically, we compare ${\int }_{0}^{\infty }{y}_{A}\left(t,s\right)\enspace \mathrm{d}s$ with stochastic simulations of $\langle {n}^{-1}{X}_{t}^{A},1\rangle$ and y_B (t), with $\langle {n}^{-1}{X}_{t}^{B},1\rangle$. As it can be verified, the approximation error vanishes in the limit. Because X_t is a Markov process, the simulation of the stochastic CRN in equation (2.1) can be carried out by adapting the Doob–Gillespie's SSA, which involves simulating two quantities at each step: (1) simulating the next reaction time; and (2) determining the reaction type. Note that, for the CRN in equation (2.1), there are (2N_A (t) + 1) different reactions possible at time t, even though there are only three types of reactions. The next reaction time can be simulated by drawing an exponential random variable with rate equal to the total hazard (the sum of the hazards corresponding to those (2N_A (t) + 1) possible reactions). The total hazard is given by $b+\langle {X}_{t}^{A},\tau \rangle +\langle {X}_{t}^{A},d\rangle$. The type of reaction is then decided by drawing a categorical random variable whose probability masses are the ratios of the individual hazards and the total hazard. This discrete event simulation algorithm is a straightforward adaptation of Doob–Gillespie's SSA for CTMCs. However, it must be noted that the simulation of a non-Markovian CRN is computationally more expensive than the CTMCs. For the sake of completeness, a pseudocode describing the above procedure is given in algorithm 2.1. An implementation in the Julia programming language [46] is also made available in [38].

Zoom In Zoom Out Reset image size

Algorithm 2.1. Pseudocode for the exact simulation of the CRN in equation (2.1).

Require n, X0, K	⊳ K: Maximum number of iterations
Ensure $\left({t}_{1},{X}_{{t}_{1}}\right),\left({t}_{2},{X}_{{t}_{2}}\right),\dots$	⊳ Timings of the reactions along with the measures
1:  Set t = 0
2:  for i = 1, 2, ..., K do	⊳ Compute the next reaction time
3:    Calculate ${\Lambda}={\left(b+\langle {X}_{{t}_{i-1}}^{A},\tau \rangle +\langle {X}_{{t}_{i-1}}^{A},d\rangle \right)}^{-1}$	⊳ Λ−1: Total hazard
4:    if 0 < Λ < ∞ then
5:       Draw an exponential random variable T with mean Λ, i.e., T ∼ Exponential(Λ)	⊳ Advance time to the next reaction time
6:       Set ti = ti−1 + T	⊳ Determine the reaction type
7:       Define π1 = Λb	⊳ Probability for the birth reaction
8:       Define ${\pi }_{j}={\Lambda}\tau \left({\sigma }_{j-1}\left({X}_{{t}_{i-1}}^{A}\right)\right)$ for j = 2, 3, ..., (NA (ti−1) + 1)	⊳ Probabilities for the transformation reaction
9:       Define ${\pi }_{j}={\Lambda}d\left({\sigma }_{j-{N}_{A}\left({t}_{i-1}\right)-1}\left({X}_{{t}_{i-1}}^{A}\right)\right)$ for j = (NA (ti−1) + 2), (NA (ti−1) + 3), ..., (2NA (ti−1) + 1)	⊳ Probabilities for the death reaction
10:       Set $\pi {:=}\left({\pi }_{1},{\pi }_{2},\dots ,{\pi }_{2{N}_{A}\left({t}_{i-1}\right)+1}\right)$
11:       Draw a categorical random variable L with probability distribution π
12:       if L = 1 then	⊳ Birth reaction
13:          ${X}_{{t}_{i}}^{A}={\delta }_{0}+{\sum }_{k=1}^{{N}_{A}\left({t}_{i-1}\right)}{\delta }_{{\sigma }_{k}\left({X}_{{t}_{i-1}}^{A}\right)+T}$	⊳ Advance ages of all A molecules by T and add an atom {0}
14:          ${X}_{{t}_{i}}^{B}={\sum }_{k=1}^{{N}_{B}\left({t}_{i-1}\right)}{\delta }_{{\sigma }_{k}\left({X}_{{t}_{i-1}}^{B}\right)+T}$	⊳ Advance ages of all B molecules by T
15:       else if L ⩽ (NA (ti−1) + 1) then	⊳ Transformation reaction
16:          ${X}_{{t}_{i}}^{A}={\sum }_{k=1}^{{N}_{A}\left({t}_{i-1}\right)}{\delta }_{{\sigma }_{k}\left({X}_{{t}_{i-1}}^{A}\right)+T}-{\delta }_{{\sigma }_{L-1}\left({X}_{{t}_{i-1}}^{A}\right)+T}$	⊳ Remove the atom $\left\{{\sigma }_{L-1}\left({X}_{{t}_{i-1}}^{A}\right)\right\}$ from the measure ${X}_{{t}_{i-1}}^{A}$ and advance ages of all other A molecules by T
17:          ${X}_{{t}_{i}}^{B}={\delta }_{0}+{\sum }_{k=1}^{{N}_{B}\left({t}_{i-1}\right)}{\delta }_{{\sigma }_{k}\left({X}_{{t}_{i-1}}^{B}\right)+T}$	⊳ Advance ages of all B molecules by T and add an atom {0}
18:       else	⊳ Death reaction
19:          ${X}_{{t}_{i}}^{A}={\sum }_{k=1}^{{N}_{A}\left({t}_{i-1}\right)}{\delta }_{{\sigma }_{k}\left({X}_{{t}_{i-1}}^{A}\right)+T}-{\delta }_{{\sigma }_{L-{N}_{A}\left({t}_{i-1}\right)-1}\left({X}_{{t}_{i-1}}^{A}\right)+T}$	⊳ Remove the atom$\left\{{\sigma }_{L-{N}_{A}\left({t}_{i-1}\right)-1}\left({X}_{{t}_{i-1}}^{A}\right)\right\}$ from the measure ${X}_{{t}_{i-1}}^{A}$ and advance ages of all other A molecules by T
20:          ${X}_{{t}_{i}}^{B}={\sum }_{k=1}^{{N}_{B}\left({t}_{i-1}\right)}{\delta }_{{\sigma }_{k}\left({X}_{{t}_{i-1}}^{B}\right)+T}$	⊳ Advance ages of all B molecules by T
21:       end if
22:    else
23:       Stop and break loop
24:    end if
25:    Set i = i + 1.
26:   end for

In section 1, we mentioned that introduction of delay into a CRN could also serve the purpose of model reduction. Indeed, the LLN limit y := (y_A , y_B ) provides a model reduction of the original non-Markovian CRN in equation (2.1). In the following, we discuss two other examples of usefulness of the LLN limit in the form of a PDE system. The first one approximates mean first passage times while the second one describes a faster simulation algorithm.

Mean first passage times are important quantities in the study of stochastic processes and dynamical systems. In the context of CRNs, they could arise in several ways [47, 48]. For instance, natural questions that could arise for the CRN in equation (2.6) are how long it takes to deplete all molecules of species A or to produce the first molecule of B. One of the benefits of the LLN limit is that it can be used to approximate FPTs when the scaling parameter n is sufficiently large. The following illustrates this point.

Suppose we are interested in the time required to produce the first molecule of B. Following the exact simulation algorithm 2.1 adapted from Doob–Gillespie's SSA, the total hazard for the production of a B molecule is $\langle {X}_{0}^{A},\tau \rangle$. In the large-volume limit, we can approximate this hazard by ${\int }_{0}^{\infty }n\tau \left(s\right){y}_{A}\left(0,s\right)\enspace \mathrm{d}s$. Therefore, for a sufficiently large n, the MFPT can be approximated by

$$\begin{equation}m={\left({\int }_{0}^{\infty }n\tau \left(s\right){y}_{A}\left(0,s\right)\enspace \mathrm{d}s\right)}^{-1},\end{equation} \tag{ 2.8 }$$

which, of course, vanishes in the limit of n → ∞. Moreover, the FPTs can be approximated by a random variable following an exponential distribution with mean m, whenever n is sufficiently large. It follows that we can use a simple likelihood function (based on the exponential distribution) for the purpose of statistical inference of the underlying parameters, provided we have observations on the FPTs. This method, called dynamic survival analysis, of estimating parameters based on timings rather than counts was recently explored in the context of epidemiology in [34]. Dynamic survival analysis of general CRNs will be discussed elsewhere.

In figure 2, we show the accuracy of this approximation when n = 100. The approximation appears to be reasonably accurate. More importantly, this suggests we might be able to devise an efficient simulation algorithm using such approximate results. We explore this idea next.

Zoom In Zoom Out Reset image size

Consider a situation when the species A is abundantly available at the beginning of the reaction. Naturally, we expect the PDE approximation to the age density of the species A to be quite accurate, even though there will be considerable stochastic fluctuations in the copy numbers of B, at least initially. However, if we approximate the age density of A by the limiting PDE, we can also approximate the initial growth of the B molecules by a Poisson process whose time-varying intensity is driven by the PDE. We use this idea to devise a hybrid simulation algorithm, which is, again, essentially an adaptation of the Doob–Gillespie's SSA in the sense that it only draws next reaction times from an exponential distribution whose mean depends on the solution to the PDE. For the sake of completeness, a pseudocode describing the idea is provided in algorithm 2.2.

Algorithm 2.2. Pseudocode for the hybrid simulation algorithm.

Require n, yA , K	⊳ K: Maximum number of reactions
Ensure t1, t2, ...	⊳ Timings of creation of B molecules
1:  Set t0 = 0
2:  for i = 1, 2, ..., K do
3:    Calculate ${\Lambda}={\left({\int }_{0}^{\infty }n\tau \left(s\right){y}_{A}\left({t}_{i-1},s\right)\enspace \mathrm{d}s\right)}^{-1}$.
4:    if 0 < Λ < ∞ then
5:       Draw an exponential random variable T with mean Λ, i.e., T ∼ Exponential(Λ)
6:       Set ti = ti−1 + T
7:    else
8:       Stop and break loop
9:    end if
10:    Set i = i + 1.
11:  end for

In figure 3, we show the accuracy of the hybrid simulation algorithm. Expectedly, the hybrid simulation is considerably faster than the full stochastic simulation of the CRN in equation (2.1). A more elaborate comparison of performance is shown in figure 4. However, it is worth noting that the hybrid simulation algorithm, by design, will underestimate the variance in the counting process for the species B. Therefore, one should use the hybrid simulation when it suffices to get the mean trajectory accurately. Alternatively, one can borrow ideas to estimate the variance correctly in other simulation algorithms [49–51]. Similar ideas to expedite simulations have been proposed previously. For instance, Ganguly et al [52] propose a jump-diffusion approximation to the stochastic CRNs and provide error analysis while others [28, 53] introduce hybrid simulation methods using a piecewise deterministic Markov process.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We assume the binding reaction depends on the age of the participating molecule of the enzyme. That is, only k₁ depends on the age of the E molecules (and not on the age of the S molecules); k₋₁ and k₂ are constants. For the species E, S, C, and P, define the measure-valued stochastic processes

$$\begin{align*}\hfill & {X}_{t}^{E}{:=}\sum _{i=1}^{{N}_{E}\left(t\right)}{\delta }_{{e}_{i}\left(t\right)},\quad {X}_{t}^{S}{:=}\sum _{i=1}^{{N}_{S}\left(t\right)}{\delta }_{{s}_{i}\left(t\right)},\quad \hfill \\ \hfill & {X}_{t}^{C}{:=}\sum _{i=1}^{{N}_{C}\left(t\right)}{\delta }_{{c}_{i}\left(t\right)},\quad {X}_{t}^{P}{:=}\sum _{i=1}^{{N}_{P}\left(t\right)}{\delta }_{{p}_{i}\left(t\right)},\hfill \end{align*}$$

where N_E , N_S , N_C , N_P denote the counts of molecules of E, S, C, and P respectively. Similarly, e_i , s_i , c_i , p_i denote the age of the ith molecule of E, S, C, and P respectively. The process X := (X^E , X^S , X^C , X^P ) is a Markov process on the space $D\left(\left[0,T\right],{\mathcal{M}}_{P}{\left({\mathbb{R}}_{+}\right)}^{4}\right)$. Please note that we need to scale the hazard function k₁ corresponding to the bimolecular reaction by n⁻¹ following the stochastic law of mass actions [1].

As before, we are interested in the large-volume limit of the scaled process n⁻¹ X_t . The scaled stochastic process n⁻¹ X_t converges to a deterministic function ${x}_{t}{:=}\left({x}_{t}^{E},{x}_{t}^{S},{x}_{t}^{C},{x}_{t}^{P}\right)$ whose components ${x}_{t}^{E},{x}_{t}^{S},{x}_{t}^{C},{x}_{t}^{P}$ are finite measures on ${\mathbb{R}}_{+}$ by virtue of the LLN.

Let y_E denote the density of the measure ${x}_{t}^{E}$ with respect to the Lebesgue measure. Also, let y_S , y_C , y_P denote the concentrations of the S, C, and P molecules. Then, we get the following limiting system:

$$\begin{equation}\begin{aligned}\hfill \left({\partial }_{t}+{\partial }_{s}\right){y}_{E}\left(t,s\right)& =-{k}_{1}\left(s\right){y}_{E}\left(t,s\right){y}_{S}\left(t\right),\hfill \\ \hfill \frac{\enspace \mathrm{d}}{\enspace \mathrm{d}t}{y}_{S}\left(t\right)& =-{y}_{S}\left(t\right){\int }_{0}^{\infty }{k}_{1}\left(s\right){y}_{E}\left(t,s\right)\enspace \mathrm{d}s\hfill \\ \hfill & \quad +{k}_{-1}{y}_{C}\left(t\right),\hfill \\ \hfill \frac{\enspace \mathrm{d}}{\enspace \mathrm{d}t}{y}_{C}\left(t\right)& ={y}_{S}\left(t\right){\int }_{0}^{\infty }{k}_{1}\left(s\right){y}_{E}\left(t,s\right)\enspace \mathrm{d}s\hfill \\ \hfill & \quad -\left({k}_{-1}+{k}_{2}\right){y}_{C}\left(t\right),\hfill \\ \hfill \frac{\enspace \mathrm{d}}{\enspace \mathrm{d}t}{y}_{P}\left(t\right)& ={k}_{2}{y}_{C}\left(t\right),\hfill \end{aligned}\end{equation} \tag{ 3.11 }$$

with the boundary condition

$$\begin{equation*}{y}_{E}\left(t,0\right)=\left({k}_{-1}+{k}_{2}\right){y}_{C}\left(t\right)\end{equation*}$$

and the initial condition y_E (0, s) = f_E (s) such that ${\int }_{0}^{\infty }{f}_{E}\left(s\right)\enspace \mathrm{d}s=\left[{E}_{0}\right]$. Appropriate initial conditions for S, C, and P are also assumed. This limiting system can now be used to study the effects of delay in the binding reaction. One interesting approximation that has been widely applied in the context of Michaelis–Menten enzyme kinetic reactions is what is known as a quasi-steady state approximation [58]. There are many forms of QSSAs, namely, standard QSSA (sQSSA), total QSSA (tQSSA), and reversible QSSA (rQSSA). Detailed analysis of any of the QSSAs is beyond the scope of the present work. For the purpose of illustration, we informally describe an analogue of the sQSSA here.

The QSSAs are a multiscale approximation of the Michaelis–Menten enzyme-kinetic reactions. The basic assumption behind the standard QSSA is that the substrate-enzyme complex C reaches its steady-state much quicker than the other species. In the deterministic set-up, the approximation is achieved by setting $\frac{\enspace \mathrm{d}}{\enspace \mathrm{d}t}{y}_{C}\left(t\right)=0$ in equation (3.11), which allows one to work with a smaller system of ODEs. Several conditions for the validity of the sQSSA have been proposed in the literature. See [57] for a detailed discussion.

Following the deterministic approach in our case, we set $\frac{\enspace \mathrm{d}}{\enspace \mathrm{d}t}{y}_{C}\left(t\right)=0$ in equation (3.11) to get a reduced PDE system that is analogous to the sQSSA. To be more precise, $\frac{\enspace \mathrm{d}}{\enspace \mathrm{d}t}{y}_{C}\left(t\right)=0$ yields

$$\begin{equation*}{y}_{C}\left(t\right)=\frac{{y}_{S}\left(t\right){\int }_{0}^{\infty }{k}_{1}\left(s\right){y}_{E}\left(t,s\right)\enspace \mathrm{d}s}{{k}_{-1}+{k}_{2}},\end{equation*}$$

which further yields an approximate system

$$\begin{equation}\begin{aligned}\hfill \frac{\enspace \mathrm{d}}{\enspace \mathrm{d}t}{y}_{S}\left(t\right)& =-\frac{{k}_{2}}{{k}_{-1}+{k}_{2}}{y}_{S}\left(t\right){\int }_{0}^{\infty }{k}_{1}\left(s\right){y}_{E}\left(t,s\right)\enspace \mathrm{d}s,\hfill \\ \hfill \frac{\enspace \mathrm{d}}{\enspace \mathrm{d}t}{y}_{P}\left(t\right)& =\frac{{k}_{2}}{{k}_{-1}+{k}_{2}}{y}_{S}\left(t\right){\int }_{0}^{\infty }{k}_{1}\left(s\right){y}_{E}\left(t,s\right)\enspace \mathrm{d}s.\hfill \end{aligned}\end{equation} \tag{ 3.12 }$$

Recall that y_E solves $\left({\partial }_{t}+{\partial }_{s}\right){y}_{E}\left(t,s\right)=-{k}_{1}\left(s\right){y}_{E}\left(t,s\right){y}_{S}\left(t\right)$ with boundary condition y_E (t, 0) = (k₋₁ + k₂)y_C (t) and initial condition y_E (0, s) = f_E (s). As a consequence, y_E is determined by y_S and y_C , and can be partially solved in terms of y_S and y_C . Therefore, the reduced system of ODEs in equation (3.12) is indeed autonomous and therefore, serves as an sQSSA of the CRN in equation (3.9).

In the stochastic set-up, the QSSAs are obtained by means of the probabilistic multiscaling techniques developed in [3, 4]. The stochastic and the deterministic QSSAs mostly agree with each other with some notable differences. Please see [57] for examples of discrepancies as well as more details on the methods. Here, for paucity of space, we do not consider the stochastic QSSAs or possible discrepancies between stochastic and deterministic methods in the present age-structured models.

First, under the above mentioned assumptions, we can show that the components of the scaled process n⁻¹ X_t do not explode (similar derivation in [37, lemma 2.6 and proposition 2.7]). Now, note that the trajectory equations for the processes ${X}_{t}^{A}$ and ${X}_{t}^{B}$ given in equations (2.4) and (2.5) are driven by PPMs. Since we have

$$\begin{align*}\hfill {f}_{t}\left(a+t-s\right)& ={f}_{s}\left(a\right)+{\int }_{s}^{t}\left(\frac{\partial }{\partial u}{f}_{u}\left(a+u-s\right)\right.\hfill \\ \hfill & \quad \left.+\frac{\partial }{\partial a}{f}_{u}\left(a+u-s\right)\right)\enspace \mathrm{d}u,\hfill \end{align*}$$

and using the compensated PPMs of the PPMs Q₁, Q₂, Q₃, we can show the processes

$$\begin{equation*}\begin{aligned}\hfill {M}_{t}^{A,f}& =\langle {n}^{-1}{X}_{t}^{A},{f}_{t}\rangle -\langle {n}^{-1}{X}_{0}^{A},{f}_{0}\rangle \hfill \\ \hfill & \quad -{\int }_{0}^{t}{\int }_{0}^{\infty }\left(\frac{\partial }{\partial a}{f}_{s}\left(a\right)+\frac{\partial }{\partial s}{f}_{s}\left(a\right)\right.\hfill \\ \hfill & \quad \left.-{f}_{s}\left(a\right)\left(\tau \left(a\right)+d\left(a\right)\right)\right){n}^{-1}{X}_{s}^{A}\left(\enspace \mathrm{d}a\right)\enspace \mathrm{d}s\hfill \end{aligned}\end{equation*}$$

$$\begin{equation*}\begin{aligned}\hfill {M}_{t}^{B,f}& =\langle {n}^{-1}{X}_{t}^{B},{f}_{t}\rangle -{\int }_{0}^{t}{\int }_{0}^{\infty }\left(\frac{\partial }{\partial a}{f}_{s}\left(a\right)+\frac{\partial }{\partial s}{f}_{s}\left(a\right)\right.\hfill \\ \hfill & \quad \left.+{f}_{s}\left(0\right)\tau \left(a\right)\right){n}^{-1}{X}_{s}^{A}\left(\enspace \mathrm{d}a\right)\enspace \mathrm{d}s\hfill \end{aligned}\end{equation*}$$

are zero mean, square integrable, càdlàg martingale processes with predictable quadratic variations of the order n⁻¹. Since we expect the predictable quadratic variations to vanish in the limit of n → ∞, the scaled process n⁻¹ X_t converges to a deterministic, continuous function x_t . The tightness of the process n⁻¹ X_t can be established by verifying a criterion due to Roelly [80] in the vague topology and the Aldous–Rebolledo criteria [81]. See [36] or [37, proposition 3.1] for similar calculations. Furthermore, thanks to the martingale representations above, we expect the limit x_t to satisfy

$$\begin{align*}\hfill \langle {x}_{t}^{A},{f}_{t}\rangle & =\langle {x}_{0}^{A},{f}_{0}\rangle +{\int }_{0}^{t}{\int }_{0}^{\infty }\left(\frac{\partial }{\partial a}{f}_{s}\left(a\right)+\frac{\partial }{\partial s}{f}_{s}\left(a\right)\right.\hfill \\ \hfill & \quad \left.-{f}_{s}\left(a\right)\left(\tau \left(a\right)+d\left(a\right)\right)\right){x}_{s}^{A}\left(\enspace \mathrm{d}a\right)\enspace \mathrm{d}s\hfill \\ \hfill \langle {x}_{t}^{B},{f}_{t}\rangle & ={\int }_{0}^{t}{\int }_{0}^{\infty }\left(\frac{\partial }{\partial a}{f}_{s}\left(a\right)+\frac{\partial }{\partial s}{f}_{s}\left(a\right)\right.\hfill \\ \hfill & \quad \left.+{f}_{s}\left(0\right)\tau \left(a\right)\right){x}_{s}^{A}\left(\enspace \mathrm{d}a\right)\enspace \mathrm{d}s.\hfill \end{align*}$$

The uniqueness of the solutions can be shown by first establishing that the solutions remain bounded on finite time intervals (recall the global jump rates are assumed bounded) and then invoking Grönwall's lemma to show the distance between two possible solutions must vanish proving the desired uniqueness.