## 2 Materials and Methods  
### Reaction-Based Modeling and Stochastic Simulation Algorithm  
Considering the stochastic formulation of chemical kinetics [16], we can define a reaction-based model (RBM) by specifying two disjoint sets:* The set \(\mathcal{S}=\{S_{1},\ldots,S_{N}\}\) of molecular species;
* The set \(\mathcal{R}=\{R_{1},\ldots,R_{M}\}\) of the biochemical reactions describing all interactions among the species in \(\mathcal{S}\). In RBMs, a biochemical reaction is defined as follows: \[R_{m}:\,\sum_{n=1}^{N}\alpha_{m,n}S_{n}\xrightarrow{c_{m}}\,\sum_{n=1}^{N}\, \beta_{m,n}S_{n},\] (1) where \(\alpha_{m,n},\beta_{m,n}\in\mathbb{N}\) indicate the stoichiometric coefficients associated with the \(n\)-th reactant and with the \(n\)-th product of the \(m\)-th reaction, respectively. The value \(c_{m}\in\mathbb{R}^{+}\), associated with reaction \(R_{m}\), is named the stochastic (kinetic) constant and encompasses all its chemical-physical properties. Although deterministic simulations of mechanistic RBMs allow for reproducing the dynamic behavior of biochemical systems [17], only stochastic simulations of RBMs can reproduce the emerging effects due to biological noise, which plays an important role when the molecular species occur in low amounts within the system [18]. The Stochastic Simulation Algorithm (SSA) [16] is a seminal procedure that allows for generating an exact realization of the temporal evolution of an RBM SSA requiring that (i ) all reactions involving the chemical species must take place in a single volume, characterized by constant physical conditions (e.g., pressure and temperature) throughout the simulation; (ii) the molecules must be uniformly distributed within the volume, which is then considered well-stirred; and (iii) the (discrete) quantity of each molecular species is denoted as an integer number \(x_{n}\in\mathbb{N}\), for each \(n=1,\ldots,N\). The state of the system at time \(t\) is denoted in SSA by the vector \(\mathbf{x}=\mathbf{x}(t)\equiv(x_{1}(t),\ldots,x_{N}(t))\), and its content is updated according to the execution of reactions. To be more precise, at each step of the simulation, SSA identifies the reaction that will take place in the time interval \([t,t+\tau)\) according to the probability of each reaction to occur in the infinitesimal time step \([t,t+dt)\), which is proportional to the so-called propensity function of each reaction. The propensity functions are calculated as \(a_{m}(\mathbf{x})=c_{m}d_{m}(\mathbf{x})\), where \(d_{m}(\mathbf{x})\) correspond to the number of distinct combinations of the reactant molecules involved in reaction \(R_{m}\), in the current state of the system \(\mathbf{x}\). Successively, SSA calculates the waiting time \(\tau\) before the reaction is executed, as follows:  
\[\tau=\frac{1}{a_{0}(\mathbf{x})}\ln\!\left(\frac{1}{\xi_{1}}\right),\]  
where \(\xi_{1}\) is a random number sampled from the uniform distribution in \((0,1)\), and \(a_{0}(\mathbf{x})=\sum_{m=1}^{M}a_{m}(\mathbf{x})\). \(m\) denotes the index of the reaction to fire and corresponds to the smallest integer in the interval \([1,M]\) such that  
\[\sum_{m^{\prime}=1}^{m}a_{m^{\prime}}(\mathbf{x})>\zeta_{2}a_{0}(\mathbf{x}),\]  
where \(\zeta_{2}\) is a random number sampled in from the uniform distribution in \([0,1]\). The interested reader is referred to [16] for more details about SSA. It is worth noting that the vector \(\mathbf{c}=(c_{1},\ldots,c_{M})\) of stochastic constants governs the behavior of the system, since it affects the calculation of the propensity functions of the reactions. Unfortunately, the values in \(\mathbf{c}\) are hard (or even impossible) to measure with specific lab experiments; nonetheless, such fundamental parameters can be estimated by exploiting computational intelligence methods [19].