Prompted by published experimental reports, which indicate that CK2 protein interacts with clock proteins in circadian rhythm (Yu et al., 2006; Szabó et al., 2013), we incorporated the molecular signaling pathways of CK2 in our model. The proposed schematic circadian–CK2 integrative model (Figure 1) is an extension of the circadian rhythm model in Drosophila, reported by Gonze et al. (2002) and Albrecht (2010). This model is based on the inhibition of a nuclear clock protein ( P N ) at the level of transcription of its gene into mRNA ( M P ) (Albrecht, 2010). mRNA is synthesized in the nucleus and is transferred to the cytosol, where it accumulates at a maximum rate of k 9 and gets ubiquitinated by enzymes (E1, E2, E3, and E4) with a rate of k 12 . In our model, the spatial temporal fluctuations in the concentrations of various forms of the nuclear clock protein ( P N ) or cytosolic regulatory protein ( P 0 , P 1 and P 2 ) is governed by the biochemical pathways. The rate of synthesis of protein P 0 is proportional to the formation of M P , and is characterized by an apparent first order rate constant k 13 . Parameters k 16 and k 22 denote the maximum rate(s) and Michaelis constant(s) of the phosphatase and kinase involved in the reversible phosphorylation of P 1 into P 2 and P 2 into P 1 , respectively. The fully phosphorylated state ( P 2 ) is degraded by CK2 and is transported into the nucleus at a rate characterized by the apparent first-order rate constant k 20 . Transport of the P N into the cytosol is characterized by the apparent first-order rate constant k 30 . It is reported that the negative feedback is exerted by the nuclear clock protein on gene transcription. CK2 phosphorylates the clock protein, which may form three complexes, namely CK 2 − P 0 , CK 2 − P 1 , and CK 2 − P 2 , due to three different forms of the available proteins. Synthesis of CK 2 − P 0 is assumed to occur with a rate constant of k 32 , and subsequently, the dissociation of this complex is assumed to occur with a rate constant of k 33 . In a similar manner, it is hypothesized that complex formation of CK 2 − P 1 occurs with a rate constant of k 34 and its dissociation follows with a rate constant of k 35 . Further, complex formation of CK2 − P2 is considered to occur with a rate constant of k 36 , and subsequently, the dissociation of this complex is assumed to occur with rate constants of k 37 . CK2 is ubiquitinated at a rate constant of k 38 . Synthesis of CK2 in the network is assumed to occur at the rate constant of k 38 . CK2 plays a significant role in various processes of the organism such as biological clocks. Table 1 provides the list of clock proteins linked to the integrated model and Table 2 presents the list of the biochemical reactions and propensity functions (probability of reaction). The rate constants associated with the proposed model are presented in Table 3.

Complex dynamical processes, governed by a set of well-defined reaction channels, are generally noise-induced stochastic processes due to random molecular interactions in the system (origin of intrinsic noise) and continuous interaction of the system with random environmental fluctuations (origin of extrinsic noise) (Doob, 1942; Gillespie, 1977; Miyazaki et al., 2004). The system with N clock protein variables, whose population vector is defined by, X → ( t ) = ( X 1 , X 2 , . .. , X N ) T , where T is the transpose of the vector, which undergo M reaction channels, R I , i = 1 , 2 , . .. , M is given by,

where { k i } , i = 1 , 2 , . .. , M is the set of classical rate constants. The stochastic rate constant c i of any of its reactions can be expressed in terms of classical rate constant k i as c i = k i V 1 − ν , where V is the system size and ν is the state change parameter of its reaction (Gillespie, 1977; Gillespie, 2000). The trajectories of the variables provided by birth and death processes due to molecular interaction given by the Equation 1 can be traced by solving the master equation, for the rate of change of configurational probability P ( X → , t ) as a function of time. The master equation is denoted as

where W and W ′ are the transition probabilities of the two configurational states X → and { X → ′}. Solving the master Equation 2 analytically for complex biological process is extremely difficult. However, numerical solutions for the master Equation 2 can be derived using stochastic simulation algorithm (SSA) (Gillespie, 1977; Gillespie, 2000) that is based on the theoretical foundations developed by Doob (1942) and initially proposed by Kendall (1950). The stochastic simulation implements a Monte Carlo algorithm that provides the exact numerical solution by considering every possible interaction in the system (Gillespie, 1977; Singh et al., 2018). This algorithm is in deed a non-spatial, individual-based analog of the master Equation 2, which is constructed on the physical basis of molecular collision in each reaction channel at a certain constant temperature.

