In recent years, single-cell and multi-omics data analysis has attracted attention in the field of molecular cell biology (Lim et al., 2023; Massimino et al., 2023; Baysoy et al., 2024). In particular, analysis using machine learning models, such as deep learning, has been active (Lin et al., 2022; Gunawan, et al., 2023; Wagle et al., 2024). However, understanding stoichiometric reactions is still considered essential to understand gene expression and protein enzyme functions in detail. This paper proposes a simulation program for biochemical reaction systems based on the natural number simulation (NNS) method. NNS offers a simplified approach for modeling biochemical reactions stoichiometrically and performing computations of the time course of the reactions. This method excels in its straightforward formulation process and is particularly suitable for complex biological systems (Blinov et al., 2004; Browning et al., 2022; Huizing et al., 2022) because of the absence of complicated reaction formulas using mathematical equations.

Various models have been considered to represent complex biological systems, including ordinary differential equations and Gillespie’s algorithm (Gillespie, 1977; Ullah and Wolkenhauer, 2011; Alon, 2019). For example, metabolic reaction systems continue to be vigorously analyzed in detail using differential equations (Khodayari and Maranas, 2016; Himeoka and Mitarai, 2022). Petri nets substantially contribute to analyzing biochemical reaction networks by integrating stochastic algorithms and/or non-parametric strategies. This progression has resulted in the emergence of specialized forms such as signaling Petri nets, large-scale metabolic models, and multilevel biological models using colored Petri nets, indicating their increasing utility in the analysis of complex systems (Ruths et al., 2008; Rohr, 2018; Brinkrolf et al., 2021; Shaikh et al., 2022; Liu et al., 2023).

However, ordinary differential equations require complex formulations of reaction rate equations. Stochastic simulations require detailed formulations of chemical reaction networks, which may limit their usefulness, particularly for large-scale simulations. Petri net-based simulations offer a comprehensive framework for formulating chemical reaction networks involving specific elements such as places, transitions, arcs, and markings. However, modeling actual biochemical reaction systems requires consideration of the definition of the number of tokens and their correspondence to the number of molecules. In addition, when introducing stochastic models, it is necessary to explicitly define the timing of firing.

The NNS method implemented in Python is very easy to model because the procedure for determining detailed time evolution is included in the computational algorithm. Requiring only the stoichiometric equation, its rate constants, and the initial number of molecules, the new algorithm, based on a probabilistic binomial distribution, immediately calculates the number of molecules after one calculation step for the elements in the model system. It is crucial to emphasize the necessity of precise rate constant determination to depict the system accurately. Nevertheless, the simplicity of the formulation has the advantage of allowing the optimization process to be easily performed.

Recent advancements in non-parametric analytical techniques also present an intriguing possibility for NNS (Ruths et al., 2008; Rohr, 2018). By meticulously setting the initial molecular counts, rate constants, and stoichiometric equations, an analysis comparable to those seen in studies of complex signaling networks may be achieved. This approach could potentially enhance the applicability and accuracy of NNS for modeling intricate biochemical systems.

Contrary to traditional methods, such as ordinary differential equations and Gillespie’s algorithm, which might not optimally represent complex biological systems (Gillespie, 1977; Alon, 2019), NNS provides a more precise simulation even for a small number of molecules in DNA-related reactions. It uses stoichiometric equations with rate constants for specific and selective processes such as transcription and translation (Watson et al., 2014). Moreover, the binomial distribution in NNS facilitates calculating informational and entropic metrics, including Shannon’s entropy (Shannon, 1948; Adriaans and van Benthem, 2008). Examples of these calculations are discussed later in the paper.

Our method calculates time-developing stoichiometric reactions using a binomial algorithm. Consider the general reaction as Equation 1: q 1 X 1 + q 2 X 2 + ⋯ + q s X s → k r 1 Y 1 + r 2 Y 2 + ⋯ + r t Y t (1)

