Kullback–Leibler Divergence of an Open-Queuing Network of a Cell-Signal-Transduction Cascade

Queuing networks (QNs) are essential models in operations research, with applications in cloud computing and healthcare systems. However, few studies have analyzed the cell’s biological signal transduction using QN theory. This study entailed the modeling of signal transduction as an open Jackson’s QN (JQN) to theoretically determine cell signal transduction, under the assumption that the signal mediator queues in the cytoplasm, and the mediator is exchanged from one signaling molecule to another through interactions between the signaling molecules. Each signaling molecule was regarded as a network node in the JQN. The JQN Kullback–Leibler divergence (KLD) was defined using the ratio of the queuing time (λ) to the exchange time (μ), λ/μ. The mitogen-activated protein kinase (MAPK) signal-cascade model was applied, and the KLD rate per signal-transduction-period was shown to be conserved when the KLD was maximized. Our experimental study on MAPK cascade supported this conclusion. This result is similar to the entropy-rate conservation of chemical kinetics and entropy coding reported in our previous studies. Thus, JQN can be used as a novel framework to analyze signal transduction.


Introduction
A cell-signal-transduction cascade forms a network of reactions that modify protein molecules. The modified protein can diffuse in the cytoplasm and modify another protein.
Finally, modified proteins translocate to the nucleus to bind to the promoter in DNA to promote the gene expression encoded in the DNA [1]. Currently, kinetic models based on partial-differential equations (PDEs), including diffusion terms, primarily simulate signal-transduction kinetics and networks. Although PDEs reduce theoretical difficulties through mathematical simplification, the stochastic factors essential to biological kinetics are generally omitted.
The queue in operations is used for analyzing the congestion phenomenon of stochastic systems. For example, queuing theory applies to service-processing design to maintain the waiting time within a specific range [2][3][4][5][6]. The queuing theory applies to a system consisting of a "server", a "waiting room", and a "customer" who arrives and stays for a certain time and accepts the "service." Kendall introduced the symbol of the queue in the form of A/B/C/D (A: the arrival of the customer, B: the distribution of the service time, C: the number of servers, D: the capacity of the system including the waiting room) [7]. Recently, this theory has been applied to analyzing the queuing of packets in routers in communication networks to model gene-expression networks [8]. Furthermore, the theory can apply to the management of the healthcare system [9] in which the patient-waiting time is analyzed when introducing electronic medical-record systems. COVID-19 pandemic infections have been modeled using this theory [10,11]. On the other hand, there have been (ATP). The signal transduction forms a network of sequential chain reactions that modify signaling molecules, a i , which constitute the network nodes. Figure 1 illustrates a modeled signal-transduction network. For simplification, the signal is mediated by a mediator P. The reaction of signaling molecules, A i , in the ith node (1 ≤ i ≤ n) is expressed as follows: where P diffuses from the outside and arrives at an arrival rate of λ oi at the ith node, a i , where protein kinase A i stays (Equation (1)). λ i is the arrival rate of A i−1 -P at the ith node (Equation (2)), where the signal can be transduced by the exchange of P from A i−1 to A i at a rate of µ i (Equation (3)).
Entropy 2023, 25, x FOR PEER REVIEW 4 of 13 node in the cascade. The arrival rate of P at each node i is λi from the outside of QN according to the Poisson process, and it is possible for P to leave the network from any node; (2) there is a customer P that is served (transferred to Ai) in the network; (3) P in the network is sufficiently high; (4) The service time follows an exponential distribution at all nodes; and (5) The service code is the arrival order of service (FCFS: first-come firstserved). In equilibrium state, given that the number of guests arriving per unit time to a node is equal to the number of departures per unit time from that node, we get the following equation, which is called the traffic equation: λi (0 ≤ i ≤ 6) represents the arrival rate of P from node ai. λoi denotes the arrival rate of P from the outside. The ladder-like mark symbolizes the queue of the signaling molecule. µi inside the circle represents the exchange ("service" in the queueing theory) rate in which P carried by ai−1 is transferred to another signaling molecule ai at the ith node. The arrow represents the orientation of irreversible signal transduction. The final node is DNA (binding histone proteins), where the transcriptional reaction occurs [7]. This is depicted in a figure, and schemes follow the same formatting. The symbol ai (Ai) above the µi represents that signal molecule Ai stays at the node.
The vertical axis represents the ratio of the phosphorylated signaling molecule ai concentration (1 ≤ i ≤ 3) to that at the steady state, ai (t)/ai = ρi (t)/ρi, where ai represents the signaling molecule Aj-P concentration at the steady state. The horizontal axis indicates the duration, where τi (1 ≤ i ≤ 3) represents the phosphorylation period for each signaling molecule ai. In this case, the signal transduction proceeds in the order of 1→2→3.
Furthermore, we set ρi = λi/µi. ρi denotes the value at the initiation of the signal event. As mentioned below, ρi is equal to ai. The following expression can be obtained using Jackson's theorem [20] and the probability that ai stays at the ith node at the steady state: and where = Since the arrival rate is sufficiently low in the entire signal-transduction cascade, λi/µi << 1 holds for all ai. In this case, from Equation (9), In the actual signaling system, the individual concentrations of the signaling molecules, A i (1 ≤ i ≤ n), are sufficiently low; a single signaling event rarely occurs. Therefore, I hypothesized that under the following three assumptions, the signaling event occurs by a Markov process: (i) the probability that a signal event will occur is constant, (ii) the probability that an ith step (node) in the cascade occurs at any time interval (t, t + ∆) does not depend on the number of events before time t where ∆ signifies the minimal interval, and (iii) the probability of an event occurring twice or more during a minute's time ∆ is negligible and is denoted as o(∆). Let the likelihood p i in the ith step (nodes) in the cascade be the probability that an event occurs A i -P (=a i ) times in time t to t + ∆, and let λ i ∆ be the probability that an event occurs once in a minute, where ∆ denotes a stable arrival rate. The arrival interval of p i becomes a random variable sequence that follows the same exponential distribution, and the arrival interval of the signal components is independently and identically distributed: Transformation of Equation (4) gives: Approaching ∆ → 0 gives a differential equation: Solving Equation (6) gives the solution that p i is reached by a Poisson distribution: Entropy 2023, 25, 326 4 of 12 Equation (7) describes a queuing model in which A i-1 -P arrival rate can be formulated as a Poisson process. Here, the service was defined as the time taken for transcription. In the queue, in accordance with the Poisson process, λ i is termed as the arrival rate.
The example of JQN in shown ( Figure 1). Thus, P is transferred by A i through the network, and the endpoint is the transcription of the nuclear DNA. The diffusion of P carried by A i is the rate-limiting step determining the signal-transduction rate. In this case, P is known as a "customer" because it comes from the outside. The exchange rate of P from A j-1 -P to A i , µ i , is known as the "service rate". In this way, we assumed that signaling molecules form an open-queue network known as JQN.
Here, the QN of signal transduction satisfies the following conditions of Jackson network: (1) There are phosphorylation "servers" A i -P with the exchange rate µ i at the ith node in the cascade. The arrival rate of P at each node i is λ i from the outside of QN according to the Poisson process, and it is possible for P to leave the network from any node; (2) there is a customer P that is served (transferred to A i ) in the network; (3) P in the network is sufficiently high; (4) The service time follows an exponential distribution at all nodes; and (5) The service code is the arrival order of service (FCFS: first-come first-served). In equilibrium state, given that the number of guests arriving per unit time to a node is equal to the number of departures per unit time from that node, we get the following equation, which is called the traffic equation: represents the arrival rate of P from node a i . λ oi denotes the arrival rate of P from the outside. The ladder-like mark symbolizes the queue of the signaling molecule. µ i inside the circle represents the exchange ("service" in the queueing theory) rate in which P carried by a i−1 is transferred to another signaling molecule a i at the ith node. The arrow represents the orientation of irreversible signal transduction. The final node is DNA (binding histone proteins), where the transcriptional reaction occurs [7]. This is depicted in a figure, and schemes follow the same formatting. The symbol a i (A i ) above the µ i represents that signal molecule A i stays at the node.
Because the signal-molecule arrival rate is sufficiently low, we assumed that the ith molecule, A i , is localized at the ith node. Here, we set {a} = {a 1 , a 2 , . . . . A n } = {a i }, which represents the JQN-node state at the steady state and {a(t)} = {a 1 (t), a 2 (t), . . . . a (t)} = {a i (t)}, which represents the JQN-node state during the signal transduction.
The vertical axis represents the ratio of the phosphorylated signaling molecule a i concentration (1 ≤ i ≤ 3) to that at the steady state, a i (t)/a i = ρ i (t)/ρ i , where a i represents the signaling molecule A j -P concentration at the steady state. The horizontal axis indicates the duration, where τ i (1 ≤ i ≤ 3) represents the phosphorylation period for each signaling molecule a i . In this case, the signal transduction proceeds in the order of 1→2→3.
Furthermore, we set ρ i = λ i /µ i . ρ i denotes the value at the initiation of the signal event. As mentioned below, ρ i is equal to a i . The following expression can be obtained using Jackson's theorem [20] and the probability that a i stays at the ith node at the steady state: and where Since the arrival rate is sufficiently low in the entire signal-transduction cascade, λ i /µ i << 1 holds for all a i . In this case, from Equation (9), The queue length of the ith node, L i , is equal to that of the independent queue. In this case, L i and the sum of L i , L, are expressed by Little's formula [24] as follows: where λ i /µ i << 1, Equations (13) and (14) were used. Likewise, we set the concentration of a i during the signal transduction: and An example of the time courses of signal transduction expressed by ρ i (t)/ρ i is shown in Figure 2. and The queue length of the ith node, Li, is equal to that of the independent queue. In this case, Li and the sum of Li, L, are expressed by Little's formula [24] as follows: where λi/µi << 1, Equations (13) and (14) were used. Likewise, we set the concentration of ai during the signal transduction: and An example of the time courses of signal transduction expressed by ρi (t)/ρi is shown in Figure 2.