This SSA is based on two crucial independent random processes, namely reaction fire and reaction time. These two processes are implemented in this algorithm by generating two statistically independent random numbers, namely r 1 and r 2 , such that the reaction time is computed using τ = − 1 a 0 ln r 1 , where a 0 ∑ i a i , a i is the i t h propensity function given by a i = h i c i , where h i is the number of possible molecular combinations of i t h reaction, and the k t h reaction will fire when it satisfies ∑ i = 1 k a i ≤ a 0 r 2 < ∑ i = 1 k + 1 a i . Intrinsic noise ( ξ ) associated with the clock protein dynamics in the system is inversely proportional to the square root of the systems size V (i.e., ξ ∝ 1 V ) (Gillespie, 2000; Singh et al., 2018; Sharma et al., 2019; Singh et al., 2021).

Permutation entropy can be used to measure the complexity of a system associated with the dynamics of the system’s variables (Bandt and Pompe, 2002). The basic algorithm for calculating permutation entropy ( H ) of a time series is as follows:

Consider a dynamical variable x ( t ) of a system given by the discrete time series = { x i } , i = 1 , 2 , . .. , N , where N is the finite total number of discrete time elements in the time series data. We define an embedding dimension ( d ), preferably d = 3 , . .. , 7 , to represent the data in consecutive patterns of the size dimension d . For a particular value of d , there are M possible permuted sequences of inequalities of sequence elements. If we take d = 3 , we will have M = 6 arrangements of permutations given by,

u 1 = { x 1 , x 2 , x 3 } , u 2 = { x 1 , x 3 , x 2 } , u 3 = { x 2 , x 1 , x 3 } , u 4 = { x 2 , x 3 , x 1 } , u 5 = { x 3 , x 1 , x 2 } , u 6 = { x 3 , x 2 , x 1 } , where, x 1 ≠ x 2 ≠ x 3 . Then, we calculate the probability ( p u ) for each single permutation sequence, which is the ratio of number of values for a particular sequence of permutations to the total number of all possible permutations for the embedding dimension, d , in the time series data. Shannon entropy for a sequence of perturbations can then be calculated by Equation 3,

and permutation entropy of embedded dimension (d) by the sum of these entropies as mentioned in Equation 4,

where 0 ≤ H ≤ 1 . The mapped permutation entropy spectrum of time series x ( t ) is indicated by H = { H 1 , H 2 , . .. , H N } , and exhibits a behavior similar to that of Lyapunov spectrum of the same time series (Bodine et al., 2001).

We first present the results of how CK 2 regulates the dynamics of clock protein ( P N , M P ). Since the population of the CK 2 protein in the system of size V (at fixed V = 200 ) is proportional to the rate of synthesis of this protein in the system ( CK 2 ∝ k 31 ), variation in CK 2 might cause changes in the interaction rate of other clock proteins in the reaction network model. Hence, we looked for variations in the dynamics of the clock proteins, driven by CK 2 via changes in the values for k 31 . First, we allowed all three phosphorylation events in the cytosolic PER proteins (well-known clock proteins), P 0 , P 1 , P 2 , with CK 2 to occur (Figure 1), and found that for small values of k 31 = 0.001 , prominent oscillation (active state) in the nuclear PER protein is exhibited with time period of oscillation, T P N ~ 23.4 ± 0.13 h and amplitude of oscillation A P N ∼ 285 ± 3 (Figure 2A upper panel). Further increase in k 31 suppressed the oscillation ( k 31 = 0.04 , 0.08 ) , allowing increased T P N and decreased A P N significantly, which could correspond to weak circadian activity (Crino, 2011). This increase in the time period was due to increased phosphorylation owing to an increase in interaction of CK2 with the PER protein ( P 0 , P 1 , P 2 ), which was evident from the experimental reports of Cao et al. (2009). Further, weak circadian activity may cause various diseases, such as aging of the brain, metabolic dysfunction, dementia, and cancer-related disorders (Rao et al., 2002). When k 31 is sufficiently large ( k 31 ~ 1.0 ), the P N dynamics showed both oscillation and amplitude death scenarios (Rao et al., 2002). This state may correspond to circadian rhythm death (Rao et al., 2002). These three important states in the dynamics of P N of the circadian rhythm induced by CK2 can be shown in two-dimensional space of parameter as ( k 31 , A P N and k 31 , T P N ) (Figure 2B). Here, the parameters A P N and T P N are mean values of P N amplitude and time period for a range of t = [ 10 − 500 ] hours. In the phase diagram presented in Figure 2B, all three circadian states are seen to be clearly demarcated.