X ₁ , X ₂ , ⋯ , X _s denote molecular elements before the reaction while Y ₁ , Y ₂ , ⋯ , Y _t represent those after the reaction. The coefficients q_i and r_i are reaction orders. The subscripts “s” and “t” are natural numbers, with no limitation on their sizes. The rate constant k implies the probability of forming Y ₁ , Y ₂ , ⋯ , Y _t . In the NNS approach, reaction orders are consistently natural numbers because they necessitate dealing with element counts as natural numbers at all times. Our methodology accommodates stoichiometric reactions without limiting specific chemical formulas, thereby ensuring the conservation of the atom number.

Initial element counts are necessary for computing the time evolution. X _i _n and Y _i _n are defined as the elemental number of X _i and Y _i ; then, the following calculation in Equations 2, 3 with Equations 4, 5 provide the decrement and increment quantities of the elements. ∆ X i _ n = − q i × R b i n o m i a l n , p (2) ∆ Y i _ n = r i × R b i n o m i a l n , p (3) n = int X i _ n q i min (4) p = ∏ i ≠ ⁡ min k ∙ X i _ n q i ∙ N 1 + k ∙ X i _ n q i ∙ N (5)

R _binomial (n, p) generates random numbers following the binomial distribution, returning natural numbers (Rohr, 2018). In this function, n denotes a trial number, each of which can be either a success or a failure (binary outcomes). The parameter p specifies the probability of success in each trial. The function provides the total number of successes in n trials. NumPy’s random binomial (n, p) achieves this. The function returns one successful natural number under n, including zero. Stochastic simulation often uses these discrete numbers (Székely and Burrage, 2014; Gholami and Ilie, 2021). The trial number n is represented with the minimum integer value among X ₁ _n/q ₁ , X ₂ _n/q ₂ , X ₃ _n/q ₃ , so on. The probability p is derived based on the rate constant k, the reaction order q _i , a normalization parameter N, and X _i _n, except for k ∙ X i _ n q i ∙ N 1 + k ∙ X i _ n q i ∙ N with X i _ n q i selected in Equation 4. In a calculation step, X _i _n becomes X _i _n + ΔX _i _n and Y _i _n becomes Y _i _n + ΔY _i _n. The rate constant k in NNS is semantically the same as the rate constant that appears in the stoichiometric equation defined by concentration. However, k in this study is defined as a constant that affects the stochastic results defined in Equation 5. The rate constant k, the number of elements, and the normalization constant N determine the probability p, as indicated by Equation 5. The k value influences the probability variation between 0 and 1 and the reaction rate, but it is not the same as the conventional rate constant.

Equations 4, 5 have some minor formulas. One is in the case of a one-element reaction as Equation 6: aX → r 1 Y 1 + r 2 Y 2 + ⋯ + r t Y t  with  a  rate constant  k (6)

X is a molecular element, and a is an order. The decrement ΔX_n and increment ΔY _i _n are derived from the following formulas as Equations 7, 8 with Equations 9, 10: ∆ X _ n = − a × R b i n o m i a l n , p (7) ∆ Y i _ n = r i × R b i n o m i a l n , p (8) n = int X _ n a (9) p = k 1 + k (10) where X_n is the number of the element X.

Another exceptional reaction type is as Equation 11: aX  → 0   zero , with  a  rate constant  k (11)

This expression indicates that element X is decomposed and disappears, so Equation 3 does not work. Thus, only X decreases with Equation 7 with Equation 9, 10.

Another type of reaction is defined as a linear decreasing and increasing reaction for one element. One needs these reactions for the mathematical formulation of linear changes via molecular addition and subtraction. The linear decreasing case is defined similarly to Equation 11 as Equation 12: aX  → 0   zero  by linear , with  a  rate constant  k (12)

Moreover, the increment of X is determined by the following Equation 13: ∆ X _ n = − a × R b i n o m i a l 1 , p (13) where p is the same as in Equation 10.

Subsequently, the linear increasing case is defined by the following reaction as Equation 14: 0   z e r o →   b Y  by linear , with  a  rate constant  k (14)

