### Model signaling pathway: simplified MAPK pathway

Utilizing the model and numerical simulation method described above in "Method" section, we simulated the activity of the simplified MAPK signaling pathway (Fig. 2), comprising 11 proteins and 14 interaction weights. In the model pathway, activated ERK is translocated across the nucleus membrane to regulate cell proliferation, migration, and anti-apoptosis effects [30]. These processes controlled by ERK are considered to be cell phenotypes in our MAPK signaling pathway. It is worth noting that all of the weights in the current model are assumed to be non-negative, representing positive interactions only. We chose this simplified MAPK pathway structure to test the current model, particularly the effect of the stochastic process on signaling transduction.

We take GRB2 as an example target protein to further explain the modeling process, as shown in Fig. 2. The GRB2, denoted by \(x_4\), is activated by three immediately upstream proteins: EGFR (\(x_1\)), ERBB2 (\(x_2\)), and MET (\(x_3\)). GRB2 can stimulate only one downstream protein, SOS (\(x_6\)). EGFR (\(x_{1}\)) has two weights, \(w_{11}\), and \(w_{12}\); ERBB2 (\(x_{2}\)) has two weights, \((w_{21}\) and \(w_{22})\); MET (\(x_{3}\)) has two weights, \((w_{31}\) and \(w_{32})\); and GRB2 \((x_{4})\) has the only one weight, \((w_{41})\), for stimulating SOS. Therefore, the activities of upstream and downstream molecules of GRB2 can be expressed as follows:

$$\begin{aligned} {\text{upstream of GRB2:}} & \;\left( {r\min \left( {w_{{11}} ,w_{{12}} } \right) + \left( {1 - r} \right)\max \left( {w_{{11}} ,w_{{12}} } \right)} \right)x_{1}, \\ & + \;\left( {r\min \left( {w_{{21}} ,w_{{22}} } \right) + \left( {1 - r} \right)\max \left( {w_{{21}} ,w_{{22}} } \right)} \right)x_{2}, \\ & + \;\left( {r\min \left( {w_{{31}} ,w_{{32}} } \right) + \left( {1 - r} \right)\max \left( {w_{{31}} ,w_{{32}} } \right)} \right)x_{3}, \\ {\text{downstream of GRB2:}} & \;\left( {r\min \left( {w_{{41}} } \right) + \left( {1 - r} \right)\max \left( {w_{{41}} } \right)} \right)x_{4}. \\ \end{aligned}$$

The self-activity of GRB2 occurred by stochastic fluctuation can be expressed as \(\sigma x_{4}(t)dW(t)\), where \(W(t)\approx N(0,\Delta t)\).

### Single stimulus simulation

We first simulated the pathway upon receiving a single stimulus. We initially simulated the temporal evolution of proteins in the pathways without stochastic self-activity and then added the stochastic activity effect in subsequent simulations.

We consider a case where the receptor EGFR is instantly activated by EGF initially (\(x_i(0)\) = 1, where \(x_i\) indicates EGFR). The activity of EGFR then decreases monotonically since no further EGF stimulation is applied. The activity of each downstream of EGFR increases sequentially according to the order of the MAPK signaling pathway (depicted in the dotted rectangle range in Fig. 3). The activity of each protein then sequentially decreases following the decrease in EGFR activity. The inactivation rates of all proteins were extremely slow (depicted as a solid rectangle in Fig. 3).

Figure 4 shows the trajectories of the activity of each protein in the MAPK signaling pathway when the stochastic self-activity was included. We set \(\sigma = 1.0\) to solely investigate the influence of the stochastic process, *W*(*t*). Following EGFR activity, the proteins in the MAPK signaling pathway are sequentially activated (depicted in the solid rectangle range in Fig. 4). GRB2, downstream of EGFR, activates SOS, which in turn increases the activity of RAS, and so on down the pathway. The signal appeared to be amplified as it moved downstream in the pathway. Interestingly, the protein shifted to an inactive state (protein activity = 0) significantly faster than observed in the non-stochastic case (Fig. 4 vs. Fig. 3). We compared the decay time of all proteins in non-stochastic and stochastic cases (Additional file 1: Fig. S1). Since the protein activities did not decrease to near 0 in a non-stochastic case, we chose a specific time point, the time to reach 10% of the maximum value (*p*-value < 0.05, Student t-test). Although the random variables of the stochastic process *W*(*t*) were sufficiently small due to \(\Delta t = 0.001\) and \(W(t)\approx N(0,\Delta t)\), its impact on the MAPK pathway was nevertheless significant.

### Influence of weight in the MEK–ERK interaction

We next examined the influence of weight on protein activity. We focused on the weight of the MEK–ERK interaction to ignore the effect of a downstream protein affecting another protein. Figure 5 shows the trajectories of the activities of MEK and ERK. We set all weights to 0.5 except for the weight of the MEK–ERK interaction. We performed the simulation in response to a single stimulus assuming three different weights of MEK and ERK: 0.25, 0.5, and 0.75 (Fig. 5). In all cases, ERK activity was sustained over a longer period (Fig. 5), and there was a slower decrease in ERK (black) than MEK (cyan) activity. When the weight of the MEK–ERK interaction was smaller than the weights of other interactions, the initial activity of ERK was lower than that of MEK (Fig. 5a), although ERK activity caught up to reach the level of MEK activity at around time step 30,000. When the weight was greater than or equal to the weights of the other interactions, the ERK activity was maintained at a higher level for the entire period of the simulation (Fig. 5b, c).

