Nonlinear signal transduction network with multistate*

Han-Yu Jiang(姜寒玉) and Jun He(何军)

School of Physics and Technology,Nanjing Normal University,Nanjing 210097,China

Keywords: signal transduction,biological network with multistate,CICR,nonlinear signaling

1. Introduction

The signal transduction is important to a variety of intracellular and intercellular processes.[1]The signal propagates between different cellular components and communicates between biological cells. It usually mediates between two receptors by signaling molecule,which changes the conformational state of receptor protein. Generally,the response of receptor is typically a nonlinear function of the local concentration of signaling molecule.[2-4]In this way,signal transduction systems can perform robust switchlike operations.[2,5]The state transition of the receptors triggers or inhibits the release of signaling molecule, which then affects other receptors’ state. With such mechanism, the signal propagates on the network of intracellular and intercellular receptors. Given that signal transduction involves a large number of receptors with complex connections by signaling molecule, it is natural to apply the network science to understand their dynamics.[6-9]The network sciences has been widely practiced in different research fields.[10-14]To keep and regulate the life activity, there exist a variety of signaling networks, such as calcium signaling network,protein-protein interaction networks,and excitationcontraction coupling in muscle tissue.[15-17]Up to now,much efforts have been paid in this direction and some key insights have been achieved using such interdisciplinary approach.[18]

In the signal transduction network,the receptors are taken as the nodes which are connected by signaling molecule.Generally speaking, the receptors are not uniformly distributed spatially, which makes the network architecture more complex. It is also well known that the intracellular space is crowded with organelles and various obstructions.[19]The diffusing signaling molecule has to navigate between complex intracellular structures. Besides, the receptors usually embed on the biomembrane,which is not flat. For example,the ryanodine receptor (RyR) distributes on the endoplasmic reticulum,which is folded in the cytoplasm,and some parts of different sheets may be very close to each other.The calcium signale transduction also happens at the ER-PM junction.[20-22]In this way,signals can be communicated between distant regions in the biomembrane.The connections between receptors can be very complex. In the words of the network science,the signal transduction should occure on some types of the networks. However, based on the current knowledge, we do not know their explicit network architectures. Thus,it is interesting to study the effect of different network architectures on signal transduction,especially whether the biological function performs stably and similarly with different network architectures.

The bistable regime can emerge from the network dynamics in a study of the propagation of nonlinear signaling on the network[4]as in electronic systems[3]where the network architecture was not considered. In their work, only two states, a closed state and an open state, were considered for a receptor.[4]In many signaling system, the receptor has more conformational states,the transition between these states may exhibit more complex patten of the signaling transduction. The Ca2+-induced Ca2+release (CICR) mechanism is a basic mechanism of the calcium signal transduction system, which is important second messenger in the biological system, and widely involves in the regulation of cell life activities, including cell membrane permeability, cell secretion,metabolism and differentiation.[23,24]In the literature,[25]a simplified mechanism that mimics “adaptation” of the RyR has been developed to reproduce experimental data from cardiac cells. In such model, the nodes of the calcium signaling network, RyRs, have four states. Moreover, among the transitions between these states, there are two fast processes and one slow process. In this study, we will study the CICR on the network to explore the behavior of signaling networks with multistate in which nodes are regulated by reaction rates nonlinearly. In the next section,we will present explicit CICR mechanism and corresponding mathematical description of the states and the networks considered. Then, the dynamics of CICR on three networks will be studied, and time revolution of fractions of states will be simulated with different network architectures. The pattern of steady states will be also studied to provide more informations about the possible bistable regime from the CICR mechanism.

2. CICR mechanism

In the current work, we take an RyR,the receptor of the CICR mechanism,as a node of network. RyR has four states as shown in Fig. 1. StatesC1andC2are closed states, in which no calcium will be released, and statesO1andO2are open states, in which RyR releases calcium. The transitions between different states of an RyR in an environment with concentration [Ca2+] are also illustrated in Fig. 1. The relevant parameters have been extracted from the experiment,and listed in Table 1.

Fig. 1. Schematic diagram of transitions among the four states of the RyR.[25]

