^{1,a),b)}, André Leier

^{2,a),b)}and Tatiana T. Marquez-Lago

^{3,b)}

### Abstract

Accurate modelling and simulation of dynamic cellular events require two main ingredients: an adequate description of key chemical reactions and simulation of such chemical events in reasonable time spans. Quite logically, posing the right model is a crucial step for any endeavour in Computational Biology. However, more often than not, it is the associated computational costs which actually limit our capabilities of representing complex cellular behaviour. In this paper, we propose a methodology aimed at representing chains of chemical reactions by much simpler, reduced models. The abridgement is achieved by generation of model-specific delay distribution functions, consecutively fed to a delay stochastic simulation algorithm. We show how such delay distributions can be analytically described whenever the system is solely composed of consecutive first-order reactions, with or without additional “backward” bypass reactions, yielding an exact reduction. For models including other types of monomolecular reactions (constitutive synthesis, degradation, or “forward” bypass reactions), we discuss why one must adopt a numerical approach for its accurate stochastic representation, and propose two alternatives for this. In these cases, the accuracy depends on the respective numerical sample size. Our model reduction methodology yields significantly lower computational costs while retaining accuracy. Quite naturally, computational costs increase alongside network size and separation of time scales. Thus, we expect our model reduction methodologies to significantly decrease computational costs in these instances. We anticipate the use of delays in model reduction will greatly alleviate some of the current restrictions in simulating large sets of chemical reactions, largely applicable in pharmaceutical and biological research.

M.B. would like to thank Carlos Marijuán for helpful discussions on eigenvalues.

I. INTRODUCTION

II. RESULTS

A. Consecutive reversible reactions

B. Large chemical systems

C. Fast reversible reactions with slow turnover

D. Additional forward and backward bypass reactions

E. Systems with additional degradation reactions

F. Systems with additional incoming reactions

G. Lumping binary reactions

H. Application: Chains of reactions in mRNA decay

III. METHODS

A. Determinant of A

B. Eigenvalues of A

C. Resolvent of A

D. First-passage time PDF

E. Backward bypass reactions

F. Numerical solution of

G. Degradation

H. Simulations

IV. DISCUSSION

### Key Topics

- Eigenvalues
- 42.0
- Chemical reactions
- 25.0
- Biochemical reactions
- 19.0
- Addition reactions
- 15.0
- Reaction rate constants
- 14.0

## Figures

