A Model of Cell Biological Signaling Predicts a Phase Transition of Signaling and Provides Mathematical Formulae

A biological signal is transmitted by interactions between signaling molecules in the cell. To date, there have been extensive studies regarding signaling pathways using numerical simulation of kinetic equations that are based on equations of continuity and Fick’s law. To obtain a mathematical formulation of cell signaling, we propose a stability kinetic model of cell biological signaling of a simple two-parameter model based on the kinetics of the diffusion-limiting step. In the present model, the signaling is regulated by the binding of a cofactor, such as ATP. Non-linearity of the kinetics is given by the diffusion fluctuation in the interaction between signaling molecules, which is different from previous works that hypothesized autocatalytic reactions. Numerical simulations showed the presence of a critical concentration of the cofactor beyond which the cell signaling molecule concentration is altered in a chaos-like oscillation with frequency, which is similar to a discontinuous phase transition in physics. Notably, we found that the frequency is given by the logarithm function of the difference of the outside cofactor concentration from the critical concentration. This implies that the outside alteration of the cofactor concentration is transformed into the oscillatory alteration of cell inner signaling. Further, mathematical stability kinetic analysis predicted a discontinuous dynamic phase transition in the critical state at which the cofactor concentration is equivalent to the critical concentration. In conclusion, the present model illustrates a unique feature of cell signaling, and the stability analysis may provide an analytical framework of the cell signaling system and a novel formulation of biological signaling.

A transient binding cofactor such as ATP/GTP, or phosphorylation of amino acid residues, controls signaling molecule interactions and the subsequent modification of signaling molecules. Resulting reaction cascades operate to transmit cellular signals [6]. In cell signaling, oscillation has been reported in many studies, with circadian rhythms being a well-known example [9]. Calcium ion signaling oscillation is another well-known phenomenon [10]. Mathematical models have been proposed to explain signaling kinetic behaviors based on a set of kinetic equations. Systems biology approaches have also been developed in recent years [1,2,11].
Systems biology can describe the kinetics of a signal pathway usingsimultaneousequations of a complex reaction network. On the other hand, there are other types of models consisting of many simultaneous reaction rate equations, including more than ten variables [12,13]; furthermore, systems biology models including more variables are also known [14]. Signaling networks frequently include non-linear reaction autocatalytic processes. To date, there have been many fascinating models of such autocatalytic reactions, enabling bifurcation and/or bi-stability in association with physical theory [5,15]. However, autocatalytic models or positive feedback are not necessarily applicable to all biological signaling pathways.
In the current study, to understand biological signaling pathways, we used the following three novel perspectives (A)-(C) based on a non-linear and non-equilibrium kinetic model, which included only two concentrations of the signaling molecule to describe the biological signaling pathway. (A) An equation of the continuity of the chemical concentration of c i (i = 1…., n) including chemical reaction items can be described using diffusion coefficients D i , kinetic coefficients k i , and concentrations of individual compounds c i as follows: In the above formula, the diffusion rate is hypothesized to obey Fick's law. In general, the diffusion items and chemical reaction items are thus separately described. On the other hand, because the biological signaling pathway network, including protein interactions, is limited by the slow diffusion rate of the signaling molecular proteins, kinetic coefficients generally depend on the diffusion coefficient. Therefore, diffusion items and reaction items cannot necessarily be separated without validation in the biological reaction. In the current model, we therefore described kinetic coefficients in a diffusion-coefficient-dependent manner.
(B) A feedback process due to non-linear self-catalytic reactions was not assumed in the current model, but instead, interactions between signaling molecules in their diffusion was assumed to give non-linearity to the model.
(C) A model system far from equilibrium due to a continuous supply of chemicals from the outside was hypothesized. The main issue is how minimal extracellular changes can be transformed into intracellular environmental changes. We aimed to evaluate the behavior of the model around the critical state by perturbation expansion using a minimal change of the supplied molecule concentration. By this mathematical evaluation, we aimed to illustrate the dynamic continuous oscillatory concentration change of signaling molecules from a static state.
In the current study, we constructed a novel model and aimed to evaluate the general intrinsic properties underlying cellular signaling based on signaling molecule interaction kinetics. Previous systems biology models have not necessarily focused on the diffusion process of proteins. Given the non-linearity during diffusion, we assumed kinetic instability of the signaling molecule interaction, and the sensitivity of the cell signaling in response to the environmental change was evaluated. The model system consists of several steps as follows: (i) the signaling molecule achieves an interaction active state by reversibly binding a cofactor that provides the signaling molecule with interaction activity; (ii) the signaling molecule has the ability to hydrolyze the cofactor; (iii) the signaling molecule interaction activity becomes lower when binding a hydrolyzed inactive cofactor compared to the signaling molecule binding an active cofactor; (iv) the signaling molecule has the ability to exchange the inactive cofactor with an active one; (v) active cofactors are supplied continuously from the outside. Thus, we set the interaction activity to be self-limiting, causing dynamic instability of the signaling molecule interaction. In the present model, we assumed that a signaling protein diffuses relatively slowly in the cytoplasm, and the whole signal transduction is a diffusion-limited reaction (assumption A.1). Following the protein interaction, one of the signaling molecules is phosphorylated or is bound to a cofactor such as GTP or ATP. In the kinetics of the protein interaction, fluctuation analysis was not fully performed in spite of the greater fluctuation in concentration relative to the solution of small molecules. Furthermore, we systematically analyzed the roles of the fluctuation in cellular signaling.