Such model was proposed in Ref.[25]to mimic the experiment.The statesC1andC2and their conversions are the basic mechanism, and such two-state model’s behavior on network was also discussed in Ref. [4]. The stateC1is dominant at a low concentration [Ca2+], for example 0.1 µM. If the [Ca2+]increases,the RyR will be activated from the closed stateC1to an open stateO1at a rate ofk+a[Ca2+]n, and deactivates back toC1state at a rate ofk−a. Such simple two-state model is not enough to describe the experimental phenomenon. As suggested in Ref. [25], to keep the plateau open probability, the open stateO1also may be activated to the second open stateO2at a rate ofk+b[Ca]mand back to the first open stateO1at a rate ofk−b. It is also important to obtain the bistable regime in the current work. To describe the adaption phenomenon,the transition between first open stateO1and the second closed stateC2is also added, which is independent on the concentration[Ca2+]. And compared with the rates of processesC1⇆O1andO1⇆O2, the processO1⇆C2should be slow, which makes the decrease of the fractions of open states be a longtime-scale process. With such model, the main experimental behavior can be reproduced. For simplification, the conversions betweenC1,C2andO2are neglected as in Ref.[25].

Table 1. RyR kinetic constants.[25] The parameters n and m were also determined as 4 and 3,respectively.

With the above mechanism,we study the dynamics of the CICR on the network. For a network withNnodes, the state of a node can be denoted as vectors,[1,0,0,0],[0,1,0,0],[0,0,1, 0],[0,0, 0, 1]forC1,O1,C2,andO2,respectively. The overall state of the network can be written as a 4×NmatrixS(t),which revolves with timet. Besides the nodes,i.e.,RyR,we need connections between them to study the dynamics of the CICR network. As in Ref.[4], we assume the concentration[Ca2+]at an RyR is determined by a general background concentrationc0and calcium released by other RyRs. Hence,it can be obtained by the states of other nodes as

wherem=1,2,3,4 corresponds to statesC1,O1,C2,andO2,respectively, andiorjare the number of the nodes. It is nature to assume a small general background concentrationc0,0.1µM,in this work,and the concentration is mainly from the calcium released from other nodes. Hence, only the nodes in two open states contribute to the concentration, which is denoted asS2jandS4j. TheAijdenotes adjacency matrix,which is different for different types of networks. The parameterrgives the strength of connection between two nodes,which is assumed as a general constant. A more appropriate parameter can be defined ass=r〈k〉/c0with the average degree〈k〉=1/N∑Ni kiandki=∑Ni Aij,which will be adopted in the following calculation (its meaning can be seen more clearly under the mean-field ansatz in Section 4).

3. Dynamics of CICR on the network

In this section,we consider the dynamics of CICR on the network. Since it is difficult to detect the explicit network architecture of the calcium signal transduction, in our study we perform the investigation with three different types of networks,Erd¨os-R´enyi(ER)network,Watts-Strogatz(WS)network, and BaraB´asi-Albert (BA) network. The networks are generated with the help of the NetworkX package in Python Language. In the ER network, the two nodes are connected with a fixed probabilityp. In the WS network, every node connects tokneighbor nodes, and two nodes are connected with fixed probabilityp. Such network may be more appropriate to describe the behavior of the CICR mechanism because usually the calcium released by an RyR induces the release of neighbor RyRs on the endoplasmic reticulum while some other distant RyRs may be also activated due to the folding of the membrane. We also consider the BA network, which is formed by adding a sequence ofmnew nodes to an existing network. By the generated network,we have the network adjacency matrixAin Eq.(1).

With overall state matrixS, it is easy to characterize the activity of the network by computing the fraction of node which is in an stateiat timetas

To study the dynamics, we also should consider an ensemble of initial conditions where activity level of network is varied. We choose initial conditions that the nodeihas a probabilityhatC1and 1−hatO1. For simplicity,we do not consider other possible initial conditions withO2andC2states,which should provide similar conclusion. Here we simulate initial conditions whereh=i/(K −1)andi=0,1,...,K −1.With such treatment, we can explore the time revolution of a range of random initial conditionsSi(0).

Fig.2. The fractions of states C1 (blue dot),O1 (magenta dot),C2 (green dot),and O2 (dark dot)on three networks as a function of time for an ensemble of initial conditions with strength of connection s=1. The panels(a1),(a2),(a3)are for the ER networks with p=0.1 at short,medium,and long time scales,respectively. The panels(b1),(b2),(b3)are for the WS networks with K=5 and p=0.1,and the panels(c1),(c2),(c3)for the BA networks with m=10. In the calculation,105 time steps are adopted.