Kullback-Leibler Divergence of JQN
Suppose that the signal-network activation in which the ith node is activated is obtained. If the complete signal-transduction process is repeated a times in the entire signal-transduction system and each signal-transduction process is repeated ai times in the ith node, the probability that the transduction system stays at the given JQN state {a} is calculated as follows: Taking the logarithm of the left-hand side, we obtain (Appendix A): Here,

Kullback-Leibler Divergence of JQN
Suppose that the signal-network activation in which the ith node is activated is obtained. If the complete signal-transduction process is repeated a times in the entire signaltransduction system and each signal-transduction process is repeated a i times in the ith node, the probability that the transduction system stays at the given JQN state {a} is calculated as follows: a! a 1 (t)!a 2 (t)! . . . a n (t)! p(a) = a! a 1 !a 2 ! . . . a n ! ρ 1 a 1 ρ 2 a 2 · · · ρ n a n Taking the logarithm of the left-hand side, we obtain (Appendix A): log a! a 1 (t)!a 2 (t)!...a n (t)! ρ 1 a 1 (t) ρ 2 a 2 (t) · · · ρ n a n (t) Here, which has the form of the KLD. The KLD has been used in sensor and imaging analyses [25], Bayesian model diagnostics [26], clinical trial assessments [27], and information science [28]. We recently performed the EGFR cascade analysis based on the KLD theoretical framework [23]. Furthermore, from Little's law in queueing theory [24], In Equation (21), τ oi denotes the duration of a single signal mediator P for a single signaling molecule a i. We hypothesized that cells signal more efficiently by eliminating redundancies in terms of energy metabolism. This assumption was also adopted in our previous study [23] and led to maximizing KLD, i.e., information gain, per the signal-event time at each step of the cell-signal-transduction cascade.
To maximize D(ρ(t)||ρ) per a given duration τ, we introduced function F with the constraint Equation (20), as given below.
where γ is an arbitrary parameter independent of node i. From Equation (22), Subsequently, Using Equation (21), Accordingly, taking the sum of both sides, Therefore, we have The negative value of γ indicates the average KLD per signal duration. In conclusion, the average KLD rate per signal duration is independent of the signal step node because γ is independent of the step-node number, i.