Protein interaction kinetics
The model scheme is shown in Figure 1. There are two types of signaling molecule, an active cofactor-binding signaling molecule (X), and an inactive cofactor-binding signaling molecule (Z). An active cofactor is non-hydrolyzed, and the inactive cofactor is the hydrolyzed type. X has the higher interactive activity and Z has the lower interactive activity. First, X can associate with oligomeric enzyme complex R consisting of X and Z, which transforms the active form X into the inactive form Z by hydrolysis of the binding cofactor and is to be released as Z irreversibly: ð2:1Þ In the above formula, we assumed that signaling proteins diffuse relatively slowly in the cytoplasm and that the dissociation rate (relating to the irreversible orientation k 21 ) of an encountering pair X and Z is significantly slower than the hydrolysis rate of the active cofactor changing into inactive cofactor (relating to k 2 ) (assumption 2; A.2). On a simple consideration of the diffusion limited step, when the kinetic rate can be described according to Fick's law using diffusion coefficients D X and D Z of X and Z, respectively: Next, Z recovers its interaction activity by exchange active cofactor P into inactive cofactor P', returning to X ( Figure 1).
Signaling molecules have the potential to hydrolyze the cofactor by interacting with identical species:

Kinetic equations of interaction signaling molecules
Here, kinetic equations were set according to the above simple reaction cascade. The equations consist of protein interactive items and an item of small-molecule-cofactor exchange. When the concentration of the protein is sufficiently small, the dependency of the diffusion coefficient on the concentration is linear [11,[13][14][15][16][17][18]. In comparison with the exchange kinetic rate of the cofactor (2.3), the rate of macromolecular protein interaction that depends on the diffusion step can be regarded as significantly smaller ((1.1), (2.4), & (2.5)) as a general (assumption 1; A.1). In this case, the whole reaction system can be regarded as diffusion-limited [4,19,20], and the diffusion rate is given according to Fick's law using a gradient of concentration. Thus, we get the protein interaction kinetics equations using the diffusion coefficient: Here, the exchange rate of the cofactor is expressed by the item k 3 ' PZ. Further, Because kinetic coefficients depend on diffusion coefficients, we set: Using (2.10), (2.8) and (2.9) are rewritten as follows: Here we set the oligomer concentration as constant because de novo asssembly is considered to be much slower than monomer interaction at the steady state (assumption 3; A.3): ð2:13Þ Here, a small letter affixed signifies values at the steady state. This assumption is based upon on the steady state in the protein assembly [21][22][23][24]. Setting the right hand sides of Eqs. (2.11) and (2.12) gives:

Fluctuation of diffusion coefficient
Subsequently, let us consider the fluctuation of participant proteins. We set: X~X e zx, Z~Z e zz ð3:1Þ In actuality, receptor R interacts with other proteins, X and Z, in the course of diffusion (Figure 1). In actual signaling pathways, signaling molecules associate with other signaling molecules and phosphorylate them or are phosphorylated by them. The diffusion coefficients can be altered in proportion to the signaling molecule concentration. By using the Gibbs-Duhem expression, the diffusion coefficient D of one macromolecule in the solution can be generally written as [13]: where T is the temperature of the solution, k B is the Boltzmann constant, and g is the frictional coefficient of a given macromolecule in solution. A is the second virial coefficient, v is the partial specific volume of protein with molecular weight M, and N A is Avogadro's number. The small letter c denotes the concentration of the solute. Further, we hypothesized that D R , D X , and D Z , the diffusion coefficients of R, X and Z, are given by extension of (3.2) to the mix solution of two macromolecules, X and Z: where v X and v Z are the partial specific volumes of X and Z with molecular weights M X and M Z , respectively. A X A X ', A X '', A Z , A Z ', and A Z '' are the second virial coefficients. In actuality, X and Z are the same molecules except with bound ATP or ADP. The fluctuation of the diffusion coefficient is given by: And therefore, we set: Here, an increase in X contributes to a decrease in diffusion coefficient D 1 in thefluctuation item, ax (a.0), because of the higher interaction activity that reducesdiffusion; in contrast, increased Z contributes to increases in the diffusion coefficient D 1 in the fluctuation item bz (b.0).

Calculus simulation of concentration oscillation
The time-course of the signaling molecule concentrations was simulated via the substitution of appropriate numerical values into (3.9). A numerical calculation was performed over a sufficiently long period to evaluate the trend of signaling protein behavior. In the current simulation, the fluctuation coefficients a, b, c, d, e, and f, of D 1 , D 4 , and D 5 are of the same order of magnitude, 10 2 (100,856). The concentrations of X and Z at the steady state are given by Eq. (2.14). On the basis of A.3, D 1 , the diffusion coefficient of the assembling rate of X to R, is significantly smaller than D 4 and D 5 , which are diffusion coefficients of the assembling rate between X and X, and X and Z. Applying the above conditions, the simulation results are shown in Figure 2. When p increases to the values that satisfy: the concentration of the signaling molecules continuously oscillates, showing nearly equivalent frequency and amplitude at oscillatory individual peaks, except for at the initial phase. Subsequently, we evaluated the mean amplitude and frequency of the fluctuation oscillation. From t = 50, around when the oscillation initiates, to t = 1000, the mean amplitude was calculated as the division of the sum of each size of the peak by the number of peaks. In addition, the frequency was estimated by dividing the number of the peak by Dt = 950.
Notably, the present simulation shows that both frequency and amplitude are nearly proportional to the logarithm of e = p2p c ( Figure 3). Namely, vf w! log e ð4:2Þ These formulae imply that the outside alteration is transformed inside into the information of cell signaling. On the basis of the above simulation results, the present signaling model system is characterized by: ( 5 -  Using the eigenvectors of L, (l 1 , l 2 ), coordinate transformation is performed as follows: The above parameters are subsequently set as: Further u is given by: (4) ð6:5Þ Therefore, using we can then obtain using (6.5) and (6.6): (4) ð6:7Þ Solving the above (6.7), the coefficients of a i in (6.5) are given (see Appendix S1: mathematica cord [#55-70]). By substitution of u that is given by n and e into f v (u,v) in (6.4), we can obtain the kinetic stability equation of fluctuation n using coefficients n i (i = 1, 2, 3, 4, 5, 6) as follows (See Appendix S1, Out [#73] in the mathematica cord): _ v v~n 1 v 2 zn 2 "vzn 3 " 2 vzn 4 v 2 zn 5 "v 2 zn 6 v 3 zO(4) ð6:8Þ To evaluate the amplitude of fluctuation, setting the right hand equal to zero, n 1 v 2 zn 2 v"zn 3 " 2 vzn 4 v 2 zn 5 "v 2 zn 6 v 3 zO(4)~0

ð6:9Þ
As a result, when e,7.0610 24 , Eq. (6.9) has two real number solutions of v other than zero, indicating the bi-stability of the fluctuation v; when e$7.0610 24 (6.9), it has only the zero solution of the fluctuation (See Appendix S1; mathematica cord Out [#74]). Therefore, a bifurcation of the fluctuation with respect to the value of e is predicted (Figure 4). Because v is simply given by the linear equation (6.3), this directly demonstrates the amplitude bifurcation of x and z.

Discussion
Here, we presented a model of cell signaling systems and performed mathematical analysis on the model in addition to numerical simulations. An increase in the supply of the cofactor near the critical concentration induces a 'phase transition' of the system, indicating that the model system has the ability to transform information on the concentration change of a cofactor outside the system into inside information, i.e., the amplitude or frequency of the concentration oscillation of the signaling molecule. This term, 'phase transition', is a metaphor implying that the model system nearly discontinuously acquires the ability to dynamically transform outside information into inside information.
The introduced non-linear kinetic equations include only two independent parameters, the active or inactive cofactor binding protein. The observed oscillation of signaling molecules in the simulation is not a chaotic behavior that requires more than two parameters [25]. However, the fluctuation of the signaling molecule concentration shows chaos-like oscillatory behavior. In fact, neither the amplitude nor frequency of every oscillation is precisely constant for a lengthy period as shown in the trajectory.
We will report mathematical validation of the result elsewhere. We will report mathematical validation of the result elsewhere.
The simulation allowed us to define the formula (5.3) mentioned above. The formula using the logarithmic function brings to light an important idea. The present simple model illustrates an essential property in which a system is relatively stable to the outside environment, because a minimal increase in the concentration of the cofactor is transformed into the logarithm of the concentration change inside the cell system. These features may be crucial in responding to transformations of the outside environment while minimizing the inside influence caused by outside alteration. Further, the present simple formulation is reminiscent of the definition of entropy in informatics. In conclusion, our model indicates that the non-linearity of a protein interaction theoretically gives an interesting cohesive behavior, such as an oscillation, of the signaling system leading to self-organization in vivo. Still, the theoretical basis of cell signaling systems for quantitative evaluation requires further formulation on the basis of experimental study.