We delineated the measure of complexity in the three derived states (active, weak activity, and rhythmic death) driven by CK 2 by calculating permutation entropy ( H P N ) of the dynamics of P N at the three respective states (Figure 2C). Upon considering the average permutation entropies of active, weak activity, and rhythmic death as H P N , A H P N W A , and H P N R D , the simulation results demonstrated that H P N A > H P N W A > H P N R D . The three circadian rhythm states (active, weak activity, and rhythmic death) could easily be classified using H P N A (Figure 2D). In this phase diagram, each point is the average of permutation spectrum for each value of k 31 with error bars. Hence, the results indicated that H P N could be used as a parameter to classify various states of the circadian rhythm.

The dynamics of M P exhibited similarity to that of P N at fixed system size and variable concentration of CK2. The dynamics of M P for various values of k 31 for constant values of system size is presented in Figure 3A. We measured the permutation entropies of MP at the three states under the situation of fixed system size and of variations of the CK2 concentration. In the case of permutation entropies of MP at the three states (Figure 3B), we made observations similar to that of P N (as depicted in Figure 2B).

We propose that H P N and H M P can be useful for clinicians and medical practitioners to classify the various states of circadian rhythm. These results suggest that CK2 plays a very dynamic role in the regulation of circadian rhythm in living systems. The change in rhythmic properties may alter the physiological processes of the organism leading to various diseases.

Figure 4 presents the various possible interactions of CK2 with each of the PER proteins ( P 0 , P 1 , and P 2 ) to form complexes CK 2 − P 0 , CK 2 − P 1 , and CK 2 − P 2 . All simulations were performed for the same values of k 31 and V . used in the previous simulation depicted in Figure 3. The results demonstrated that if CK 2 interacts individually with any one of the PER proteins via formation of a complex, P 2 dynamics exhibits both active (for smaller values of k 31 ) and weak activity (for larger values of k 31 ) states (Figure 4A). However, the system needs significantly large values of k 31 to obtain the rhythmic death state (data not shown). Further, configurational interaction of C K 2 with P 2 is much more sensitive in driving the circadian rhythm states, as compared to those with P 1 or P 0 . We then calculated A P N and T P N for each time series when CK2 was allowed to interact with each one of P 0 , P 1 , and P 2 in the circadian–CK2 model (Figure 4B). Each curve in A P N and T P N (Figure 4B) illustrates that the three circadian rhythm states could be observed. We also found that interaction of CK2 with P 2 is more sensitive than that with either P 0 or P 1 because the three circadian states are found at smaller values of k 31 in the case of CK 2 − P 2 complex .

The three models, namely model1, model2, and model3, as presented in Figure 4, are patterns of weak activity when CK2 systematically interacts with any one of the PER proteins ( P 0 , P 1 , and P 2 ); these three patterns can also be considered as pathological states of the corresponding configuration. This rationale is acceptable because the patterns of CK2 activities (namely amplitude and time period) are significantly changed due to drastic increase in phosphorylation of PER protein(s) with CK2 (Cao et al., 2009). This drastic change in circadian rhythm may cause various diseases ranging from socio-psychological diseases (Wulff et al., 2010) and metabolic syndrome (Maury et al., 2010), to various types of cancer (Savvidis and Koutsilieris, 2012). Excess phosphorylation of CK2 with PER protein switches the pathological state to rhythmic death pattern which could be apoptosis signature.

We then calculated permutation entropy ( H P N ) of the time series of P N as a function of k 31 for individual interactions of CK2 with P 0 , P 1 , and P 2 (Figure 4C). We found that the three circadian states can be detected distinctly for each time series, as in the previous simulation, and H P N can be used as a parameter to detect these states which can be useful to clinicians and medical practitioners (Figure 4D).