Linear sequence of unimolecular reactions (example 1). (Left) Histogram for the number of *S* _{9} molecules at time T = 80 obtained from 10^{4} SSA (blue) and DSSA (red) simulations where delays were drawn from the 8-exponential iCDF with parameters (the eigenvalues of the system's rate matrix *A* for species *S* _{1}–*S* _{8}). (Right) Average time evolution of species *S* _{9} in steps of Δt = 2 until T = 80 for SSA and DSSA.

Linear sequence of unimolecular reactions (example 1). (Left) Histogram for the number of *S* _{9} molecules at time T = 80 obtained from 10^{4} SSA (blue) and DSSA (red) simulations where delays were drawn from the 8-exponential iCDF with parameters (the eigenvalues of the system's rate matrix *A* for species *S* _{1}–*S* _{8}). (Right) Average time evolution of species *S* _{9} in steps of Δt = 2 until T = 80 for SSA and DSSA.

Linear sequence of reactions with time-scale separation (example 2). (Left) Histogram for the number of *S* _{9} molecules at time T = 1000 for 10^{6} SSA simulations (blue) and 10^{4} DSSA simulations where delays were drawn from the corresponding 8-exponential iCDF (red) or 1-exponential iCDF for the smallest absolute eigenvalue (green). (Right) Average time evolution of species *S* _{9} in steps of Δt = 10 for SSA (blue) and DSSA (red/green). The blue and red trajectories are undistinguishable while the green trajectory is above the other two.

Linear sequence of reactions with time-scale separation (example 2). (Left) Histogram for the number of *S* _{9} molecules at time T = 1000 for 10^{6} SSA simulations (blue) and 10^{4} DSSA simulations where delays were drawn from the corresponding 8-exponential iCDF (red) or 1-exponential iCDF for the smallest absolute eigenvalue (green). (Right) Average time evolution of species *S* _{9} in steps of Δt = 10 for SSA (blue) and DSSA (red/green). The blue and red trajectories are undistinguishable while the green trajectory is above the other two.

Sequence of unimolecular reactions with backward bypass reactions (example 3). (Left) Histogram for the number of *S* _{9} molecules at time T = 120 obtained from SSA simulations (blue) and DSSA (red) where delays were drawn from the corresponding 8-exponential iCDF. (Right) Average time evolution of species *S* _{9} in steps of Δt = 10 using SSA and DSSA.

Sequence of unimolecular reactions with backward bypass reactions (example 3). (Left) Histogram for the number of *S* _{9} molecules at time T = 120 obtained from SSA simulations (blue) and DSSA (red) where delays were drawn from the corresponding 8-exponential iCDF. (Right) Average time evolution of species *S* _{9} in steps of Δt = 10 using SSA and DSSA.

Sequence of unimolecular reactions with forward bypass reactions (example 4). (Left) Histogram for the number of *S* _{9} molecules at time T = 30 for 10^{6} SSA simulations (blue) and 10^{4} DSSA simulations (red, green) where delays were obtained from first-passage time distributions based on 10^{2} (red) and 10^{3} (green) sample times. (Right) Average time evolution of species *S* _{9} in steps of Δt = 10 for SSA and DSSA.

Sequence of unimolecular reactions with forward bypass reactions (example 4). (Left) Histogram for the number of *S* _{9} molecules at time T = 30 for 10^{6} SSA simulations (blue) and 10^{4} DSSA simulations (red, green) where delays were obtained from first-passage time distributions based on 10^{2} (red) and 10^{3} (green) sample times. (Right) Average time evolution of species *S* _{9} in steps of Δt = 10 for SSA and DSSA.

Sequence of unimolecular reactions with backward and forward bypass reactions (example 5). (Left) Histogram for the number of *S* _{9} molecules at time T = 50 for 10^{6} SSA simulations (blue) and 10^{4} DSSA simulations where delays were either drawn from first-passage time distributions based on 10^{2} (black x) or 10^{3} (green x) samples or obtained (via inverse sampling) from the numerical solution of the CDF evaluated at 6401 time points in the interval [0, 640] (black dot). (Right) Average time evolution of species *S* _{9} in steps of Δt = 10 for SSA and DSSA.

Sequence of unimolecular reactions with backward and forward bypass reactions (example 5). (Left) Histogram for the number of *S* _{9} molecules at time T = 50 for 10^{6} SSA simulations (blue) and 10^{4} DSSA simulations where delays were either drawn from first-passage time distributions based on 10^{2} (black x) or 10^{3} (green x) samples or obtained (via inverse sampling) from the numerical solution of the CDF evaluated at 6401 time points in the interval [0, 640] (black dot). (Right) Average time evolution of species *S* _{9} in steps of Δt = 10 for SSA and DSSA.

Sequence of unimolecular reactions with additional degradation reactions (example 6). (Left) Histogram for the number of *S* _{9} molecules at time T = 80 for 10^{4} SSA (blue) and DSSA (black dots) simulations where delays were drawn from the inverse of the numerical solution of the CDF (using matrix exponentials). (Right) Average time evolution of species *S* _{9} in steps of Δt = 10 for SSA and DSSA. For this scenario, we computed *w* = 0.194.

Sequence of unimolecular reactions with additional degradation reactions (example 6). (Left) Histogram for the number of *S* _{9} molecules at time T = 80 for 10^{4} SSA (blue) and DSSA (black dots) simulations where delays were drawn from the inverse of the numerical solution of the CDF (using matrix exponentials). (Right) Average time evolution of species *S* _{9} in steps of Δt = 10 for SSA and DSSA. For this scenario, we computed *w* = 0.194.

Histogram (normalized) of species *S* _{9} at time T = 80 from 10 000 simulations of the full (SSA, blue) and the abridged system (DSSA, red).

Histogram (normalized) of species *S* _{9} at time T = 80 from 10 000 simulations of the full (SSA, blue) and the abridged system (DSSA, red).

Histograms of species *S* and *P* at time T = 10 from 10 000 simulations of the full (SSA, blue) and the abridged system (DSSA, red) for the parameter set *c* _{1} = 1, *c* _{2} = 10, *c* _{3} = 10 and initial conditions *S* _{0} = 10^{5}, *E* _{0} = 10^{2}, *ES* _{0} = 0, *P* = 0. Here, we use *k* = 10^{15}.

Histograms of species *S* and *P* at time T = 10 from 10 000 simulations of the full (SSA, blue) and the abridged system (DSSA, red) for the parameter set *c* _{1} = 1, *c* _{2} = 10, *c* _{3} = 10 and initial conditions *S* _{0} = 10^{5}, *E* _{0} = 10^{2}, *ES* _{0} = 0, *P* = 0. Here, we use *k* = 10^{15}.

(a) and (c) Complete model and kinetic parameters used in modelling MFA2pG-mRNA degradation as stated in Ref. ^{ 12 } . 3^{′} fragments (L + I1 + I2) and 5^{′} fragments (G + M) are highlighted in grey, all other species are considered full length mRNA. (b) Abridged model. The dotted line in (a) refers to the process that is represented by a delay distribution in our abridged model. Moreover, reactions C → E and D → F are lumped into a single degradation reaction. The probability for such a degradation is 1 − *w* (cf. Sec. III G ).

(a) and (c) Complete model and kinetic parameters used in modelling MFA2pG-mRNA degradation as stated in Ref. ^{ 12 } . 3^{′} fragments (L + I1 + I2) and 5^{′} fragments (G + M) are highlighted in grey, all other species are considered full length mRNA. (b) Abridged model. The dotted line in (a) refers to the process that is represented by a delay distribution in our abridged model. Moreover, reactions C → E and D → F are lumped into a single degradation reaction. The probability for such a degradation is 1 − *w* (cf. Sec. III G ).

(Left) Histogram for the number of I2 molecules at time T = 10 000 for SSA simulations (blue) and DSSA (black) where delays were drawn from the inverse of the numerical solution of the CDF (using matrix exponentials). (Right) Average time evolution of species I2 for SSA (blue line) and DSSA (black dots). For this scenario, we computed *w* = 0.7745, i.e., a 77.45% probability that the former full length mRNA is degraded via this pathway.

(Left) Histogram for the number of I2 molecules at time T = 10 000 for SSA simulations (blue) and DSSA (black) where delays were drawn from the inverse of the numerical solution of the CDF (using matrix exponentials). (Right) Average time evolution of species I2 for SSA (blue line) and DSSA (black dots). For this scenario, we computed *w* = 0.7745, i.e., a 77.45% probability that the former full length mRNA is degraded via this pathway.

## Tables

Computational savings in terms of average numbers of SSA reactions per single DSSA reaction and speed-up (runtime of SSA over runtime of DSSA runs). Simulations ran until T = 100. Mean values are calculated over 100 simulations.

Computational savings in terms of average numbers of SSA reactions per single DSSA reaction and speed-up (runtime of SSA over runtime of DSSA runs). Simulations ran until T = 100. Mean values are calculated over 100 simulations.

Article metrics loading...

Full text loading...

Commenting has been disabled for this content