Oscillations Governed by the Incoherent Dynamics in Necroptotic Signaling

Emerging evidences have suggested that oscillation is important for the induction of cell death. However, whether and how oscillation behavior is involved and required for necroptosis remain elusive. To address this question, a minimal necroptotic circuit is proposed based on the CNS pathway. Stochastic parameter analysis demonstrates that the essential structure for oscillation of the CNS circuit is constituted by a paradoxical component embedded with positive feedback among the three protein nodes, i.e., RIP1, caspase-8, and RIP3. Distribution characteristics of all parameters in the CNS circuit with stable oscillation are investigated as well, and a unidirectional bias with fast and slow dynamics that are required for high occurrence probability of oscillation is identified. Four types of oscillation behaviors are classified and their robustness is further explored, implying that the fast oscillation behavior is more robust than the slow behavior. In addition, bifurcation analysis and landscape approach are employed to study stochastic dynamics and global stability of the circuit oscillations, revealing the possible switching strategies among different behaviors. Taken together, our study provides a natural and physical bases for understanding the occurrence of oscillations in the necroptotic network, advancing our knowledge of oscillations in regulating the various cell death signaling.


INTRODUCTION
Oscillatory signaling is ubiquitous in vital movement, and precisely controls the sequential behaviors of many important physiological processes [1][2][3], such as circadian rhythm [4], embryonic development [5], neuron electrophysiology [6], and cell cycle [7]. Oscillation behavior of molecules, for instance, calcium ions, tumor suppressor gene p53, p38 MAPK (mitogenactivated protein kinases), and NF-κB (Nuclear Factor-kappa B) in inflammatory response, determines the fate of cells [8][9][10]. Organisms can identify, encode, and transmit different biological information, and perform different responses with different modulation methods through oscillatory signaling, including amplitude modulation or frequency modulation [11,12]. The biological systems are nonlinear and quite complicated, and the dynamic behavior of oscillation in individual cells might be annihilated in population average dynamic behavior [9]. Thus, to accurately understand how biological oscillation controls biological functions is still a challenge.
Oscillations seem to be an important prerequisite for the processing of the cell death signaling [13][14][15][16][17][18]. The bistable, excitability, and oscillation are combined with spatial coupling mechanism like diffusion, which can propagate as trigger waves [19]. Recently, the spread pattern of self-regenerating trigger waves had been confirmed in the process of apoptosis [20]. Besides, reports have shown that the ferroptosis signal spreads stably in the cell population in a wave-like pattern [21,22]. Necroptosis, a passive and irreversible inflammatory form of cell death, is quite different from apoptosis and ferroptosis. Its physiological manifestations are cell volume expansion, cell membrane rupture, intracellular material overflow, incomplete DNA degradation, and local severe inflammatory reaction. Necroptosis is involved in many diseases, such as acute pancreatitis ischemia, cardio-cerebrovascular disease, and so on [23,24]. Whether and how the cellular necroptosis signaling (CNS) could spread and induce cell death in the form of oscillation-induced trigger waves remains unclear.
Mathematical modeling with qualitative or quantitative analysis has become an important approach for dissecting the complex regulatory mechanisms of biological system [25]. Numerous studies have explored the design principle of signaling circuits with various biological functions. The relationships between network motifs [26] (such as paradoxical components [27][28][29], feed-forward loop [30], and feedback loop [31][32][33]) and biochemical oscillator [34][35][36], adaptation [37], robustness [38][39][40][41], and noise attenuation [42][43][44] have been revealed successively. Designing and developing a signaling circuit system with regulatory function is of guiding significance for the further understanding of how oscillation signals regulate biological processes [13,45]. For instance, the minimum free-energy cost of biological oscillation to overcome the large fluctuations from the environment and maintain coherence has been successively revealed [46][47][48]. The design principle of a biological oscillation network with high accuracy and sensitivity under limited energy cost was further explained [49,50].
In our study, we propose a CNS circuit to systematically analyze the robust oscillation dynamics in the CNS pathway. The essential structure of the circuit for generating the occurrence of necroptotic oscillation is determined. Moreover, how each component in the network mediates the probability of oscillation occurrence is investigated. The two different incoherent dynamics, fast and slow dynamics within a unidirectional bias, are demonstrated to be the key mechanism for the high probability of oscillation occurrence. Ultimately, both the bifurcation analysis and potential landscape theory are performed to explore the switching strategies among different types of oscillation through regulating the interaction within the unidirectional bias.