The dynamics of MP (Figure 5A) exhibited similarity to that of P N (as seen in Figure 4A). We also measured the permutation entropies of MP at fixed CK2 concentration and variation in system size. The observations from permutation entropies of MP at the three dynamical states of MP (Figure 5B) were similar to those of P N .

Noise is an inherent property associated with the dynamics of any natural system (Glass, 2001; Bhadana et al., 2019, 2021). We studied the dynamics of P N for four possible configurational interactions of CK2 with the cytosolic PER proteins ( P 1 , P 2 , a n d [ P 0 + P 1 + P 2 ] ) , as we did in the previous simulation by keeping k 31 fixed at k 31 = 0.01 (which is associated with regular pattern of active circadian rhythm), and changing the strength of the noise ξ ∝ 1 V with four values for V as 80 , 100 , 200 , and 500 (Figure 6). In all the configurational interactions of CK 2 with the PER protein, we observed that an increase in noise (i.e., decrease in the value of V ) can alter the P N dynamics in each of the three circadian states. This indicates that noise can regulate the dynamics of circadian rhythm.

Generally, the size of an animal cell (of both the proliferation and non-proliferation types) can change due to various reasons (Lloyd, 2013), such as cell cycle progression and cell growth (Conlon and Raff, 2003), pathological states in muscle cells and starvation (Bodine et al., 2001), tumor or cancer progression (Lum et al., 2005), defects in synaptic wiring/rewiring in neurons (Crino, 2011), and manipulation of extracellular signals to prevent apoptosis (Lum et al., 2005). Animal cells can have size variability of up to ~ 59 % of the normal cell size (Lum et al., 2005), and this change in the cell size can drastically affect molecular crowding within the cell (Feig et al., 2017). This variation in the system’s size is reflected in the dynamics of the system’s variables as internal noise fluctuation ( ξ ∝ 1 V ) (Gillespie, 2000). Since the circadian rhythm system is governed by well-defined reaction channels (Table 2), this change in molecular crowding could cause two significant impacts on the system. Firstly, it allows changes in the rate of interaction of the clock proteins in the system. This change in molecular interaction leads to changes in the internal noise associated with the system reflected in the dynamics of the constituting variables P N , M P , and of other clock proteins in the system. Hence, excess noise may destroy the signal associated with the system variables and may lead to external collapse of the system itself. Secondly, this molecular crowding may trigger changes in molecular traffic in the system; traffic jam due to excess molecular crowding may make the system unable to perform normal functions. In such a scenario, the amplitude is minimized, and the system faces death with infinitely large time period (Figure 6). Similar behavior can be observed in the other three configurational interaction models of CK2 with each of P 0 , P 1 , and P 2 , particularly, in the interaction of CK2 with P 2 wherein the noise was slightly more sensitive (Figure 6A).

Results of calculations to derive the amplitudes A P N and time periods T P N of P N dynamics as a function of V for four different configurational interaction of CK 2 with P 0 , P 1 , and P 2 are presented in Figure 6B. The results indicated that noise, measured by V , could drive the system to the three distinct states, namely active, weak activity, and rhythmic death. Hence, it can be inferred that noise is an important parameter that can trigger the system at various rhythmic states and can regulate the dynamics of the system. Permutation entropy spectra of P N corresponding to all the four models as a function of V values are presented in Figure 6C. We delineated the measure of complexity in three different derived states driven by CK 2 by calculating permutation entropy H P N of the dynamics of P N for all possible configurational interaction of C K 2 with cytosolic PER proteins ( P 0 , P 1 , P 2 , and P 0 + P 1 + P 2 ) . Further, it can be observed that the measure of H P N as a function of V can distinctly classify the three different circadian rhythm states (Figure 6C). Upon denoting H P N , A H P N W A , and H P N R D as permutation entropies corresponding to the circadian states of active, weak activity, and rhythmic death, it can be inferred from the results that H P N A > H P N W A > H P N R D (Figure 6D).

A similar behavior can be seen in the case of M P (Figure 7A) and H MP (Figure 7B) dynamics. The transition from sustained oscillation fluctuations (active state) to amplitude death (rhythmic death) via the weak activity state can be clearly seen in the two-dimensional plot ( P N , M P ) as presented in Figure 7C.

