Mathematical model of the cell signaling pathway based on the extended Boolean network model with a stochastic process

Background In cell signaling pathways, proteins interact with each other to determine cell fate in response to either cell-extrinsic (micro-environmental) or intrinsic cues. One of the well-studied pathways, the mitogen-activated protein kinase (MAPK) signaling pathway, regulates cell processes such as differentiation, proliferation, apoptosis, and survival in response to various micro-environmental stimuli in eukaryotes. Upon micro-environmental stimulus, receptors on the cell membrane become activated. Activated receptors initiate a cascade of protein activation in the MAPK pathway. This activation involves protein binding, creating scaffold proteins, which are known to facilitate effective MAPK signaling transduction. Results This paper presents a novel mathematical model of a cell signaling pathway coordinated by protein scaffolding. The model is based on the extended Boolean network approach with stochastic processes. Protein production or decay in a cell was modeled considering the stochastic process, whereas the protein–protein interactions were modeled based on the extended Boolean network approach. Our model fills a gap in the binary set applied to previous models. The model simultaneously considers the stochastic process directly. Using the model, we simulated a simplified mitogen-activated protein kinase (MAPK) signaling pathway upon stimulation of both a single receptor at the initial time and multiple receptors at several time points. Our simulations showed that the signal is amplified as it travels down to the pathway from the receptor, generating substantially amplified downstream ERK activity. The noise generated by the stochastic process of protein self-activity in the model was also amplified as the signaling propagated through the pathway. Conclusions The signaling transduction in a simplified MAPK signaling pathway could be explained by a mathematical model based on the extended Boolean network model with a stochastic process. The model simulations demonstrated signaling amplifications when it travels downstream, which was already observed in experimental settings. We also highlight the importance of stochastic activity in regulating protein inactivation. Supplementary Information The online version contains supplementary material available at 10.1186/s12859-022-05077-z.


Ito-Taylor expansion
We briefly explain how we derived the numerical integration scheme for our stochastic differential equation from Ito-Taylor expansion [1].
First, we call a stochastic differential equation (SDE): where f and g are sufficiently smooth. Let F be the continuously differentiable function, then from Ito's lemma, we have ∂F (X(t)) ∂X ∂F (X(t)) ∂X dW (t). (3) Next, we define operators as follows: Then, Eq.
Taking integral, Eq. (5) is derived as: where R is the remaining terms which include the double integral terms: Selecting Eq. (4), we obtain where new remainderR. Note that the double integral is evaluated to be t t 0 Substituting Eq. (9) into Eq. (8), we can produce the numerical integration scheme for the SDE from Ito-Taylor expansion with a time discretization 0 = t 0 < t 1 < · · · < t N = T of a time interval [0, T ] as follows: where △t = t i+1 − t i and △W i = W i+1 − W i for i = 0, 1, 2, · · · , N − 1 with the initial condition X(t 0 ) = X 0 . The random variable △W i are independent N (0, △t) normally distributed random variables.

Algorithms
To explain the algorithms, we take an example of GRB2 and GRB2 upstream proteins such as EGFR, ERBB2, and MET (Fig.2). First, we calculate each Boolean value of EGFR, ERBB2, and MET based on a Boolean logic function (defined as void Boolean gate X).
# Function o f X (X: EGFR, ERBB2, MET) void B o o l e a n g a t e X ( · · ·){ int j ; Second, we estimate upstream flowing into GRB2 using upstream function (defined as void upstream gate grb2). The effect of stochastic factors on protein inactivation We compared inactivation time for all proteins in non-stochastic and stochastic cases (See Figures 3-7). Because the protein activities did not decrease to 0 in a non-stochastic case, we compare the time to decrease to a certain value (e.g., 10% of the maximum). The protein inactivation in the non-stochastic case is significantly delayed compared to the inactivation in the stochastic case.  Figure S1: The time to decay of all proteins in stochastic (left) and non-stochastic (right) cases. p-value < 0.05.

Effect of model parameters
We performed additional simulations with different sets of parameter sets of k u , k d , σ, and r. First, we present an average protein activity over 1000 simulations when the three parameters are set to be one (See Figures 6-7 and Figure S2). Effect of parameters k d and k u Figure S3: Temporal evolution of average protein activity over 1000 simulations when the downstream coefficient is decreased to 50% compared to Fig. S2.
To compare the effect of the downstream coefficient (k d ) on the signaling pathway, we considered two cases of k d = 0.5 (Fig. S3) and k d = 0.1 (Fig. S4). A decreased downstream coefficient makes the magnitude of each target protein activation larger. Although the signaling magnitude is different from Fig. S2, the order of protein amplification is maintained. In other words, the protein activity is sequentially amplified and sustained longer as  Similarly, we performed 1000 simulations with a different upstream coefficient k u (=0.5,0.1) in Fig. S5-Fig. S6. The upstream coefficient (k u ) affects the magnitude of the protein activation. A smaller coefficient makes the target protein activity smaller since an activation rate by an upstream protein Figure S6: Temporal evolution of average protein activity over 1000 simulations when the upstream coefficient is decreased to 10% compared to Fig. S2.
of the target protein is smaller. The changes in the upstream coefficient did not affect the decay of protein activation nor the amplification of protein activity.  We considered two different stochastic coefficients σ = 0.5, 0.1. A smaller stochastic coefficient did not affect the magnitude of protein activation but delayed the protein inactivation ( Fig. S7 -Fig. S8).

Effect of parameter r
To evaluate the effect of the parameter r, we considered four different values of r (r = 0.0, 0.25, 0.75 and r = 1.0). Figure S9 illustrates the temporal evolution of all proteins for each parameter r. Suppose the r value is one in the Boolean function, a minimum weight among the weights between the target protein and upstream proteins will be selected. Thus, a little activation of a target protein is expected by upstream proteins. Similarly, a minimum weight between a target and downstream proteins is selected, leading to a slower decay of a target protein. To quantify the effect of the parameter r on the signaling transduction, we calculated the time to almost full inactivation (e.g., a time point of each protein to reach less than 0.000001). We compare the time to this full inactivation of all proteins. The full inactivation time assuming r = 1.0 is not statistically different from the time assuming r = 0.75 (Fig. S10, Fig. S9a vs Fig. S9b). However, the time is significantly different from r = 0.0 (Fig. S11, Fig. S9a vs Fig. S9d). Lastly, the inactivation time assuming r = 0.75 is significantly different from r = 0.25 (Fig. S12, Fig. S9b