With the network adjacency matrixAand the initial state matrixSi(0)obtained as above,we can simulate the time revolution of the dynamics of the CICR on the network under the mechanism in Fig.1.In the simulation,we adopt network withN=500 nodes, and chooseK=6 random initial conditions.We consider the time revolutions of the fractions of four states on three networks considered. The parameters are chosen asp=0.1 for ER network,k=5 andp=0.1 for WS network,andm=10 for BA network. Since we need synchronize the time for all nodes in the network to give the explicit time revolution of fractions of four states of node, the Gillespie algorithm is not practical,so we perform this work by fixing time step to simulate the transitions between the four states of the RyR in the current work.

In Ref. [4], the nonlinear network with a closed and an open state exhibits bistable regime with a smallηin their work. If we only considerC1andO1state, an explicit calculation suggests that there are no bistable regime emerging because the values of ratesk+aandk−ashown in Table 1,which were extracted from experimental data leads to a largeη. It is interesting to see that the bistable regime recovers after stateO2is added. However, the inclusion of the slow processO1⇆C2breaks such bistable regime,which can be seen in the results at different strengths of connectionsillustrated in Figs.2-5.

In Fig. 2, the results with small strength of connections=1 are presented. We provide the results at both short time and long time scales to show the effects of fast processes and slow process. The results at medium time scale are also given as a reference of general picture of the time revolution. The time revolution of the fractions of four states on three networks exhibit very similar behaviors. The open stateO1deactivates within a relaxation timeτ～0.1 s(see Figs.2(a1),2(b1),and 2(c1)). At such time scale,the stateO2is nearly not activated,and the fraction of stateC2increases very slowly. Such pattern keeps up to a time about 1 s (see Figs. 2(a2), 2(b2), and 2(c2)). At the scale of 10 s, the slow process takes effect.The fractions of statesC2andC1increase and decrease a little,respectively,and approach to steady states(see Figs.2(a3),2(b3),and 2(c3)).

Fig.3. The fractions of four states. The same as Fig.2,but with the strength of connection s=5.

In Fig.3,the results withs=5 are presented. Compared with the results with strengths=1, the behaviors of the activity of the network become more complex. At short time scale,the stateO1deactivates more rapidly than the case withs=1. The fast deactivation is due to the largers. The larger connection between the nodes makes the propagation of the calcium signal faster,which consequently leads to fast deactivation. However,different from the case withs=1,at a time about 0.02 s,the deactivation ofO1becomes relatively slower,which means a steady state,but disturbed by the slow processO1⇆C2. With the time revolution, the fraction of stateO1approaches to zero finally as shown in the mediate panels. In fact,if we neglect the slow process,the activity of the network enters a bistable regime at long time scale. However, the effect of the slow processO1⇆C2makes the time not enough to exhibit such behavior,which will be discussed later.At long time scale,the tendencies of fractions of statesO1,C1,andO2at short time scale are smeared. A stable pattern as in the case withs=1 emerges at timet>20 s. Besides, with strengths=5,three networks also exhibit analogous behaviors.

Fig.4. The fractions of four states. The same as Fig.2,but with strength of connection s=10.

In Fig. 4, the results withs=10 are presented. In this case,the quasi-bistable regime becomes more clearly.At short time scale,the quasi-steady states can be seen clearly(here and hereafter, we will call the approximate steady state and approximate bistable regime at short time scale as“quasi-steady state” and “quasi-bistable regime” to distinguish them from the real steady state at long time scale and bistable regime only with fast processes). With the increase of the strength of connections, the calcium signal becomes easier to propagate from a node to another one. The relaxation time of the fast process becomes smaller, about 0.01 s with connections=5. However, as shown in Fig. 1, the slow process is independent of thes, which is the reason why the quasi-steady states here are more obvious than those with small strength.One can observe the transitions from a quasi-steady state to another quasi-steady state. For example,on WS network,the stateO1has two quasi-steady states, one is about 0 and the other about 0.2. At a time about 0.05 s,the network leaves the former and jumps to the latter. The stateO2also may jump from a sate about 0 to the state about 0.8. The values of the quasi-steady states for ER and WS networks are almost the same. However,the BA network has a larger fraction for stateO2. At a time scale of 0.1 s,the fraction of stateO2decreases until reaching zero.

Fig.5. The fractions of four states. The same as Fig.2,but with strength of connection s=45.