Chemical Potential in JQN
Here, we defined the chemical potential of a i at the pre-signal event state as and the signal state as where k B denotes the Boltzmann coefficient and T represents the temperature of the system. The thermodynamic entropy changes, ∆s i , between the pre-signal and signal event status are given by Therefore, the entropy changes can be related to KLD [29][30][31]. From Equation (30), we can obtain the relationship between the thermodynamic entropy change and KLD: Accordingly, γ has the dimension of a negative entropy-production rate.

Application of KLD Theory to Signal-Transduction Analysis
The mitogen-activated protein kinase (MAPK) signaling cascade (consisting of ASK1, MKK4, JNK, and HSF1) induced by cell stress has been analyzed [32][33][34][35][36][37][38]. The MAPK cascades include the ERK pathway, p38 pathway, and JNK cascade, which are associated with the MAPK stress-response pathway. This pathway is activated by a range of extracellular stress stimuli, including heat shock, oxidation, and exposure to ultraviolet light or radiation [39][40][41]. A MAPK signal cascade in which signal transduction is represented as a sequential activation of the signaling molecules can be described as follows (see Figure 1): where phospho represents the phosphorylated status of signaling molecules, and P denotes the released phosphate and adenosine diphosphate, respectively. In this cascade, phosphate is the mediator. With reference to Equation (32), the average KLD rate of the ASK1-MKK4-JNK-HSF1 cascade is calculated as follows [29]: In the above equation, we used ρ i (t) = a i (t) and ρ i = a i . The integrals in the third term in Equation (33) were calculated using the integral of the plot (Figure 3) (0 ≤ t ≤ τ i ). In the experimental study, we cultured the A431 skin-cancer cell line in a serum-free medium and stimulated it with EGFR. In the starved A431 cells, the MAPK-related signal cascade (ASK1-MKK4-JNK-HSF1) was tentatively activated with starvation and minimal EGFR stimuli with EGF. As a result, the KLD rates were around 3.0, which are similar to each other in terms of the Cohen's factor d ≤ 0.5 (Table 1), indicating that Equations (25)