Overview of the Core Circuit in CNS Pathway
The diagram of TNF-induced necroptotic signaling is shown in Figure 1A. Upon TNF stimulation, TNF binds TNFR1 on cell membrane and promotes oligomerization of TNFR1 to form homotrimer. Oligomerization of TNFR1 induces the conformational changes in its intracellular domain, leading to the recruitment of TRADD (TNFR1 associated death domain protein) and RIP1 (Receptor-interacting protein kinase 1) to form Complex-Ⅰ [51]. TRADD and RIP1 compete for the binding sites on TNFR1 in complex I to mediate the downstream signaling [52]. Then, the modulated RIP1 dissociates from Complex Ⅰ and combines with the RIP3 through RHIM-domain to form a "Necrosome" complex [53]. RIP3 is phosphorylated by RIP1 in necrosome [54,55]. Besides being activated by RIP1, RIP3 can also be autophosphorylated in the kinase domain of Ser227 (human)/Thr231-Ser232 (mouse) [54,56,57]. RIP1 also binds FADD through the death domain, and further recruit caspase-8 into necrosome [58][59][60]. TRADD can also drop from complex I to recruit FADD and caspase-8, which is an important process for the activation of caspase-8 [61]. Under normal circumstances, cleavage of RIP1 and RIP3 by caspase-8 limit the occurrence of necroptosis [62,63]. Besides, a very recently study reported that phosphorylated RIP3 suppresses caspase-8 activity through the recruitment of RSK [64]. Necrosome subsequently recruits and phosphorylates the downstream protein MLKL [65], which will form oligomers and transfer from cytoplasm to the cell membrane, resulting in the induction of necroptosis.
To intuitively describe the interaction of the core module, coarse-grained method [66] is employed to simplify the CNS pathway into circuit, and the schematic diagram is shown in Figure 1B. The circuit includes four components, i.e., TRADD, RIP1, caspase-8, and RIP3. The necroptotic dynamics are mainly determined by these four components in necrosome, while MLKL is not considered in the CNS circuit. This is because MLKL mediates cell necroptosis as a downstream substrate of RIP3. Studies have shown that knocking down MLKL expression by RNA interference will not affect the RIP3 or the other upstream signaling [67,68]. In the circuit, there are various interaction among the four proteins ( Figure 1C), such as negative feedback components, positive feedback components, and the paradoxical components [27][28][29], which can produce abundant dynamic behaviors [36,[69][70][71].

Searching the Essential Structure for Oscillation Within the CNS Circuit
To discuss the oscillation dynamics in the CNS pathway, we construct a set of self-evolving ODEs model based on the circuit shown in Figure 1B. The model has four variables and ten interaction terms. ODEs and the detailed description of the parameters in the interaction terms can be found in Supplementary Information S1.1. In the model, the components are switched between active and inactive forms by phosphorylation, dephosphorylation, or cleavage. For systematically analysis, Latin Hypercube Sampling method [72] is employed to randomly scan the network parameters in the widely parameter space for evaluating the stable oscillation behavior of the CNS circuit (Supplementary Information S1.2). For sampling uniformity, the exponential sampling range of k_i is −1 to 1, and the exponential sampling range of j_i is −3 to 2. The Hill coefficient n_i is a random integer sampling with 1-4. All the parameters k_i, j_i, and the Hill coefficient n_i of the model are sampled within the proper ranges based the previous studies [31,37]. We randomly select more than 30 million sets of parameter combinations, in other words, we analyze 30 million dynamic systems to test whether the system can produce robust oscillation. The probability of robust oscillation of pRIP3 (phosphorylated RIP3) is taken as the evaluation object and the analysis procedure is shown in Figure 2A.
We firstly intend to identify the key protein nodes and essential interaction in the CNS circuit ( Figure 2B) that required for oscillation generation. In the numerical simulation, 500,000 groups of random samples are taken for each simulation test, and each simulation test is repeated 5 times. For the four protein nodes, i.e., TRADD, RIP1, caspase-8 and RIP3, we respectively set their total amounts to 0 one by one and count the probability of oscillation behavior of pRIP3. This total amount is determined by the natural generation rate and degradation rate of protein, which is considered as unit 1 in the normalized model [9]. Thus, the steady states of the system could be regulated by the total amount of the proteins. As the result shown in the upper panel of Figure 2C, TRADD is not required and the probability of circuit oscillation without considering TRADD exhibits nearly 3-fold compared with the full model (FM). While RIP1, caspase-8, and RIP3 are essential for the circuit to oscillate. We next explore the key interaction for the generation of oscillation by limiting the corresponding parameters. There are 10 interaction terms among TRADD, RIP1, caspase-8, and RIP3 in the circuit, which are mainly described by the 10 parameters (k 1 -k 10 ). The 10 parameters in these terms are in turn fixed at 0, respectively. As the result shown in the lower panel of Figure 2C, when the four parameters of k 6 , k 7 , k 8 , and k 10 are limited, the probability of pRIP3 oscillation is zero, indicating that the interaction terms characterized by these four parameters are necessary for producing oscillation behavior. Moreover, two interaction parameters are simultaneously fixed to zero in turn in our further analysis (Supplementary Figure S2), supporting the results shown in Figure 2C. The statistical results of simultaneously fixing three or more parameters to zero are similar. Therefore, based on these statistical analysis, the essential structure for generating oscillation of the CNS circuit is determined and highlighted in Figure 2D. The structure is constituted by a paradoxical component embedded with a positive feedback among three protein nodes (i.e., RIP1, caspase-8, and RIP3). Despite the parameters of k 3 and k 4 are not necessary for the oscillation behavior, the probability of the system is only one-third of FM when k 3 or k 4 is set to 0 ( Figure 2C, down panel), indicating these two interaction terms are also important for the oscillation induction.

Oscillation Determined by the Fast and Slow Dynamics Within a Unidirectional Bias
Having identified the essential structure for generating robust oscillation, we further explore the control mechanism of how the behavior is generated and regulated by the key protein nodes and essential interaction. Figure 3A presents how the network oscillation probability is regulated by the total amounts of the three key nodes, i.e., RIP1, caspase-8 and RIP3. As the results suggested, the oscillation probability increases with the increase of the total amount of RIP1 or caspase-8. A 10-fold increase of their total amounts will improve the oscillation probability by nearly 7-fold. However, variation of RIP3 amount has no significant effects on the probability. Therefore, high amount of RIP1 or caspase-8 will accelerate the occurrence of oscillation in the CNS circuit.
We next select more widely the parameter combinations that can generate oscillation. By collecting 33 million sets of parameters, we obtain 14,906 sets of oscillation parameter combinations that meet the appraisal standard (Supplementary Information S1.1). We cluster the 14,906 sets of parameters, however, no matter what clustering algorithm (such as k-means clustering, spectral clustering, etc.) is adopted, it is unable to distinguish their parameter features. The result of spectral clustering is shown in Supplementary Figure S3A, The red mark highlights the essential structure of the CNS circuit for generating oscillation.
Frontiers in Physics | www.frontiersin.org August 2021 | Volume 9 | Article 726638 meaning that the parameter set that produces oscillation behavior of the system only follows one distribution characteristic. Thus, we could project the high-dimensional parameter space to each dimension to investigate the distribution characteristic of a single parameter. The distribution characteristics of the four identified essential interaction terms ( Figure 2D) in the parameter space are investigated and the corresponding results are presented in Figure 3B, indicating that k 6 or k 10 should be large, k 7 is at an intermediate value, and k 8 trends to be small to promote the occurrence of oscillation.
If an interaction term is constrained to be beneficial to improve the oscillation probability of the system, then the term is functionally significant [35]. We constrain the interaction terms in most probable The probability of RIP3 oscillation is substantially improved when two interaction terms are constrained during the random parameter search. The "+" indicates that the term is constrained, and the "−" indicates that the term is unconstrained. (D) The heat maps that reflect the variation of probability when two-parameters are changed. (E) Schematic diagram of the fast and slow dynamics within the essential structure. Thick red lines indicate fast dynamics and thin blue lines indicate slow dynamics.
Frontiers in Physics | www.frontiersin.org August 2021 | Volume 9 | Article 726638 5 strength range according to the statistical results in Figure 3B. In numerical simulation, we first constrain a single interaction term, and then constrain two groups of interaction terms simultaneously ( Figure 3C). The result shows that when the interaction term of RIP3 inhibition by caspase-8 (k 10 ) is restrained alone, the oscillation probability is improved most significantly, and the probability is increased about 7.3-fold compared with all parameters randomly sampled (yellow bar in Figure 3C). While restricting the two interaction terms characterized by k 10 and k 7 (inhibition of caspase-8 by RIP3) at the same time significantly increases the oscillation probability 19.5-fold compared with all parameters randomly sampled.
We further systematically analyze the relationship between oscillation probability and two-parameters variations. As the four heat maps shown in Figure 3D, the high oscillation probability is always in the top left of the diagonal line, indicating the larger value of k 6 or k 10 compared to k 7 or k 8 . Thus, the most probable dynamics of the essential structure are shown in Figure 3E, indicating that RIP1 activates caspase-8 with a fast dynamic (represented by k 6 ), while RIP1 phosphorylates RIP3 with a relatively slow dynamic (represented by k 8 ). Moreover, caspase-8 blocks RIP3 (represented by k 10 ) with a relatively faster dynamics compared with the suppression of RIP3 by caspase-8 (represented by k 7 ). Therefore, above analysis indicates that the fast and slow dynamics within the essential structure are functionally significant for the occurrence of oscillation.
Aside from the essential structure, how the rest interaction terms within the CNS circuit mediate the oscillation probability are studied as well to identify the most probable structure for robust oscillation. In addition to the distribution characteristics of four parameters shown in Figure 3B, the distribution characteristics of other parameters are shown in Figure 4A and Supplementary Figure S3B. As the results indicated, the strength of this interaction term of RIP1 activated by RIP3 (represented by k 3 ) trends to be strong to promote the occurrence of oscillation, while the other terms prefer a weak or at an intermediate strength. Hence, the most probable structure for oscillation of the CNS circuit is proposed and presented in Figure 4B. Apparently within the structure, a unidirectional bias is captured by the two feedback loops formed by RIP1, caspase-8, and RIP3 ( Figure 4C). That is, the negative feedback loop interaction formed by k 3 , k 6 , and k 10 , and the positive feedback loop interaction formed by k 8 , k 4 and k 7 . The negative feedback loop presents fast dynamics, while the positive feedback loop presents slow dynamics, which provides a general underlying mechanism for the oscillation behavior of the CNS circuit.

Classification of Oscillation Behaviors by Amplitude and Period
Our above analysis determines the essential structure and further reveals the unidirectional bias with fast and slow dynamics for oscillation of the CNS circuit. The distribution of various oscillation behaviors can be characterized by amplitude and period, as shown in Figure 5A. We thus divide the total oscillation groups in detail according to amplitude and period. As a fact, experimental data are lacking to set the boundary of amplitude and period of different oscillation behaviors. In our study, the boundary of period is set to 100 min based on the osicllation period of NF-κB, which is also an important cell death regulator [10,73]. While the boundary of amplitude is set to 0.4, which is the median value of all the counted amplitudes. The boundary of amplitude divides the oscillation behavior into low amplitude and high amplitude. The boundary of the period divides the behavior into fast and slow oscillation. The results in Figure 5A show that the low amplitude accounts for about 37.74% and the slow period accounts for about 13.42% of the total oscillation groups. As shown in Figure 5B Figure 5C.
To investigate the discrepancy among different types of oscillation, the link between the parameter characteristics within the unidirectional bias ( Figure 4C) and the four types of oscillation are further explored. We take the top two high Frontiers in Physics | www.frontiersin.org August 2021 | Volume 9 | Article 726638 8 probability of oscillation, LF and HF as an example. The spatial parameter distribution of the LF and HF are presented in the top panel of Figure 5D (k 4 and k 10 ) and Supplementary Figure S4. We further calculate the spatial distribution difference of parameters under the two oscillation forms ( Supplementary  Information S1.2). The results of the probability distribution difference between LF and HF for each subspace are shown in the bottom panel of Figure 5D. Then, the statistical diagram of the total difference of distributions for the six parameters within the unidirectional bias ( Figure 4C) is shown in Figure 5E, indicating that parameters of k 4 and k 10 are the top two significant differences between LF and HF oscillation. Hence, from the point of dynamics mechanism, the suppression of RIP1 (represented by k 4 ) and RIP3 (represented by k 10 ) by caspase-8 are important for determining the amplitude behavior of fast oscillation.

Oscillation Behaviors Switching Controlled by the Unidirectional Bias
Having classified the four typical oscillation behaviors of the CNS circuit, we next seek to explore their switching dynamics. For simplification, we select a set of deterministic parameters in a reasonable biological range (Supplementary Information S2.1) as an example to investigate the switching of oscillation behaviors by tuning the interaction terms within the unidirectional bias. Dynamic behavior of the deterministic model upon TNF stimulation is shown by the bifurcation diagram in Figure 6A. The solid and dashed black lines respectively represent the stable and unstable equilibrium points of the system, and the colored lines represent the maximum and minimum amplitudes of the oscillations. A more detailed oscillation behavior (amplitude and period) is shown in the sub-figure, indicating that the stimulation strength of TNF for the deterministic system to oscillate ranges from 0.11 to 0.14.
To systematically study the stochastic properties of the CNS oscillations, we further employed the recently developed potential landscape theory [74][75][76][77] that describes the global dynamic behavior of the CNS system in phase space (Supplementary Information S1.2). The dimensionless potential function (U) and steady state probability distribution (P) of the system is given by the Boltzmann relation, that is, U −log(P). The system exhibits HS oscillation upon TNF 0.13 and the corresponding potential landscape of the system that mapped onto the caspase8-RIP3 phase space is shown in Figure 6B. As a result, the system evolves into a unique annular valley from any initial values, indicating that there is a unique limit cycle in the twodimensional phase space. Taking two typical evolution trajectories as examples ( Figure 6B, down panel), it is shown that the system will eventually evolve into the ring trap in any initial state. Similar results of potential landscapes that mapped onto other two-dimensional phase spaces are shown in Supplementary Figure S5.
Bifurcation analysis of all the six sets parameter within the unidirectional bias of the system is performed and the results shown in Figure 6C and Supplementary Figure S6, confirming that variations of these parameters can efficiently switch the oscillation behaviors besides k3. As an example presented in Figure 7A, the dynamic behavior of the deterministic system is switched from HS oscillation (green line) to LS oscillation (black line) by tuning the interaction term of caspase-8 activated by RIP1. However, no matter how single interaction term is tuned of the system, we could not realize the transition from fast to slow oscillation behavior, which is supported by the bifurcation diagram of all parameters shown in Figure 6C and Supplementary Figure S6. To further quantitative analysis the regulation mechanism of interaction parameter k_i in oscillation behavior, the phase diagrams in two-dimensional parameter space are shown in Figure 6D and Supplementary Figure S7. Two panels on the left of Figure 6D respectively show the phase diagrams of amplitude (upper panel) and period (lower panel) in the phase space of k 8 -k 10 , in which the amplitude is characterized by oscillation coefficient. The four typical representative time series of pRIP3 are displayed in the right panel, and they correspond to four parameter combinations in the left phase diagram. Here, the quantified values of amplitude and period of oscillation behavior are marked in the figure, also only indicating that oscillation has been switched between HS and LS.
Next, we try to tune the interaction of two terms at the same time to obtain the switching of the system from HS oscillation to LF or HF oscillation. As the simulation results shown in Figures  7B,C, by simultaneously tuning the two interaction terms of RIP1 suppressed by caspase-8 and RIP1 phosphorylated by RIP3, the HS oscillation is converted into HF (yellow line) or LF oscillation (red line). We further perform the phase diagram of the system with these two terms (Supplementary Figure S8). There are two discontinuous oscillation regions, one slow oscillation region and one fast oscillation region, revealing that the CNS oscillation circuit has strong period robustness and relatively weak amplitude robustness. In addition, the landscapes of the system after corresponding switching are shown in Figures  7D-F. The characteristics of the oscillation behavior can be displayed globally by quantifying the landscape topography. The size of the potential well on the two-dimensional phase plane quantitatively characterizes the amplitude. While the depth and breadth of the potential well reflect the stability and attraction domain of the limit cycle attractor [78].

Strong Period and Weak Amplitude Robustness of the CNS Oscillation
In the previous section, we have explored the switch mechanism between different oscillation behaviors using a deterministic system, implying that the fast oscillation behavior of the system has strong robustness. We next comprehensively investigate the robustness of the CNS oscillation network. We randomly select 100 systems from each of the four types oscillation behaviors. The value of the parameters in the unidirectional bias are tuned sequentially and continuously in the whole range, and the maximum value of oscillation amplitude change and the maximum value of oscillation period change are recorded. The violin plots in Figure 8 show the effects of the six parameters within the unidirectional bias on the amplitude and period of the four types of oscillation. The wide area in the violin Frontiers in Physics | www.frontiersin.org August 2021 | Volume 9 | Article 726638 9 plot means that the data density in this region is higher, and the white dots in the box diagram embedded in the violin plot are the median of the data.
The statistical results indicate that the period robustness of the fast oscillation systems (LF and HF) is stronger than that of the slow oscillation systems (LS and HS), but the amplitude robustness is weaker. For HF oscillation system, the amplitude robustness is the weakest, especially for the statistical results of amplitude sensitivity of k3, k4, and k6, which present two obvious high-density regions. This suggests that most HF oscillation systems are easy to switch to LF oscillation by tuning these interactions. However, owing to the strong period robustness, it is difficult to switch HF oscillation to HS or LS. For LS oscillation system, the amplitude robustness is the strongest among the four types, but its period robustness is the weakest. Thus, we can switch LS oscillation to LF by tuning a single interaction term, but it is difficult to switch to HS or HF. The results of other relative parameters are shown in Supplementary Figure S9, supporting the conclusion that the robustness of fast oscillation systems (LF and HF) is stronger than that of the slow oscillation systems (LS and HS). Therefore, these results provide potential guidance for efficiently switching oscillation behaviors among various types.

DISCUSSION
Oscillations are important characteristics for regulating many cellular physiological processes, which are conducive to stable diffusion of biological signals in vivo. In this study, we propose a circuit of the CNS pathway and further reveal the essential components and pivotal interaction terms for oscillation. Stochastic parameter analysis indicates that the essential structure for oscillation is constituted by a paradoxical component embedded with positive feedback among the three protein nodes, i.e., RIP1, caspase-8, and RIP3 ( Figure 2D). The negative feedback of caspase-8 cleavage by pRIP3 could provide an indirect time delay effect within the paradoxical component, which is consistent with the previous study that a paradoxical component with time delay could induce oscillation [27]. Most recently, we have found that the paradoxical component is also necessary for the biphasic roles of RIP1 in necroptosis [66], indicating that the paradoxical component is a core module for the CNS pathway to induce various dynamic behaviors.
We also explore how each component and interaction term within the circuit mediate the oscillation occurrence probability. Our results suggest that high levels of RIP1 and caspase-8 accelerate the occurrence of oscillation, while high level of TRADD reduce the probability. Thus, the oscillation of CNS pathway likely occurs in the cell types with high expression levels of RIP1 and caspase-8, and low expression level of TRADD. For efficiently inducing the occurrence of oscillation, we show that a unidirectional bias with incoherent fast and slow dynamics is required. Actually, besides the cell death type of necroptosis, oscillation has been proved to be important for propagating the death signaling of apoptosis and ferroptosis. Our results provide a valuable guidance for the finding of oscillation behaviors in CNS pathway, which will further advance our Biological oscillation signals contain the information of amplitude and period. Different types of oscillation signals selectively activate the downstream signaling and thus performing different biological functions. Our recent study has uncovered that the cytosolic calcium regulates apoptosis mainly through the oscillation amplitude, rather than period [14]. In this study, four types of oscillation behaviors have been classified within the CNS circuit. We have determined their probability distributions and discussed the robustness of these behaviors. The amplitude and period robustness of different oscillation behaviors are significantly different. We also evaluate the possible switching strategies among different oscillation behaviors. Our results of potential landscape provide a physical and quantitative explanation for the mechanism of modeswitching behavior of the CNS oscillation network. An urgent question is that how different oscillation behaviors can activate the downstream signaling, such as MLKL, thus determining the diverse cell fates. Salazar et al. proposed that oscillation signal is more effective than stabilization signal when kinase has a low affinity [79]. Hansen et al. found that the information transmitted via amplitude modulation signal is more reliable than that transmitted by frequency modulation signal in stress-responsive of yeast transcription factors Msn2 [80]. We suspect that RIP3 oscillation dynamics may transmit necrosis signals downstream more effectively, and MLKL can decode the information according to the amplitude and period of RIP3. However, modeling analysis of the effects of different oscillation behaviors on biological functions is still a challenge because of the lacking of effective experimental data. With the combination of further experimental observations, a more comprehensive model can be considered to dissect this urgent issue in the future. Moreover, emerging evidences have suggested the important role of oscillatory behavior in drug delivery [81,82]. Our study provides clues for efficiently switching different oscillation behaviors, which is relevant to the development of more effective strategy for medical application.
The spatial structures of oscillation model might play important roles in the cell necroptosis signal pathway. In 2018, Ferrell et al. observed that apoptosis of Xenopus laevis eggs spread through trigger waves [20]. Specifically, intracellular caspase family proteins are activated and then come into contact with adjacent caspase by floating and activating them. The apoptotic trigger waves with the speed of ∼30 μm per minute or even faster, while the average diameter of the normal cell is 10-20 μm. Thus, the spatial spread time of the proteins signaling is quite fast compared with the time from the initiation of programmed cell death to the end of ∼12 h [83]. In our study, the necroptosis signal pathway is highly correlated with the apoptosis signal pathway [61,84]. The occurrence time from the initiation to the end of necroptosis is about 6-10 h. Thus, the spatial structure of intracellular molecules can be approximated by the meanfield relative to the time scale of the cell death signaling. Actually, most protein-protein interaction network models were modeled in 2dimensions, which can make great achievements in explaining experimental phenomena and predicting biological functions [9,14]. This possibly due to the different time scales between proteins activation time and the spatial spread time. In our future work, we will decide to consider the spatial structures of the model and hope to find more important results.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
FX and XL conceived the study. ZY, LZ, and JJ participated in the modeling and provided useful suggestions. FX performed numerical simulations and QH provided technical support. FX and XL wrote the manuscript. JS gave a lot of guidance and support. All authors discussed the results of the manuscript.