### Multiple simulations

To determine whether the stochastic effect on target protein activity significantly changes the overall behavior of the MAPK pathway, we conducted 1000 simulations for a single stimulus in which different random weights were assigned following a uniform distribution \(w_{ij}^s \approx U[0,1]\), where \(1 \le s \le 1000\) refers to each simulation trial. The temporal evolution of each protein appeared to be similar in all simulations (Fig. 6): all proteins were sequentially activated and then inactivated (converging to zero). We observed amplified signaling activity as the signal traveled downstream in the pathway (e.g., from EGFR to MEK), consistent with the observations for the single-simulation case (Fig. 4). Interestingly, the noise of protein activity also appeared to be increased when moving downstream in the pathway (RAF, MEK, ERK) compared to that upstream (EGFR, GRB2, SOS, and RAS) (Fig. 6, with a larger standard variation in (e) and (f), and Fig. 7.)

We also conducted additional simulations with different parameter sets of \(k_{u},k_{d}\), and \(\sigma\) (Additional file 1: Figs. S2–S8). Comparing with Figs. 6, 7, simulation results are qualitatively consistent. In particular, signaling amplification occurs in all parameter sets, although the amplitude of each protein activity rather changes with parameter sets. The smaller \(k_d\) that regulates the downstream magnitude, the greater the target protein activation. The smaller the \(k_u\), which governs the magnitude of the upstream, the smaller the activation of the target protein. A smaller stochastic factor (\(\sigma\)) delays inactivation.

The parameter *r* is a coefficient of restriction in the Waller–Kraft operator. A value close to one imposes strong restrictions similar to the Boolean logic“AND”operation; on the other hand, a small value close to zero indicates less restriction, similar to a“OR”operation. In this study, we choose the “AND” logic of the Waller–Kraft operator to perform numerical simulations because the signal transduction between proteins simultaneously occurs within a sufficiently fast time (\(10^{-9}\)(ns)–\(10^{-3}\)(ms)). Even though we chose the “AND” logic, we compared the difference between “AND” and “OR” logic operators (Additional file 1: Figs. S9–12). It is worth noting that the simulation results with some large *r* values (e.g., \(r \ge 0.75\), “AND” logic) did not show a significant difference. The results with “AND” logic are significantly different with some small *r* values (e.g., \(r \le 0.25\), “OR” logic).

We next compared overall protein activities from stimulation to an inactive state. For each protein, we first calculated the area under the curve of the temporal profile. We then calculated the average over the areas and the coefficient of variation (Fig. 8). The coefficient of variation (CV) represents the ratio to the average value of the standard deviation; that is, if the sample standard deviation, *s*, and mean, \(\mu\), then the CV is

$$\begin{aligned} \text {CV} = \frac{s}{\mu }. \end{aligned}$$

The actual value of the CV is dimensionless because it is independent of the unit in which the measurement was performed. It is useful to use coefficients of variation instead of standard deviations to compare data sets with different units or means. We observed an increasing tendency of overall protein activities from GRB2 to ERK. As the activity of the proteins increased, the range of the coefficient of variation due to the stochastic process also increased.

### Repeated stimulation of all receptors

Finally, we simulated the signaling pathway when the three receptors (EGFR, ERBB2, MET) were stimulated repeatedly. We conducted 1000 simulations while varying weights between proteins in the pathway. We set the receptors are fully activated initially (\(x_{i}(0) = 1\), where \(x_i\) stands for EGFR, ERBB2, or MET). After activation, the activity of the three proteins gradually decayed, returning to almost the inactive state (protein activity near zero). An additional stimulus was then applied to each receptor. Since the decrease rate of the activity differed for each receptor, the timing of additional activation also varied among the receptors (Fig. 9), resulting in different peak times for EGFR, ERBB2, and MET.

We observed various temporal evolution patterns of MAPK signaling pathway activity. From the representative temporal changes of the proteins, we observed that the activities of receptors increased quickly upon stimulation but also decreased rapidly unless additional stimulation was applied. The activity of GRB2 exhibited a time-delayed increase and a slower decrease compared with those of its upstream proteins (a representative example of the temporal evolution of protein activities is shown in Fig. 9a). Stimulation by additional EGFR activation induced a spike in GRB2 activity around the time step around 4000. Following additional stimulation by ERBB2 and MET around time step 4000 and 5000, GRB2 was again reactivated and subsequently decayed rather quickly (Fig. 9a). Interestingly, the activity of GRB2 exceeded one because ERBB2 and MET stimulated GRB2 almost simultaneously. To compare the timing of maximum activation, we compared the peak time (Fig. 9b) in the time interval [6000, 9000]. The peak time of the two receptors almost perfectly coincided with the peak time of GRB2.

Interestingly, the average activities of MEK and ERK were maintained at a high level throughout the simulation time, even though the receptors were fully inactivated (Fig. 10a). Note that even though the stimulus repeatedly occurred, the average protein activity was maintained in a narrow range. Likewise, we compared overall protein activities over a fixed time scale, demonstrating an increasing tendency of overall protein activities from GRB2 to ERK. The range of coefficient of variation also appeared to be amplified, implying that the noise generated by the stochastic process propagates through the pathway when the signals travel downstream (Fig. 10b).