The decrease of fraction of stateO2can be seen at time scale of 1 s. The fraction of stateC2increases first,but at long time scale its behaviors are different for different networks.On ER network,the fraction of stateC2tends toward a stable value about 0.1 and the one of stateC1toward a value about 0.9. The nodes of the network are almost occupied by these two states at long time scale,the fractions of statesO1andO2are about zero. On WS network, the stable value for stateC2is about 0.5, higher than the one for stateC1, about 0.4. Besides,the system approaches steady states at time smaller than 10 s, which is different from the pattern on the ER network.For BA network,the tendencies of fractions of four states are analogous to the case on ER network. However,the statistical transitions between the two steady states and stateO1by fast processes make the situation more complicated because statesC1andC2connected byO1,which is usually very low at long time scale.

Now,we turn to a large strength of connections=45,and the results are presented in Fig.5. With such large connection,the calcium signal propagates on the network very fast.Hence,at short time scale, one can find that after a very small relaxation time, the states of the network concentrates to the stateO2. The quasi-bistable regime in previous figures disappears in the large strength.With time revolution,the fraction of stateO2decreases,and reaches a value about 0.2 or smaller at larger time scale. The stateC2becomes dominant. The fraction of stateO1keeps a low level at all time scales,while the stateC1is observable on WS scale at long time scale.

4. Pattern of steady states on network

In the above discussion,we present the time revolution of the fractions of four states on three networks. One can find that at small strength of connections, the network is nearly not activated. With increasing of the strength of connection,quasi-bistable regime emerges at short time scale, where the fast processes play the dominant role. However, the low process breaks the quasi-steady states with the time revolution,and network revolves to the steady states at long time scale.And with a very large strength of the connection, the quasibistable regime disappears. Since the steady states are more interesting in biological function. In this section, we discuss the dependence of the quasi-steady states on networks.

To understand the basic feature of the network dynamics,we first provide the mean field(MF)ansatz of the system considering homogeneous.[4]Under this ansatz,we assume all nodes revolve with time in the same manner, which leads to the state functionSi j=Si0withiandjdenoting the four states andNnodes, respectively. Hence the calcium concentration can be written as

withpC1+pO1+pC2+pO2=1.

Fig. 6. The steady states on networks with different parameters with the variation of the strength of connection s. The lines are for the MF ansatz with the process O1 ↔C2 turning off(Eq.8). The symbols“·”and“+”are for the results at small time scale with full model,and results after turning off the O1 ↔C2 process,respectively.

Based on the parameters in Table 1, the processes in Equations from Eqs. (4)-(7) are fast processes while the last one is slow process.If we consider all three processes,we will reach the steady states at long time scale. Hence,first,we consider the case with the last processO1↔C2turning off. With such treatment we can obtain steady states at a long time scale which correspond to results of full model at short time scale.The steady states corresponds to the stationary points, which satisfy the algebraic condition as

wherepC1+pO1+pO2=1,and we still use a low backgroundc0=0.1 µM and consider the variation of the parameters,which reflects the strength of the connections between nodes on the network.

In Fig.6,we present the steady states obtained with above conditions. Under the MF ansatz, the difference of networks is smeared, and the effect of network is absorbed into the strength parameters. For the results with low strength, the network can not be activated as in Fig. 2, which leads to a meaningless steady state withC1=1. With the strengthening of the connection between nodes,a bistable regime of fraction of stateC1emerges ats～5.Besides the steady state atC1～0,there exists another steady state at about 0.2, which decrease rapidly with the increase of the strengths. The steady state at about 1 for sateC1keeps until ansabout 20, and disappears there. For stateO1,there exists a steady state at smalls,and a new steady state about 0.5 appears atsabout 5 and decreases to zero slower than stateC1. At ansabout 20,the steady state about 0 for stateO1disappears and only one steady state keeps.The stateO2also exhibits bistable regime with a steady state about 5 and a steady state about 0. The higher one increases to about 1 with the increase ofs, and the steady state about 0 disappears atsabout 20 also. Hence, under MF ansatz, for all three states (theC2state is removed here), the fractions of states are monostable ats5, bistable at 20s5, and monostable again ins20.

Fig. 7. The steady states on networks with different parameters with the variation of the strength of connection s. The lines are for the MF ansatz of full model(Eq.(9))and the symbols“·”is for the results of the simulation of full model at long time scale.