Discussion
A primary theoretical approach for analyzing signal transduction in systems biology is the kinetics consisting of continuous differential equations for the time-evolution of the concentration of signaling molecules. Based on some experimental facts, a further numerical simulation was performed by substituting measurable-reaction kinetic coefficients and other parameters. However, in general, the number of signaling molecule proteins in the cytoplasm was small, so concentration fluctuations were expected to be large. Therefore, an analysis using differential equations assuming continuous variables is an appropriate approximation method. In practice, to think of the dynamics of the signaling molecule as a discrete quantity, it is necessary to use the dynamics of discrete parameters, i.e., the stochastic dynamics framework. However, such kinetics are complex

Discussion
A primary theoretical approach for analyzing signal transduction in systems biology is the kinetics consisting of continuous differential equations for the time-evolution of the concentration of signaling molecules. Based on some experimental facts, a further numerical simulation was performed by substituting measurable-reaction kinetic coefficients and other parameters. However, in general, the number of signaling molecule proteins in the cytoplasm was small, so concentration fluctuations were expected to be large. Therefore, an analysis using differential equations assuming continuous variables is an appropriate approximation method. In practice, to think of the dynamics of the signaling molecule as a discrete quantity, it is necessary to use the dynamics of discrete parameters, i.e., the stochastic dynamics framework. However, such kinetics are complex and inevitably challenging to apply to complex systems such as cell signaling, which involve non-linear interactions of many signaling molecules, as is the case in systems biology. In contrast, QN theory has a simple framework that can analyze the reaction of discrete quantities. Furthermore, the queues are one-directional, which is consistent with the one-way signaling mechanism in the intracellular biochemical-reaction network, in which the intracellular information is transmitted through the processing of information in the cell membrane or cytoplasm to the nucleus.
If the customer is the signaling molecule itself, a closed QN can be applied where the concentration of the signaling molecule is constant during the signaling period. However, in this case, there are multiple customers, and the QN is more complicated. Besides, each signaling molecule concentration is unstable owing to the fluctuation as aforementioned. To understand the essence of signal transduction and simplify the model, the signal mediators ATP or phosphate were treated as customers. Since these molecules are in large quantities in cells and are continuously and constantly supplied to the signal-transduction system from the outside, we applied an open QN. The open QN, where customers arrive externally according to the Poisson process, better reflects the intracellular reaction.
In the application of JQN, when ρ i < 1 holds for all ρ i (1 ≤ i ≤ n), the solution for the probability at the equilibrium distribution is given as a product-form solution. This solution was the basis of this study, shown in Equation (9), and the KLD of the entire system was definable in Equation (19). Since the stationary distribution was in product form and the lengths of the queues of each node at any given time were independent of each other, the process of leaving each node became a Poisson process. It was the first significant development in QN theory, and applying Jackson's theorem to search for similar product-form solutions in other networks has been the subject of much research. Essentially, a cell-signaling network consists of a considerable number of nodes or phosphorylated (the active form of) signaling molecules, where each service rate (phosphorylation of other signaling molecules) has different values. Thus, JQN represents the signal-transduction system. We applied Little's formula [2] and Jackson's theorem [20] to queuing theory. We obtained Equations (26) and (27) in which the JQN KLD per transduction time is constant when the KLD is maximized. As a result, the conclusion was drawn that the KLD rate is constant regardless of the node. Intracellular signal transduction is a complicated process consisting of a series of modifications and demodification reactions of intracellular proteins. Due to the simplicity of this chain, we introduced the Markov chain (M), i.e., the assumption that the reaction step depends only on the preceding and succeeding steps and that the reaction process formed a model of a queue of M/M/A i . Furthermore, it was necessary to add thermodynamic information to the model as signal transduction is interpreted as the simultaneous transduction and conversion of biological information. In order to understand the latter, it was necessary to think about fluctuation theorem and entropy coding, and KLD represents the information gain before and after cell stimulation. These concepts are closely related to each other, forming a unified framework of an information theory of intracellular signal transduction. As a result, in this study, an important conclusion was obtained: the KLD per transduction time (code length) in the same cascade signal-transduction step is constant. Importantly, this conclusion was verified by measuring the modification and demodification of proteins before and after cell stimulation.
By introducing quantitative evaluation in this research, it was possible to estimate the transduction time in the signal-conversion step from the fluctuation before and after the whole signal network. Additionally, based on the result that the KLD rate is constant in each step in the known cascade, conversely, we may find a new cascade by comparing the KLD rate of each signaling molecule. For example, when KLD rates of signaling molecules A and B are similar in the given signal event, a signal cascade may be formed between A and B. Further experimental research will be required to validate this idea.
The model of the QN is valid for the following reasons. First, it facilitates the consideration of signal transformations that are irreversible reactions. In this model, the rate-limiting step is a long, reversible process in which a signaling molecule diffuses to the next signaling molecule that is its substrate. The phosphorylation process in service is a short process, which is irreversible in this model. Therefore, the signal never passes through the same node again without considering the feedback mechanism. Such a model framework is useful for signal transduction. Second, protein-protein phosphorylation is a complex non-linear reaction that is difficult to target for kinetic analysis. Therefore, in the application of network theory, we can ignore elementary processes that can be analyzed with simpler models. Third, protein molecules have large concentration fluctuations and are not subject to simple kinetics. This framework overcomes these problems and enables a unified analysis.
We focused on a stress-induced ASK1-MKK4-JNK-HSF1 signal-transduction cascade and found that the average information gain, i.e., KLD, per signal duration is consistent in the pathway, indicating the conservation of KLD rate. Conversely, the cell signaltransduction system proceeds in circumstances in which the rate of signal transduction is maximized. By observing the temporal variation of the phosphorylation of individual signaling molecules, plotting the increase against time, and taking the logarithm of the integral value, the KLD increase rate per hour can be comprehensively measured. As a result, signals can be transmitted between signaling molecules with similar KLD growth rates. This is one of the methods that makes it possible to discover unknown cascades between so many signaling molecules. Furthermore, by changing the dose of radiation that produces reactive oxygen, it is possible to identify signaling molecules that show an increase in the KLD rate correlated with the radiation dose and to more accurately identify cascades specific to the type of stress.
Previously, we applied entropy encoding to a code-string model of the signal cascade and found that the entropy-production rate is conserved [42], which corresponds with the conclusion of this study where the KLD rate for the queue is consistent regardless of the node, a i . Also, we reported that the non-equilibrium kinetics of signal transduction indicates that the entropy-production rate is consistent in the signal-transduction cascade [22,23]. According to Equations (30) and (31), KLD can be expressed by a chemical potential change and thermodynamic entropy ∆S in the entire signal transduction. In summary, the entropy-production rate conservation holds in the JQN, entropy-coding theory, and nonequilibrium thermodynamics. This indicates an intriguing relationship between queuing, entropy encoding, and chemical potential.
A limitation of the present study is that we used continuous derivative operations in deriving Equations (26) and (27) rather than discrete operations; therefore, this mathematical evaluation should be further verified. A detailed study, replacing the derivative operation with a discrete operation, may be explored in future work.
In conclusion, the JQN can be used for a quantitative analysis of signal transduction.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Calculation for Equation (18) in the Text
The following is the calculation deriving for Equation (18).