We then identified cellular pathways that could be triggered by PER mutants (P₀, P_1, and P₂) in the circadian rhythm model (Figure 8). These cellular pathways can be grouped into three categories, cancer pathways (Savvidis and Koutsilieris, 2012; Papagiannakopoulos et al., 2016; Lee et al., 2019), socio–psychological pathways (Foster et al., 2013), and metabolic pathways (Maury et al., 2010). PER genes with their own rhythms generally regulate expressions of other genes participating in various important cellular pathways, significant changes of which may lead to various diseases. Since disruptions in the dynamics of P 0 , P 1 , P 2 , and P N in circadian rhythm may lead to drastic changes in the dynamics and mechanisms of these closely interacting pathways and to various pathological states, they can cause various diseases. Dysfunctional PER proteins disturb key biological functions, such as cell proliferation, DNA damage, cell cycle, and apoptosis, resulting in various type of cancers (Deibel et al., 2015).

According to various studies, there may exist a tissue-dependent relationship between insulin resistance (or type 2 diabetes, T2D) and dysregulated molecular clock activity (Hart, 2013; Gabriel et al., 2021). In a time-course experiment (Jakubowicz et al., 2017) with white adipose tissue biopsies from people with healthy weight or obesity or T2D, no variation in the rhythm and amplitude of core-clock (PER1, PER2, PER3, DBP, BMAL1, and CRY2), metabolic (PGC1), and clock-related (REVERB) genes could be observed in biopsy tissues (Jakubowicz et al., 2017). On the other hand, when the sleep–wake cycle and dietary factors were controlled, it was seen that the amplitude oscillations of core-clock genes and a number of rhythmic genes are decreased in adipose tissues from patients with T2D, as compared to lean and healthy individuals (Stenvers et al., 2019). The mRNA expression of BMAL1, PER1, PER2, and PER3 in leukocytes was observed to be lower in non-diabetic individuals, as compared to those with diabetes (Hart, 2013; Stenvers et al., 2019). Additionally, the expression levels of PER1, PER3, and BMAL1 molecular-clock genes in leukocytes obtained from patients with T2D was seen inversely correlated with hemoglobin A1C (HbA1c) levels, indicating an association between insulin resistance and T2D (Gabriel et al., 2021). Expression of PER2, PER3, and CRY2 mRNA were also significantly associated with plasma HbA1c levels and islet insulin content in pancreatic islets from healthy people and from those with T2D (Gabriel et al., 2021). Uncertainties still exist regarding the fundamental mechanisms that regulate metabolic rhythmicity, especially regarding whether rhythmicity is lost in T2D. Mutations in clock genes were initially associated with glucose homeostasis (Rudic et al., 2004) and later with hyperinsulinemia, hyperglycaemia, and obesity in murine models (Turek et al., 2005). Numerous single nucleotide polymorphisms (SNPs) in human clock genes, including rs18012602, rs4580704, rs4864584, rs3749474 and rs1464490, have also been associated with obesity (Sookoian et al., 2008), hyperglycaemia, and a higher prevalence of T2D (Sookoian et al., 2008).

Aberrations in the PER protein disrupts multiple biological functions, such as age-related hydroxy methylation, which causes increased risk for socio–psychological diseases (e.g., learning and memory impairment, and Alzheimer’s disease) (Liu and Chang, 2017). Significant changes in the dynamics of P 0 , P 1 , P 2 and P N in the circadian rhythm and noise in the system may lead to drastic changes in the dynamics and mechanisms of pathways and to various pathological states, causing diseases corresponding to these pathways. Therefore, it is crucial to identify pathways that are significantly affected by changes in the P 0 , P 1 , P 2 , and P N dynamics.

On the other hand, patients having any one or more of the diseases mentioned in Figure 8 can become healthy by maintaining a proper circadian rhythm (Foster et al., 2013), and the process could possibly be termed as a reversal of pathological state. However, even though this may be possible for some diseases, one needs to categorically study these possibilities from both the mathematical modeling and experimental points of view. In accordance with these postulations, previous reports have also recommended studies on genes/proteins in the circadian rhythm pathways as potential drug targets (Wood et al., 2009).