The following Equation 15 also determines the increment of Y: ∆ Y _ n = b × R b i n o m i a l 1 , p (15) where p is the same as in Equation 10.

A biological system has numerous incredible reactions (Nelson and Cox, 2021). One can easily define each reaction if one knows the elements before and after the reaction, their orders, and rate constants.

Executing an input file in Spyder, an Integrated Development Environment requires setting “fName = “inp_file.txt”” in the main program (i.e., binomial_v016.py) and then utilizing the “Run” command. For command line execution, the command “$ python binomial_v016.py inp_file.txt” yields similar outcomes to Spyder’s run. Therefore, it is essential to position the input file in the program’s directory or one level above it.

This file needs to define calculation time, elements, reactions, and plots (Table 1). Currently, model systems are defined in text files. Currently, a Python program is available on GitHub that converts xml files that conform to Systems Biology Markup Language into text files for NNS. Users can perform calculations via NNS by setting the appropriate initial number of elements and rate constants.

The *Reaction statement has numerous definitions (Supplementary Tables S1, S2), including the above formulas used to calculate reaction increments.

The *ElementInOut statement implemented definition statements to express the inflow and outflow of molecules in and out of the system (Supplementary Table S3). The author developed the Python program used for this algorithm, and it is available under the MIT license on GitHub at https://github.com/taka-b/binomial/tree/binomial_v016_01 (Sato, 2023). This algorithm is publicly accessible and can be used by anyone without any charge. Figure 1A shows the calculation for schematic flow, while Figure 1B depicts the detailed calculation algorithm. The simplicity of the proposed procedure is demonstrated through six examples using their corresponding input files.

Some example calculations are shown below. Note that the user only needs to prepare a text file for the calculations to obtain the results. The important point is to carefully formulate the stoichiometric equation and give appropriate rate constants.

The following results are obtained for calculations based on the text file shown in each figure. The start and end times are defined by *Time; for example, zero, 10,000, displays the results of calculations for every step from zero to 10,000. If you want to compare these results with actual experimental results, you will need units and rate constants that correspond to the experimental results. Here, however, the units of calculation are denoted in Steps, and figures show changes in the number of elements.

Entropy and information are concepts related to information theory (Baez and Pollard, 2016), and stochastic processes with probability parameters can enrich our understanding of reaction systems. In our approach, parameters n and p are determined via reaction conditions and are described in the Methods section. The information on reaction is defined as Equation 20: I _ r e a c t i o n t , X t ; n t , p t = − ⁡ log 2 b i n o m i a l X t ; n t , p t (20) where binomial ( X ; n , p ) represents a binomial distribution function with a stochastic variable X (Successes) determined using n and p for each reaction step t (Chanda et al., 2020), as illustrated in Supplementary Figure S9.

The I_reaction provides insight into the reaction dynamics, and its value increases as reactions proceed.

Drawing from the Shannon entropy (Ben-Naim, 2012), the reaction entropy (RE) can be defined as Equation 21: R E t ; n t , p t = − ∑ X t = 0 n t b i n o m i a l X t ; n t , p t log 2 b i n o m i a l X t ; n t , p t (21)

The RE is illustrated in Supplementary Figure S10 as a function of n and p, and it captures the inherent randomness in reaction and activity of reactions (Baez and Pollard, 2016; Roach, 2020; Uda, 2020).

Using a DNA-type reaction example from Figure 7 and Supplementary Figure S4, the accumulated information on the reaction is calculated using Equation 22. I R a c c u m u l a t e t , X t ; n t , p t = − ∑ t k = 0 t log 2 b i n o m i a l X t k ; n t k , p t k (22)

Figure 10 shows the results of Equation 22, and Figure 11 presents the RE values. Here, the RE values represent the results obtained from Equation 21 using Figure 7 and Supplementary Figure S4. Information on reaction and RE offers fundamental tools for dissecting reaction properties and understanding the complexity of biochemical systems.