For the realistic networks,we still consider the three networks in previous calculation with more parameters. With theO1↔C2process turning off,the steady states can be reached at time long enough. As shown in Fig.6,one can find that the model only with fast processes produce a result(symbol“+”),which fit the curves under MF ansatz very well on the ER and WS networks. The result on BA network deviates from the MF results a little larger. Besides, the results with different parameters of each network are similar to each other. We also provide the results with full model(symbol“·”). Because the pattern of the steady states will be broken by the slow process,we choose time pointst=0.4×0.1sto extract the values of fractions where the quasi-steady states can be seen obviously but the effect of the slow process is still very small. It can be found that the quasi-steady states are close to these only with the fast processes at most values ofs, except about 6.As shown in Fig. 3, with a strength about 6, the effect of the slow process becomes obvious within the relaxation time of the quasi-steady states. Generally speaking,the simulation results are close to these under MF ansatz.

Now, we consider the cases with slow process, the stationary points should satisfy the algebraic conditions as

withpC1+pO1+pC2+pO2=1.

In Fig. 7, we present the steady states obtained from the above conditions. At first sight,one can find that after including the slow process,the pattern of steady states changes completely. It is also the reason of the different pattern for time revolution at short and long time scales shown in Figs. 2-5.For thesfrom 0 to about 50, no bistable regime emerges for all four states on three networks. The MF results suggests that with weak connection,smalls,the fraction of stateC1is about 1 and decreases rapidly at about 30,and to a very small value at about 50. The fraction of stateC2exhibits almost reverse behavior, about 0.1 atsabout 0, and increases fast atsabout 30,and up to about 0.8 atsabout 50. The network is dominant by these two closed state,which satisfy the requirement of the biological function.

The simulation results of the full model are also presented in the Fig.7. Different from results at short time scale,the results for different networks and parameters are different, and deviate from MF results obviously. For ER network,the cross point ofC1andC2is lower than MF curves. With the increase ofp, the simulation results become closer to the MF results.Such situations can be found also for WS and BA networks.

5. Discussion and summary

Signal transduction is an important mechanism in living things to regulate cell life activities. The signal usually mediates from a receptor to another receptor by releasing signaling molecular which can change the states of the receptor. It suggests that the signal propagates on a biological network with receptors as nodes. Hence,to study the mechanism of the signal transduction,we need consider the network architecture.

The signal transduction on the network may exhibit different behaviors from these without considering the network architecture. In Ref. [4], the authors studied the dynamics of a two-state transition model on networks. The bistable regime was observed in their calculation. However,if we apply their model to the important CICR mechanism in the calcium signaling transduction, the two-state transition model, statesC1andC2in the current work, does not produce any bistable regime in our calculation. It is interesting to see that the bistable regime is reproduced if we include second open stateO2to statesC1andO1. Such results suggest that the network with more states provides more biological patterns than in the two-states model.

In the CICR mechanism, the transitions between statesC1,O1,andO2are all fast processes. The inclusion of stateC2introduces a slow process,with which the network deactivates with time revolution. The dynamics of the CICR exhibits different patterns at different time scales. Though with this slow process the system tends to closed states in any case,with different strengths of the connection the system tends to different closed states. Weak connection between nodes leads to a large fraction ofC1state while strong connection to a large fraction ofC2state.

In this study we consider three types of networks,Erd¨os-R´enyi network,Watts-Strogatz network,and BaraB´asi-Albert network, with different parameters. The pattern at short time scale is not sensitive to the network architecture,and close to MF results. At long time scale,only monostable regime is observed. The difference of network architectures affects the results more seriously,and deviates from the MF results at most cases.

In summary,our study shows that nonlinear network with multistate exhibits rich patterns,especially after the processes with different time scales are included. The second open state of CICR is found essential to reproduce the bistable regime at short time scale. Our finding identifies features of biological signaling networks that may underlie their biological function. The results are also helpful to understand other polymorphic networks. Based on the model in the current work,more complex phenomenon of calcium signal, such as intracellular calcium concentration oscillates,can be investigated in future study after including more factors, such as leak and the SERCA pump for a closed cell, and additional PMCA pump and Ca2+influx for open cell. With abundant experimental information, the mechanism of the calcium signal transduction on the network can be further understood.Vise versa,the further explicit analysis of experimental data may be helpful to unveil explicit network architecture of the calcium signal transduction.