Fractional-Order Dynamical Modeling and Stability Analysis of the Core Circadian Clock in Arabidopsis thaliana

1. Introduction

Arabidopsis thaliana is a widely used model organism in plant molecular biology, genetics, and systems biology. Owing to its short life cycle, abundant genetic resources, and well-established experimental system, it has been extensively used to investigate plant growth and development, light signaling, environmental adaptation, and circadian regulation [1]. The plant circadian system is an endogenous timekeeping mechanism that coordinates physiological processes such as photosynthesis, flowering time, hypocotyl elongation, root growth, metabolic regulation, and responses to environmental stresses over an approximately 24 h cycle, thereby maintaining appropriate temporal relationships with external light–dark cycles [2,3,4]. In A. thaliana, the core circadian clock is driven by a transcription–translation feedback network composed of multiple clock genes and their encoded proteins. Core components, including CCA1/LHY, PRR9/PRR7, PRR5/TOC1, and the evening complex formed by ELF3, ELF4, and LUX, are expressed at different times of day and sustain rhythmic oscillations through interconnected activation, repression, and feedback regulation [2,4].

The plant circadian network involves multiple layers of feedback regulation, light input, transcription and translation, protein degradation, and downstream physiological responses. Qualitative descriptions alone are insufficient to fully reveal the dynamic relationships among these regulatory processes. Mathematical modeling has therefore become an important approach for investigating the mechanisms underlying plant circadian rhythms. Locke et al. developed early ordinary differential equation models of the A. thaliana circadian clock by integrating feedback regulation among core clock genes with experimental expression data, thereby establishing a foundation for the quantitative investigation of plant circadian rhythms [5,6]. Subsequently, Pokhilko et al. developed a plant circadian clock model containing a repressilator architecture and multiple feedback loops to explain the temporal expression patterns and phase-response characteristics of core clock components [7]. Fogelmark and Troein reanalyzed extensive publicly available time-series datasets for the A. thaliana core clock and emphasized the role of transcriptional repression in sustaining rhythmic oscillations of the core clock [8]. These models progressively improved the representation of the molecular regulatory mechanisms of the plant circadian clock. However, as the number of model components and regulatory interactions increased, the numbers of state variables and parameters to be estimated also increased, substantially increasing the complexity of parameter fitting and dynamical analysis.

To balance biological interpretability and computational tractability, De Caluwé et al. [9] proposed a compact model of the A. thaliana core circadian clock. This model groups core clock genes with similar expression patterns and regulatory functions into composite variables and uses a relatively small number of state variables to describe the mRNA and protein dynamics of the CCA1/LHY, PRR9/PRR7, PRR5/TOC1, and ELF4/LUX modules. Compared with more complex high-dimensional plant circadian clock models, the compact model has a simpler structure and fewer parameters, while retaining the ability to reproduce the major rhythmic characteristics of the wild type under entrained and free-running conditions, as well as changes in period, phase, and expression level in several mutants [9]. A recent systematic review of mathematical models of the A. thaliana circadian clock indicated that the field has gradually shifted from model expansion through the addition of regulatory components toward model reduction and the identification of core dynamical modules. Improving model analyzability while preserving the principal biological mechanisms remains an important research direction [10].

Building on the compact model, subsequent studies have incorporated spatial organization, light-quality responses, and downstream growth phenotypes. Greenwood et al. embedded the De Caluwé model within a spatial representation of seedlings to investigate how tissue-specific light sensitivity and local intercellular coupling affect circadian period, phase, and spatial expression waves [11]. Pay et al. [12] incorporated red- and blue-light receptors and a COP1-related light-signaling module into a core clock model and developed a model for describing hypocotyl growth under different light qualities and photoperiods. In addition to extending models by introducing additional state variables and regulatory pathways, delay-based formulations provide another important approach for representing temporal accumulation in biological regulatory processes. Tokuda et al. [13] introduced distributed delays into an A. thaliana circadian clock model, reducing parameter complexity in some model formulations while preserving the core regulatory mechanisms and principal dynamical features. These studies indicate that plant circadian clock models are evolving from descriptions focused mainly on core gene-expression oscillations toward integrated dynamical frameworks that account for spatial coordination, environmental inputs, time delays, and downstream phenotypic outputs.

Despite these advances, most existing plant circadian clock models are still formulated using integer-order ordinary differential equations. Integer-order differentiation is local in time, such that the instantaneous rate of change of a system is generally determined by its current state and current external inputs. By contrast, biological processes such as transcriptional regulation, protein accumulation and degradation, and responses to light–dark signals may involve multistep reactions, cumulative effects, and delayed responses. Explicit and distributed delays describe such historical influences through specified time delays or memory kernels, whereas fractional calculus continuously weights the past trajectory of a system through the nonlocal integral structure of a non-integer-order derivative [13,14]. In particular, the Caputo fractional derivative is nonlocal and permits initial conditions of the same form as those used for integer-order differential equations, making it a feasible phenomenological tool for modeling biological systems with historical dependence [14].

Fractional-order approaches have been applied to a wide range of biological dynamical systems. In epidemic dynamics, a fractional-order MSEIR model has been used to represent historical dependence in infection, latency, recovery, and immunity processes [15]. In ecological modeling, fractional-order predator–prey models have been used to investigate memory effects associated with Allee effects, group defense, and interactions between populations [16,17]. Fractional-order models have also been applied to tumor–immune systems and HIV infection dynamics [18], as well as to a broader range of biological and bioengineering systems [19]. In circadian-rhythm research, Escalante-Martínez et al. introduced the Caputo fractional derivative into coupled van der Pol oscillators to analyze the synchronization of circadian oscillators [20]. De Cezaro et al. subsequently developed a fractional-order model of coupled oscillators with different memory levels and investigated the effects of memory distribution on circadian synchronization [21]. These studies indicate that the fractional order can serve as an additional parameter governing the effective strength of historical dependence and the dynamic response characteristics of a system. However, to the best of our knowledge, existing fractional-order studies of circadian rhythms have mainly focused on abstract oscillators and their synchronization. Studies that directly address the transcription–translation feedback network of the A. thaliana core clock and combine expression time-series and mutant data for parameter estimation and validation remain limited.

Motivated by this gap, we extend the compact model of the A. thaliana core circadian clock proposed by De Caluwé et al. by replacing the integer-order derivatives with Caputo fractional derivatives, thereby constructing a fractional-order dynamical model of the core circadian clock. We first present the mathematical formulation of the model and establish the local existence and uniqueness, non-negativity, and boundedness of its solutions under continuous darkness. We then fit the model parameters to wild-type mRNA expression data under a 12 h light/12 h dark photoperiod (12L:12D) and, without refitting the parameters, validate the model using expression data obtained under other photoperiods and from mutants of major clock genes. Next, we analyze the local stability of the equilibrium under continuous darkness, investigate the effects of the fractional order on system stability, oscillation amplitude, and the phases of the core clock modules, and examine the joint effects of the fractional order and protein degradation rates on the stability region. This study aims to evaluate the applicability of fractional-order dynamics to modeling the A. thaliana core circadian clock and to provide a new dynamical framework for plant circadian systems with potential nonlocal historical dependence. The results show that the fractional-order formulation improves the description of selected post-peak decay and tough phases, yields lower condition-specific MSE values than the original ODE model under all examined wild-type light conditions and selected mutant conditions, and reveals an order-dependent transition in equilibrium stability.

2. Preliminaries and Model Formulation

Fractional integrals and derivatives provide important mathematical tools for describing dynamical systems with nonlocal dependence on their past states. To facilitate the formulation and subsequent dynamical analysis of the fractional-order core circadian clock model, we first introduce the Riemann–Liouville fractional integral, the Caputo fractional derivative, and the Mittag–Leffler function. We then formulate a fractional-order dynamical model of the Arabidopsis thaliana core circadian clock based on the compact model proposed by De Caluwé et al. [9].

2.1. Preliminaries

Definition 1([22]). Let $a<T$ , $q>0$ , and $u\in L^1(a,T)$ . The Riemann–Liouville fractional integral of $u$ of order $q$ is defined by

$$\left(I_{a+}^{q}u\right)\left(t\right)={\displaystyle \frac{1}{\Gamma \left(q\right)}}{\int }_a^t(t-s{)}^{q-1}u\left(s\right)ds,fort\in (a,T).$$

Here,  $\Gamma (\cdot )$  denotes the Gamma function, defined for  $z\in \mathbb{C}$  with  $Re\left(z\right)>0$  by

$$\Gamma \left(z\right)={\int }_0^{\infty }{r}^{z-1}{e}^{-r}dr.$$

In the initial-value problems considered in this study, the initial time is set to $a=0$; therefore, the notation $I_{0+}^q$ is used throughout the remainder of this paper.

Definition 2([23]). Let $q>0$ , and let $n=\lceil q\rceil$ , so that $n-1<q\le n$ . If $n-1<q<n$ and $u\in A{C}^n\left(\right[a,T\left]\right)$ , then

$${}_a{}^C{D}_t^{q}u\left(t\right)={\displaystyle \frac{1}{\Gamma \left(n-q\right)}}{\int }_a^t(t-s{)}^{n-q-1}{u}^{\left(n\right)}\left(s\right)ds,fort\in (a,T),$$

is called the Caputo fractional derivative of  $u\left(t\right)$  of order  $q$  , where  $A{C}^n\left(\right[a,T\left]\right)$  means that  $u^{\left(n−1\right)}\left(t\right)$  is absolutely continuous on  $[a,T]$.

When  $q=n$  is a positive integer, the Caputo fractional derivative reduces to the classical integer-order derivative, namely,

$${}_a{}^C{D}_t^{q}u\left(t\right)=u^{\left(n\right)}\left(t\right).$$

In this study, we mainly consider the case  $0<q\le 1$  , for which  $n=1$.

When  $0<q<1$  and  $u\in AC\left(\right[0,T\left]\right)$  , the Caputo fractional derivative can be written as

$${}_0{}^C{D}_t^{q}u\left(t\right)={\displaystyle \frac{1}{\Gamma (1-q)}}{\int }_0^t(t-s{)}^{-q}{u}^{\text{'}}\left(s\right)ds,0<t\le T.$$

When $q=1$ , the Caputo fractional derivative reduces to the classical first-order derivative, namely, ${}_0{}^C{D}_t^{1}u\left(t\right)={\displaystyle \frac{du\left(t\right)}{dt}}$

It follows from the above definition that the Caputo fractional derivative at time $t$ depends on historical information over the interval $[0,t]$ and is therefore nonlocal. Compared with the Riemann–Liouville fractional derivative, the Caputo fractional derivative permits initial conditions of the same form as those used for integer-order differential equations and is therefore more convenient for modeling biological dynamical systems with explicitly specified initial states [23].

Definition 3([24]). Let $\alpha >0$ , $\beta >0$ , and $z\in \mathbb{C}$ . The function

$$E_{\alpha ,\beta }\left(z\right)={\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \frac{z^k}{\Gamma (\alpha k+\beta )}}$$

is called the two-parameter Mittag–Leffler function. In particular, when  $\beta =1$  , the two-parameter Mittag–Leffler function reduces to the one-parameter Mittag–Leffler function. When  $\alpha =\beta =1$  , we have

$$E_{\mathrm{1,1}}\left(z\right)=e^z.$$

Therefore, the Mittag–Leffler function can be regarded as a generalization of the exponential function in fractional-order dynamical systems. The Mittag–Leffler function in the form $E_q(-b{t}^q)$ will be used later in the boundedness estimates of the model solutions.

Definition 4. 

Consider the autonomous Caputo fractional-order system

$${}_0{}^C{D}_t^{q}X\left(t\right)=F\left(X\left(t\right)\right),0<q<1,$$

where  $X\left(t\right)\in {\mathbb{R}}^n$  and  $F:{\mathbb{R}}^n\to {\mathbb{R}}^n$  . If there exists a constant vector  $X^{*}\in {\mathbb{R}}^n$  such that

$$F\left(X^{*}\right)=0,$$

then  $X^{*}$  is called an equilibrium point of the system.

2.2. Model Formulation

Based on the regulatory interactions included in the compact model of the Arabidopsis thaliana core circadian clock proposed by De Caluwé et al. [9], we constructed a schematic regulatory network of the composite core clock gene modules, as shown in Figure 1. The compact model consists of nine ordinary differential equations describing the dynamics of four composite core clock gene modules and one light-signal-related variable [9]. It groups core clock genes with similar expression patterns and regulatory functions. Specifically, CL represents the dawn-phased CCA1/LHY module, P97 represents the morning-phased PRR9/PRR7 module, P51 represents the evening-phased PRR5/TOC1 module, and EL represents the evening-complex-related ELF4/LUX module. The variable $P$represents a light-signal-related regulatory component that accumulates during darkness and contributes to the light-induced transcription of CCA1/LHY and PRR9/PRR7 following a dark-to-light transition [9].

As shown in Figure 1, the core of the De Caluwé compact model describes the transcription, translation, degradation, and mutual regulation of the four composite core clock gene pairs. For $i=\mathrm{1,2},\mathrm{3,4}$, let $R_i$ denote the mRNA level of the $i$-th composite module and let $P_i$ denote the corresponding protein level. More precisely,

$$R_1=[\mathrm{C}\mathrm{L}{]}_m,P_1=[\mathrm{C}\mathrm{L}{]}_p,R_2=[\mathrm{P}97{]}_m,P_2=[\mathrm{P}97{]}_p,$$

$$R_3=[\mathrm{P}51{]}_m,P_3=[\mathrm{P}51{]}_p,R_4=[\mathrm{E}\mathrm{L}{]}_m,P_4=[\mathrm{E}\mathrm{L}{]}_p.$$

Here, CL, P97, P51, and EL correspond to the CCA1/LHY, PRR9/PRR7, PRR5/TOC1, and ELF4/LUX gene pairs, respectively. The subscripts $m$ and $p$ denote the mRNA and protein levels, respectively. The ninth variable $P_L$ denotes the activity of the dark-accumulated protein, a stand-in for PIL1, PIF3 and other proteins, which induce the transcription of CCA1, LHY, and PRR9 upon dark-to-light transition.

To represent the potential historical dependence and cumulative effects involved in plant circadian regulation, we improve the De Caluwé compact model [9] by replacing the integer-order derivatives with Caputo fractional derivatives. The resultant fractional-order differential equation (FDE) model is given by

$${}_0{}^C{D}_t^q{R}_1=\left(v_1+v_{1L}L\left(t\right)P_L\right){\displaystyle \frac{1}{1+{\left(\frac{P_2}{K_1}\right)}^2+{\left(\frac{P_3}{K_2}\right)}^2}}-\left(k_{1L}L\left(t\right)+k_{1D}D\left(t\right)\right)R_1,$$

$${}_0{}^C{D}_t^q{P}_1=\left(p_1+p_{1L}L\left(t\right)\right)R_1-d_1{P}_1,$$

$${}_0{}^C{D}_t^q{R}_2=\left(v_{2L}L(t)P_L+v_{2A}+v_{2B}{\displaystyle \frac{P_1^2}{K_3^2+P_1^2}}\right){\displaystyle \frac{1}{1+{\left(\frac{P_3}{K_4}\right)}^2+{\left(\frac{P_4}{K_5}\right)}^2}}-k_2{R}_2,$$

$${}_0{}^C{D}_t^q{P}_2=p_2{R}_2-\left(d_{2D}D\left(t\right)+d_{2L}L\left(t\right)\right)P_2,$$

$${}_0{}^C{D}_t^q{R}_3=v_3{\displaystyle \frac{1}{1+{\left(\frac{P_1}{K_6}\right)}^2+{\left(\frac{P_3}{K_7}\right)}^2}}-k_3{R}_3,(2.1)$$

$${}_0{}^C{D}_t^q{P}_3=p_3{R}_3-\left(d_{3D}D\left(t\right)+d_{3L}L\left(t\right)\right)P_3,$$

$${}_0{}^C{D}_t^q{R}_4=L\left(t\right)v_4{\displaystyle \frac{1}{1+{\left(\frac{P_1}{K_8}\right)}^2+{\left(\frac{P_3}{K_9}\right)}^2+{\left(\frac{P_4}{K_{10}}\right)}^2}}-k_4{R}_4,$$

$${}_0{}^C{D}_t^q{P}_4=p_4{R}_4-\left(d_{4D}D\left(t\right)+d_{4L}L\left(t\right)\right)P_4,$$

$${}_0{}^C{D}_t^q{P}_L=0.3D\left(t\right)(1-P_L)-L\left(t\right)P_L,$$

where ${}_0{}^C{D}_t^q$ denotes the Caputo fractional derivative and $\left.q\in (\mathrm{0,1}\right]$ is the fractional order. When $q=1$, the Caputo derivative reduces to the usual first-order derivative, and system (2.1) reduces to the corresponding integer-order model with the same regulatory structure.

The $v$-parameters represent basal, light-induced, or regulated transcription rates; the $p$-parameters represent translation rates; the $k$-parameters represent mRNA degradation rates; the $d$-parameters represent protein degradation rates; and $K$-parameters are the half-saturation constants or regulatory thresholds in the Hill-type regulatory functions. Both activation and repression are described using Hill-type functions. The Hill coefficient of each activation or repression function is fixed at 2, as in the original compact model, and is not treated as an estimated parameter. All model parameters are positive. Their biological meanings and numerical values are provided in Appendix A.

The functions $L\left(t\right)$ and $D\left(t\right)$ represent the external light and dark inputs, respectively, rather than system state variables. During the light phase,$L\left(t\right)=1,D\left(t\right)=0,$ whereas during the dark phase, $L\left(t\right)=0,D\left(t\right)=1.$

Thus, different photoperiods act on the system through the external input functions $L\left(t\right)$ and $D\left(t\right)$. Specifically, the light–dark inputs regulate the transcription of the CCA1/LHY, PRR9/PRR7, and ELF4/LUX mRNA modules, the degradation of CCA1/LHY mRNA, the translation of the CCA1/LHY protein, and the degradation of the PRR9/PRR7, PRR5/TOC1, and ELF4/LUX proteins.

Under a general light–dark cycle, the system is driven by the periodic external inputs $L\left(t\right)$ and $D\left(t\right)$ and is therefore nonautonomous. Under continuous darkness, $L\left(t\right)=0$ and $D\left(t\right)=1$, and the system reduces to an autonomous fractional-order system. The basic properties of the model solutions under continuous darkness are analyzed in the following section.

3. Materials and Methods

3.1. Data Sources and Preprocessing

The mRNA time-series expression data used for parameter estimation and model validation were obtained primarily from the published results associated with the compact circadian clock model developed by De Caluwé et al. [9]. Because the original article and its supplementary materials did not provide complete numerical tables for all expression profiles, data points were digitized from the published figures using the PlotDigitizer 3.1.6. Specifically, the time and relative-expression axes were calibrated separately for each figure panel, and the normalized relative mRNA expression levels of individual genes were subsequently extracted at the sampling points shown in the published figures.

Wild-type mRNA expression data under the 12L:12D photoperiod were used as the calibration dataset for parameter estimation. Data obtained under other photoperiods were used for model validation without parameter refitting. The validation results under the 8L:16D and 16L:8D photoperiods are presented in the main text, whereas the results under continuous light (LL) are provided in Appendix B.

Only the mRNA state variables were used as direct fitting targets. The protein state variables and the light-signal-related variable were not included directly in the error calculation; instead, they were treated as internal model variables that influence the dynamics of the core mRNA modules through the governing equations. To ensure comparability between the model outputs and the published expression data, the model outputs were normalized by their maximum values before error calculation. Model predictions at the sampling times of the published data were then obtained by interpolation.

4.2. Parameter Search Space and Prescribed Fractional Order

To preserve the biological scale of the model parameters and prevent unconstrained optimization from producing parameter values that deviate substantially from those of the original model, the parameter search space was defined relative to the original integer-order parameter set reported by De Caluwé et al. [9]. A sufficiently broad search interval was retained to accommodate changes in the dynamical behavior induced by the fractional-order formulation. A total of 34 core clock parameters were included in the parameter-estimation procedure. Let ${\theta }_j^0$ denote the value of the $j$-th parameter in the original De Caluwé model, and let ${\theta }_j$ denote the corresponding parameter in the fractional-order model. The search interval was defined as

$$0.5{\theta }_j^0\le {\theta }_j\le 2{\theta }_j^0,j=\mathrm{1,2},\dots ,34.$$

For parameter estimation, the fractional order was prescribed as $q=0.95$. The same parameter bounds and numerical settings were used throughout the optimization procedure. The composite objective function consisted of an expression-fitting error term, an oscillation-amplitude constraint, and a low-expression-level constraint. The complete mathematical definitions of these components and their corresponding weights are provided in Appendix C.

3.3. Numerical Implementation

All numerical simulations, parameter-estimation procedures, and dynamical analyses were implemented in MATLAB R2024b. With the fractional order fixed at $q=0.95$ during parameter estimation, the fractional-order circadian-clock model was solved using the FDE12 routine, which implements a PECE predictor–corrector scheme based on the fractional Adams–Bashforth–Moulton method for Caputo fractional initial-value problems [31]. The model output was calculated under the 12L:12D photoperiod using a time step of $\mathsf{\Delta }t=0.1$h.

To reduce the influence of the initial conditions on parameter estimation, the model was simulated continuously for 384 h under the 12L:12D photoperiod, and only the final 48 h of the simulation were used for error calculation. Model predictions at the sampling times of the published expression data were obtained by interpolation.

Because both the published expression data and the model outputs were expressed as relative expression levels, the output of each mRNA module was normalized by its maximum value over the final 48 h simulation interval:

$${\widehat{R}}_j\left(t\right)={\displaystyle \frac{R_j\left(t\right)}{\underset{s\in \left[336,384\right]}{\max }{R}_j\left(s\right)}},j=\mathrm{1,2},\mathrm{3,4},$$

where $R_j\left(t\right)$ denotes the simulated expression level of the $j$-th core mRNA module and ${\widehat{R}}_j\left(t\right)$ denotes the corresponding normalized output. The four directly fitted mRNA modules corresponded to CCA1/LHY, PRR9/PRR7, PRR5/TOC1, and ELF4/LUX.

3.4. Parameter Estimation

With $q$ fixed at 0.95, parameter estimation was performed using a combination of random screening and local optimization. Candidate parameter vectors were first generated randomly within the parameter bounds specified in Section 4.2. For each candidate vector, the fractional-order model was numerically simulated under the 12L:12D photoperiod, and the composite objective function defined in Appendix C was evaluated using the normalized mRNA expression profiles at the experimental sampling times.

Candidate parameter vectors with relatively small objective-function values were subsequently used as initial points for local optimization with the Nelder–Mead simplex method [32]. The locally optimized parameter sets were further examined through extended numerical simulations to exclude solutions that provided a low fitting error over the calibration interval but failed to maintain stable rhythmic behavior over longer time scales.

For each prescribed fractional order, the parameter set with a small composite objective-function value and stable long-term oscillatory behavior was retained as the corresponding candidate parameter set.

3.5. Model Evaluation and Mutant Simulation Protocol

To quantitatively compare the fitted FDE model with the original ODE model, the mean squared error (MSE) for each light condition $c$was defined as

$${\mathrm{M}\mathrm{S}\mathrm{E}}_c={\displaystyle \frac{1}{4}}{\displaystyle \sum _{i=1}^4}{\displaystyle \frac{1}{N_{i,c}}}{\displaystyle \sum _{j=1}^{N_{i,c}}}{\left(R_{ij,c}^{\mathrm{s}\mathrm{i}\mathrm{m}}-R_{ij,c}^{\mathrm{e}\mathrm{x}\mathrm{p}}\right)}^2,$$

where $c$denotes the light condition, and $i=\mathrm{1,2},\mathrm{3,4}$ corresponds to the CCA1/LHY, PRR9/PRR7, PRR5/TOC1, and ELF4/LUX modules, respectively. $N_{i,c}$ denotes the number of valid experimental observations for the $i$-th module under condition $c$, and $j$ indexes the observations within that module. The quantities $R_{ij,c}^{\mathrm{s}\mathrm{i}\mathrm{m}}$ and $R_{ij,c}^{\mathrm{e}\mathrm{x}\mathrm{p}}$ denote the corresponding normalized model prediction and experimental expression value, respectively.

Mutant simulations were performed following the procedure used in the original compact model [9]. For each mutant condition, the synthesis-rate parameter associated with the corresponding gene module was reduced to 10% of its wild-type value, while all other model parameters were kept unchanged. This treatment approximated a strong loss-of-function or double-mutant condition rather than complete elimination of gene expression. For mutant-profile comparisons, the simulated and experimental expression levels were normalized relative to the maximum value of the corresponding wild-type expression profile. The MSE was calculated separately for each mutant condition using only the experimental observations associated with that condition.

4. Model Fitting and Validation Results

4.1. Parameter Selection and Model Fitting Under the 12L:12D Photoperiod

Parameter estimation was performed with the fractional order fixed at $q=0.95$, using the parameter bounds, composite objective function, and optimization procedure described in Section 4. The resulting fitted parameter set was retained for model validation and subsequent dynamical analyses. When the effect of the fractional order was investigated in Section 6, all other model parameters were kept at their fitted values, and no additional parameter estimation was performed. The original ODE parameter values, the fitted FDE parameter values, and their biological interpretations are listed in Appendix A.

Under the 12L:12D photoperiod, the core circadian clock components exhibit rhythmic expression entrained to the external 24 h light–dark cycle. The wild-type mRNA expression data obtained under this condition were used as the calibration data for parameter estimation. Using $q=0.95$ and the fitted parameter set, the expression dynamics of the CCA1/LHY, PRR9/PRR7, PRR5/TOC1, and ELF4/LUX composite modules were simulated. Because the model groups genes with similar expression patterns and regulatory functions into composite variables, each model output was compared with the two corresponding experimental gene-expression profiles, as shown in Figure 2.

The FDE model reproduced the principal rhythmic expression patterns of the four core clock modules and largely preserved the temporal ordering of their expression peaks. The CCA1/LHY peak occurred near the onset of light, followed by a rapid decrease and a prolonged tough phase during most of the light period. PRR9/PRR7 subsequently reached its peak during the early light phase, whereas the peaks of PRR5/TOC1 and ELF4/LUX occurred later, primarily during the latter part of the light phase or near the light-to-dark transition. Overall, the model reproduced the characteristic temporal sequence of the dawn-, morning-, and evening-phased modules.

The FDE and original ODE models produced generally similar expression patterns for the CCA1/LHY, PRR9/PRR7, and ELF4/LUX modules. The most pronounced difference was observed for PRR5/TOC1. The original ODE model produced a broader expression peak, a slower post-peak decline, and relatively high simulated expression during the expected tough phase. In contrast, the FDE model shortened the duration of the high-expression state and provided a closer representation of the post-peak decline and tough phase. Thus, under the 12L:12D calibration condition, the FDE model preserved the principal rhythmic characteristics and phase relationships of the core clock modules while reducing selected expression-profile discrepancies.

4.2. Model Validation Under Different Light Conditions

To examine the applicability of the FDE model under photoperiods not used for parameter calibration, the fitted parameter set obtained with $q=0.95$ was applied without parameter re-estimation to simulate the wild-type mRNA expression dynamics under the 8L:16D and 16L:8D photoperiods. The corresponding results are presented in Figure 3 and Figure 4, respectively. The validation results under continuous light are provided in Appendix B.

Under both the 8L:16D and 16L:8D photoperiods, the FDE model generated stable rhythmic behavior entrained to the external 24 h light–dark cycle and preserved the principal expression sequence of the four core clock modules. CCA1/LHY expression remained primarily associated with the onset of light, followed by PRR9/PRR7 expression during the early light phase, whereas the peaks of PRR5/TOC1 and ELF4/LUX occurred later. Changes in photoperiod length were accompanied by corresponding adjustments in the peak positions and durations of elevated expression, indicating that the model responded to changes in the external light–dark schedule.

Both the original ODE model and the FDE model reproduced the principal periodic expression patterns of the core clock modules under the different photoperiods, although differences were observed during specific expression phases. These differences were mainly associated with peak width, the post-peak decline, and the simulated expression level during low-expression intervals. Overall, the FDE model preserved the principal phase relationships among the clock modules while providing closer agreement with the experimental expression data over several time intervals. Nevertheless, because each composite model variable represents two related gene-expression profiles, the degree of agreement was not identical for the two genes within every module.

To further compare the simulation performance of the two models across the calibration and validation conditions, the pointwise mean squared error (MSE) was calculated according to the procedure described in Section 4.5. The MSE was calculated separately for each light condition using all valid wild-type observations under the 8L:16D, 12L:12D, 16L:8D, and continuous-light (LL) conditions to examine whether the difference between the two models was consistently observed across different photoperiods. The condition-specific results are presented in Figure 5.

As shown in Figure 5, the FDE model yielded a lower MSE under each of the four light conditions. Specifically, the MSE decreased from 0.0863 to 0.0513 under 8L:16D, from 0.0687 to 0.0477 under 12L:12D, from 0.0497 to 0.0431 under 16L:8D, and from 0.2031 to 0.1484 under LL. The relative reduction was largest under 8L:16D, at approximately 40.6%, and smallest under 16L:8D, at approximately 13.3%. Both models exhibited their largest MSE under LL, indicating that the continuous-light expression profiles remained the most difficult to reproduce within the present modeling framework. Overall, the fitted FDE model showed lower pointwise discrepancies under all four light regimes, although the magnitude of improvement varied among conditions.

4.3. Validation Under Mutant Conditions

The fitted FDE model was further evaluated using the mutant expression data described in Section 4.5, without additional parameter estimation. Figure 6 compares the simulated expression profiles with the corresponding experimental observations for four major clock-gene mutants.

Compared with the original ODE model, the FDE model reduced several prominent discrepancies in the mutant simulations. In the prr9/prr7 mutant, the excessively high simulated CCA1/LHY peak produced by the ODE model was substantially reduced, resulting in closer agreement with the observed mutant expression level. In the prr5/toc1 mutant, the excessive phase advance of the ODE mutant profile was reduced, and the simulated peak occurred closer to the experimentally observed timing. The FDE model also retained the principal expression changes observed in the cca1/lhy mutant, including an elevated expression level during the tough phase and a reduced oscillation amplitude.

However, the improvement was not consistent across all mutant conditions. In particular, the FDE model did not fully reproduce the experimental expression pattern of the elf4 mutant. Overall, the FDE model improved the description of selected mutant responses, particularly the peak magnitude in the prr9/prr7 mutant and the peak timing in the prr5/toc1 mutant, while mutant-specific discrepancies remained.

To quantify the mutant simulations, MSE values were calculated using only the mutant experimental observations. The condition-specific MSE results are presented in Figure 7.

The largest improvement occurred for the prr9/prr7 mutant, for which the MSE decreased from 1.5507 to 0.3499. The MSE for the cca1/lhy mutant also decreased, from 0.1027 to 0.0689. By contrast, the MSE increased slightly from 0.0077 to 0.0119 for the prr5/toc1 mutant, although the absolute errors of both models remained small. For the elf4 mutant condition, the MSE increased from 0.1359 to 0.1827. Therefore, the FDE model showed a clear advantage in simulating the prr9/prr7 and cca1/lhy mutants but did not exhibit a uniform advantage under every mutant condition.

5. Intrinsic Properties of the Free-Running Solution

Under periodic light conditions, the plant circadian clock can be entrained by the external light–dark cycle. Under continuous darkness (DD) conditions, however, the external light input is absent, and the system is no longer driven by periodic light–dark transitions. This condition therefore provides a suitable setting for investigating the.

Under DD conditions, $L\left(t\right)=0$ and $D\left(t\right)=1$, system (2.1), reduces to the following so-called “free running” system:

$${}_0{}^C{D}_t^q{R}_1=v_1{\displaystyle \frac{1}{1+{\left(\frac{P_2}{K_1}\right)}^2+{\left(\frac{P_3}{K_2}\right)}^2}}-k_{1D}{R}_1,$$

$${}_0{}^C{D}_t^q{P}_1=p_1{R}_1-d_1{P}_1,$$

$${}_0{}^C{D}_t^q{R}_2=\left(v_{2A}+v_{2B}{\displaystyle \frac{P_1^2}{K_3^2+P_1^2}}\right){\displaystyle \frac{1}{1+{\left(\frac{P_3}{K_4}\right)}^2+{\left(\frac{P_4}{K_5}\right)}^2}}-k_2{R}_2,$$

$${}_0{}^C{D}_t^q{P}_2=p_2{R}_2-d_{2D}{P}_2,$$

$${}_0{}^C{D}_t^q{R}_3=v_3{\displaystyle \frac{1}{1+{\left(\frac{P_1}{K_6}\right)}^2+{\left(\frac{P_3}{K_7}\right)}^2}}-k_3{R}_3,(5.1)$$

$${}_0{}^C{D}_t^q{P}_3=p_3{R}_3-d_{3D}{P}_3,$$

$${}_0{}^C{D}_t^q{R}_4=-k_4{R}_4,$$

$${}_0{}^C{D}_t^q{P}_4=p_4{R}_4-d_{4D}{P}_4,$$

$${}_0{}^C{D}_t^q{P}_L=0.3(1-P_L).$$

Equation (5.1) shows that, under DD conditions, the light-induced transcription term in the equation for the Evening Complex-associated mRNA variable $R_4$ vanishes, leaving $R_4$ governed solely by its degradation term. Meanwhile, the light-signaling variable $P_L$ evolves toward the dark-adapted equilibrium $P_L=1$.

This formulation facilitates the analysis of the existence, non-negativity, and boundedness of solutions, as well as the stability of equilibria. The following analysis mainly focuses on the fractional-order case $0<q<1$. When $q=1$, the Caputo fractional derivative reduces to the classical first-order derivative, and the model consequently reduces to an integer-order system of ordinary differential equations.

To ensure the biological meaningfulness of the model, the solutions should remain nonnegative. Their boundedness must also be established. Therefore, we first investigate the local existence and uniqueness, non-negativity, and boundedness of the solutions, and then analyze the stability of the equilibria.

For the subsequent analysis, we denote the state vector of system (5.1) as

$$X\left(t\right)={(R_1\left(t\right),P_2(t),R_2\left(t\right),P_2(t),R_3(t),P_3(t),R_4(t),P_4(t),P_L(t\left)\right)}^{\mathrm{T}},$$

and prescribe the initial condition for system (5.1) as

$$X\left(0\right)=X_0={\left(R_1\right(0),P_1(0),R_2(0),P_2(0),R_3(0),P_3(0),R_4(0),P_4(0),P_L(0\left)\right)}^{\mathrm{T}}\in {\mathbb{R}}_{\ge 0}^9.$$

where

$${\mathbb{R}}_{\ge 0}^9=\left\{x\in {\mathbb{R}}^9:x_i\ge 0,i=1,\dots ,9\right\},$$

3.1. Local Existence and Uniqueness of Solutions

To establish the local existence and uniqueness of solutions to the model under continuous darkness, we first recall the following standard result for Caputo fractional initial-value problems.

Lemma 1([25]). Let $0<q<1$ , let $U\subset {\mathbb{R}}^n$ be an open neighborhood of $X_0$ , and let

$$F:U\to {\mathbb{R}}^n$$

be locally Lipschitz continuous on  $U$  . Then there exists a constant  $T>0$  such that the initial-value problem

$$\left\{\begin{array}{c}{}_0{}^C{D}_t^{q}X\left(t\right)=F\left(X\right(t\left)\right),\\ X\left(0\right)=X_0,\end{array}\right.$$

admits a unique continuous solution  $X\in C\left(\right[0,T];{\mathbb{R}}^n)$.

Theorem 2. 

For any nonnegative initial value, the system (5.1) admits a unique local solution on some time interval  $[0,T].$

Proof. Since all model parameters are positive, every denominator appearing in the Hill-type regulatory functions is strictly positive for all $X\in {\mathbb{R}}^9$. Therefore, the right-hand side of system (5.1) has no singularities in ${\mathbb{R}}^9$.

Each component of $F$ is composed of linear functions and rational functions with nonvanishing denominators. Hence,

$$F\in C^1({\mathbb{R}}^9,{\mathbb{R}}^9).$$

Therefore, $F$ is locally Lipschitz continuous in ${\mathbb{R}}^9$. In particular, for every $X_0\in {\mathbb{R}}_{\ge 0}^9$, there exist a neighborhood $U$ of $X_0$ and a constant $L>0$ such that

$$\parallel F\left(X\right)-F\left(Y\right)\parallel \le L\parallel X-Y\parallel ,$$

for all $X,Y\in U$.

Thus, all assumptions of Lemma 1 are satisfied. Consequently, system (5.1) admits a unique local solution

$$X\in C\left([0,T];{\mathbb{R}}^9\right)$$

for some $T>0$. This completes the proof.

Remark 1. 

Theorem 2 establishes only the local existence and uniqueness of solutions. Global existence will be established after proving the non-negativity and boundedness of solutions in the following sections. By the continuation theory for Caputo fractional differential equations [26], the local solution obtained in Theorem 2 admits a unique maximal continuation.

Let

$$T_{max}=\mathit{sup}\left\{T>0:\mathrm{system}\left(3.1\right)\mathrm{admits}\mathrm{a}\mathrm{unique}\mathrm{solution}\mathrm{on}[0,T]\right\}.$$

Then  $[0,T_{max})$  is the maximal interval Since Theorem 2 guarantees the existence of a unique local solution on some interval  $[0,T]$  with  $T>0$  , we have

$$0<T\le T_{max}\le +\infty .$$

In the following analysis, let

$$\tau \in (0,T_{max})$$

be arbitrary, and consider the solution on the interval $[0,\tau ]$.

3.2. Non-Negativity of Solutions

Lemma 3([27]). Let $0<q<1$ , $u\in AC\left(\right[0,T\left]\right)$ , and suppose that

$${}_0{}^C{D}_t^{q}u\left(t\right)$$

exists. If there exists $\left.t_0\in (0,T\right]$ such that $u\left(t\right)$ attains its minimum at $t_0$ on the interval $[0,t_0]$ , that is,

$$u\left(t\right)\ge u\left(t_0\right),0\le t\le t_0,$$

then

$${}_0{}^C{D}_t^{q}u\left(t_0\right)\le 0.$$

In particular, if

$$u\left(t\right)\ge \mathrm{0,0}\le t\le t_0,u\left(t_0\right)=0,$$

then

$${}_0{}^C{D}_t^{q}u\left(t_0\right)\le 0.$$

Lemma 4([23,28]). Let $0<q<1$ and $m\in \mathbb{N}$ . Denote

$$X\left(t\right)=(X_1\left(t\right),\dots ,X_m\left(t\right){)}^T,f={\left(f_1,\dots ,f_m\right)}^T.$$

For each  $i=1,\dots ,m$  , suppose that

$$f_i:[0,T]\times {\mathbb{R}}^m\to \mathbb{R}$$

is continuous and satisfies the local Lipschitz condition. Consider the following two systems of Caputo fractional differential equations:

$$\left\{\begin{array}{c}{}_0{}^C{D}_t^q{X}^{\left(k\right)}\left(t\right)=f(t,X^{\left(k\right)}(t\left)\right)+{\displaystyle \frac{1}{k}}\overrightarrow{1},\\ X^{\left(k\right)}\left(0\right)=X_0+{\displaystyle \frac{1}{k}}\overrightarrow{1},\end{array}\right.$$

and

$$\left\{\begin{array}{c}{}_0{}^C{D}_t^{q}X\left(t\right)=f(t,X(t\left)\right),\\ X\left(0\right)=X_0,\end{array}\right.$$

where $k$ is a positive integer and $\overrightarrow{1}$ is the m-dimensional vector where each entry takes value 1. Then, their solutions satisfy

$$X^{\left(k\right)}\left(t\right)\to X\left(t\right),k\to \infty ,\forall t\in [0,T].$$

Theorem 5. 

Let  $X\left(t\right)$ be the solution of system (5.1) on  $[0,\tau ]$  , where  $0<\tau <T_{max}$  . For every nonnegative initial value,

$$X\left(t\right)\in {\mathbb{R}}_{\ge 0}^9,t\in [0,\tau ].$$

Proof. To prove the non-negativity of the solution, consider the following auxiliary system:

$${}_0{}^C{D}_t^q{X}^{\left(k\right)}\left(t\right)=F\left(X^{\left(k\right)}\right(t\left)\right)+{\displaystyle \frac{1}{k}}\overrightarrow{1},$$

$$X^{\left(k\right)}\left(0\right)=X_0+{\displaystyle \frac{1}{k}}\overrightarrow{1},$$

where

$$\overrightarrow{1}=(\mathrm{1,1},\dots ,1{)}^{\mathrm{T}}\in {\mathbb{R}}^9,k\in \mathbb{N}.$$

We now prove that the solution of the auxiliary system satisfies

$$X^{\left(k\right)}\left(t\right)\in {\mathbb{R}}_{>0}^9.$$

Assume that the conclusion is not true, which means that here exists a first time $t_0>0$ at which the solution touches the boundary of the nonnegative region, such that

$$X^{\left(k\right)}\left(t\right)\in {\mathbb{R}}_{>0}^9,0\le t<t_0,$$

and at $t=t_0$, at least one component vanishes. Without loss of generality, we suppose that the first component satisfies

$$R_1^{\left(k\right)}\left(t_0\right)=0.$$

That is,

$$t_0:=\mathrm{i}\mathrm{n}\mathrm{f}\left\{t>0:X^{\left(k\right)}\left(t\right)\notin {\mathbb{R}}_{>0}^9\right\}.$$

Therefore, $R_1^{\left(k\right)}\left(t\right)$ is nonnegative on the interval $[0,t_0]$ and attains its minimum at $t_0$. Moreover,

$$R_1^{\left(k\right)}\left(0\right)=R_1\left(0\right)+{\displaystyle \frac{1}{k}}>0.$$

By Lemma 3, we obtain

$${}_0{}^C{D}_t^q{R}_1^{\left(k\right)}\left(t_0\right)\le 0.$$

On the other hand, according to the specific structure of system (5.1), we have

$${}_0{}^C{D}_t^q{R}_1^{\left(k\right)}\left(t_0\right)=v_1{\displaystyle \frac{1}{1+{\left(\frac{P_2\left(t_0\right)}{K_1}\right)}^2+{\left(\frac{P_3\left(t_0\right)}{K_2}\right)}^2}}-k_{1D}{R}_1^{\left(k\right)}\left(t_0\right)+{\displaystyle \frac{1}{k}}$$

$$=v_1{\displaystyle \frac{1}{1+{\left(\frac{P_2\left(t_0\right)}{K_1}\right)}^2+{\left(\frac{P_3\left(t_0\right)}{K_2}\right)}^2}}+{\displaystyle \frac{1}{k}}>0$$

This contradicts

$${}_0{}^C{D}_t^q{R}_1^{\left(k\right)}\left(t_0\right)\le 0.$$

The remaining components can be treated analogously.

Hence, for any $t\in [0,\tau ]$,

$$X^{\left(k\right)}\left(t\right)\in {\mathbb{R}}_{>0}^9.$$

By Lemma 4,

$$X^{\left(k\right)}\left(t\right)\to X\left(t\right)(k\to \mathrm{\infty }).$$

Therefore,

$$X\left(t\right)\in {\mathbb{R}}_{\ge 0}^9$$

holds for all $t\in [0,\tau ]$. Thus, the solution of system (5.1) always remains nonnegative. $\square$

3.3. Boundedness of Solutions

Lemma 6. 

Let  $0<q<1$. Suppose that  $u\left(t\right)$  is continuous on  $[0,T]$  and that its Caputo fractional derivative  ${}_0{}^C{D}_t^{q}u\left(t\right)$  exists. If there exist constants  $a\ge 0$  and  $b>0$  such that

$${}_0{}^C{D}_t^{q}u\left(t\right)\le a-bu\left(t\right),$$

then, by the fractional-order comparison principle [29], we have

$$u\left(t\right)\le u\left(0\right)E_q(-b{t}^q)+{\displaystyle \frac{a}{b}}\left[1-E_q(-b{t}^q)\right].$$

Here,  $E_q(\cdot )$  denotes the Mittag–Leffler function. Since, for  $0<q<1$,

$$0<E_q(-b{t}^q)\le 1,$$

it follows that

$$u\left(t\right)\le \mathit{max}\left\{u\left(0\right),{\displaystyle \frac{a}{b}}\right\}.$$

Proof. Let $v$ be the solution of

$${}_0{}^C{D}_t^{q}v\left(t\right)=a-bv\left(t\right),v\left(0\right)=u\left(0\right).$$

By the comparison principle for Caputo fractional differential equations ,

$$u\left(t\right)\le v\left(t\right).$$

The solution of the linear problem is [23]

$$v\left(t\right)=u\left(0\right)E_q\left(-b{t}^q\right)+{\displaystyle \frac{a}{b}}\left[1-E_q\left(-b{t}^q\right)\right].$$

Since $E_q\left(-x\right)$ is completely monotone for $0<q<1$[30],

$$0<E_q\left(-b{t}^q\right)\le 1.$$

Therefore,

$$u\left(t\right)\le \max \left\{u\right(0),{\displaystyle \frac{a}{b}}\}.\square$$

Theorem 7. 

Let  $X\left(t\right)$ be the solution of system (5.1) on $[0,\tau ]$ , where $0<\tau <T_{max}$ . Then there exists a constant $M>0$ , independent of $\tau$ , such that

$$0\le R_i\left(t\right),P_i\left(t\right),P_L\left(t\right)\le M,t\in [0,\tau ],i=1,\dots ,4.$$

Proof. First, since

$$0\le {\displaystyle \frac{1}{1+{\left(\frac{P_2}{K_1}\right)}^2+{\left(\frac{P_3}{K_2}\right)}^2}}\le 1,$$

we have

$${}_0{}^C{D}_t^q{R}_1\le v_1-k_{1D}{R}_1.$$

By the comparison principle, we obtain

$$R_1\left(t\right)\le M_{R_1},$$

where

$$M_{R_1}=\mathrm{m}\mathrm{a}\mathrm{x}\left\{R_1\left(0\right),{\displaystyle \frac{v_1}{k_{1D}}}\right\}.$$

Then,

$${}_0{}^C{D}_t^q{P}_1\le p_1{M}_{R_1}-d_1{P}_1.$$

By the comparison principle, we obtain

$$P_1\left(t\right)\le M_{P_1},$$

where

$$M_{P_1}=\mathrm{m}\mathrm{a}\mathrm{x}\left\{P_1\left(0\right),{\displaystyle \frac{p_1{M}_{R_1}}{d_1}}\right\}.$$

Since the equations for the other six variables have the same structure as those for $R_1$ and $P_1$, namely, their right-hand sides can all be estimated in the form of a bounded nonnegative input term minus a linear degradation term, their boundedness can be proved in exactly the same manner by using the fractional-order comparison principle.

define

$$M_{R_2}=\mathit{max}\left\{R_2\right(0),{\displaystyle \frac{v_{2A}+v_{2B}}{k_2}}\},$$

$$M_{P_2}=\mathit{max}\{P_2\left(0\right),{\displaystyle \frac{p_2{M}_{R_2}}{d_{2D}}}\},$$

$$M_{R_3}=\mathit{max}\{R_3\left(0\right),{\displaystyle \frac{v_3}{k_3}}\},$$

$$M_{P_3}=\mathit{max}\{P_3\left(0\right),{\displaystyle \frac{p_3{M}_{R_3}}{d_{3D}}}\},$$

$$M_{R_4}=R_4\left(0\right),$$

$$M_{P_4}=\mathit{max}\{P_4\left(0\right),{\displaystyle \frac{p_4{M}_{R_4}}{d_{4D}}}\}.$$

Finally, for the light-signal-related variable $PL$, since

$${}_0{}^C{D}_t^q{P}_L=0.3(1-P_L)=0.3-0.3P_L,$$

we obtain

$$P_L\left(t\right)\le M_{P_L},$$

where

$$M_{P_L}=\mathrm{m}\mathrm{a}\mathrm{x}\left\{P_L\right(0),1\}.$$

Let

$$M=\mathrm{m}\mathrm{a}\mathrm{x}\left\{M_{R_1},M_{P_1},\dots ,M_{R_4},M_{P_4},M_{P_L}\right\}.$$

Then,

$$0\le R_i\left(t\right),P_i\left(t\right),P_L\left(t\right)\le M,t\in [0,\tau ],i=1,\dots ,4.$$

Therefore, the solution of system (5.1) is bounded over its interval of existence. This completes the proof. $\square$

Corollary 8 For every$X_0\in {\mathbb{R}}_{\ge 0}^9$ , system (5.1) admits a unique global nonnegative and bounded solution.

Proof. By Theorems 2,5,7, the unique local solution remains nonnegative and uniformly bounded throughout its maximal interval of existence. The continuation theorem for Caputo fractional differential equations [26] therefore implies that $T_{\mathrm{m}\mathrm{a}\mathrm{x}}=\infty$. Hence, the solution is global. $\square$

6. Long-Term Behavior of the Free-Running System

6.1. Equilibrium Point and Local Stability Under Continuous Darkness

Under DD continuous when the external light input is absent, and the system assumes $L=0$ and $D=1$. Under this condition, the system is no longer driven by periodic light–dark inputs, and its dynamics are governed primarily by the intrinsic feedback interactions among the core clock genes. Therefore, the system under continuous darkness is suitable for investigating the intrinsic dynamical properties of the model.

Although the equilibrium values of the light-related variable $P_L$ and the ELF4/LUX-related variables $R_4$ and $P_4$ can be determined directly from the structure of the equations under DD conditions, the remaining core variables $R_1,P_1,R_2,P_2,R_3,$and $P_3$ remain coupled through multiple Hill-type activation and repression terms. The equilibrium equations obtained by setting the right-hand sides of the system to zero contain rational nonlinear terms and regulatory terms with a Hill coefficient of 2. Consequently, it is difficult to obtain explicit analytical expressions for the equilibrium point in terms of the model parameters. Therefore, the equilibrium point of the fractional-order model under DD conditions was determined numerically. The Jacobian matrix and its eigenvalues were then calculated at the equilibrium point to assess its local stability.

Theorem 9([33]). Consider the Caputo fractional-order linear system

$${}_0{}^C{D}_t^{q}X\left(t\right)=AX\left(t\right),0<q<1,$$

where  $A\in {\mathbb{R}}^{n\times n}$  , and let its eigenvalues be denoted by  ${\lambda }_i$  ,  $i=\mathrm{1,2},\dots ,n$  . If all eigenvalues satisfy

$$∣arg\left({\lambda }_i\right)∣>{\displaystyle \frac{q\pi }{2}},i=\mathrm{1,2},\dots ,n,$$

then the linear fractional-order system is asymptotically stable. If there exists at least one eigenvalue satisfying

$$∣arg\left({\lambda }_i\right)∣\le {\displaystyle \frac{q\pi }{2}},$$

then the system does not satisfy the asymptotic stability condition.

Remark 2. Theorem 9 provides a stability criterion for Caputo fractional-order linear systems. For the nonlinear autonomous Caputo fractional-order system

$${}_0{}^C{D}_t^{q}X=F\left(X\right),0<q<1,$$

suppose that  $E$  is an equilibrium point of the system and that  $F\left(X\right)$  is continuously differentiable in a neighborhood of  $E$  . The system can then be linearized at the equilibrium point as

$${}_0{}^C{D}_t^q\xi =J\left(E\right)\xi ,$$

where  $J\left(E\right)$  denotes the Jacobian matrix evaluated at  $E$  . According to the linearization theory for nonlinear Caputo fractional-order systems, if all eigenvalues of  $J\left(E\right)$  satisfy

$$∣arg\left({\lambda }_i\right)∣>{\displaystyle \frac{q\pi }{2}},$$

then the equilibrium point $E$ is locally asymptotically stable [34]. If at least one eigenvalue satisfies

$$∣\mathit{arg}\left({\lambda }_i\right)∣<{\displaystyle \frac{q\pi }{2}},$$

then the equilibrium point $E$ is unstable [35].

Therefore, the local stability of the equilibrium point is determined by calculating the eigenvalues of the Jacobian matrix at the equilibrium point and examining whether they satisfy the above argument conditions.

Under DD conditions, the equation for the light-related variable $P_L$ reduces to

$${}_0{}^C{D}_t^q{P}_L=0.3\left(1-P_L\right).$$

The right-hand side of this equation is independent of the other state variables, and $P_L$ does not appear in the equations for the other core clock variables. Therefore, the equation for $P_L$ is decoupled from the core circadian subsystem. Setting the right-hand side of the above equation to zero yields

$${P_L}^{\mathrm{*}}=1.$$

The corresponding eigenvalue of the linearized equation is $-0.3$, whose argument is $\pi$, and hence it satisfies the local asymptotic stability condition for fractional-order systems. Therefore, the light-related variable $P_L$ does not affect the stability determination of the system equilibrium point under DD conditions.

Furthermore, under DD conditions, the equation for the ELF4/LUX mRNA variable $R_4$ is given by

$${}_0{}^C{D}_t^q{R}_4=-k_4{R}_4,$$

and its equilibrium value is

$$R_4^{*}=0.$$

From the equation for the ELF4/LUX protein variable $P_4$,

$${}_0{}^C{D}_t^q{P}_4=p_4{R}_4-d_{4D}{P}_4,$$

it follows that $P_4^{*}=0$ when $R_4^{*}=0$. Therefore, under DD conditions, the last three components of the system equilibrium point are

$$R_4^{*}=0,P_4^{*}=0,{P_L}^{*}=1.$$

Using the parameter values fitted for the present model, the equilibrium equations under DD conditions were solved numerically, yielding the following system equilibrium point:

$$E^{*}=(\mathrm{1.99201,1.38860,0.532908,0.773259,0.192181,0.610521,0},\mathrm{0,1}).$$

To determine the local stability of the equilibrium point $E^{*}$, the Jacobian matrix of the right-hand side of system (5.1) was evaluated at $E^{*}$. The resulting eigenvalues were

$${\lambda }_1=-1.17804,$$

$${\lambda }_{\mathrm{2,3}}=-1.10361\pm 0.262304i,$$

$${\lambda }_4=-0.760893,$$

$${\lambda }_{\mathrm{5,6}}=-0.515225\pm 0.541352i,$$

$${\lambda }_7=-0.3,$$

$${\lambda }_{\mathrm{8,9}}=0.0509705\pm 0.274437i.$$

The fractional order considered in this study is $q=0.95$; therefore, the critical angle in the stability criterion is

$${\displaystyle \frac{q\pi }{2}}={\displaystyle \frac{0.95\pi }{2}}=1.49226.$$

The eigenvalue calculations show that the Jacobian matrix has a pair of complex-conjugate eigenvalues,

$${\lambda }_{\mathrm{8,9}}=0.0509705\pm 0.274437i,$$

whose arguments satisfy

$$∣\mathrm{a}\mathrm{r}\mathrm{g}\left({\lambda }_{\mathrm{8,9}}\right)∣=1.38716.$$

Thus,

$$∣\mathrm{a}\mathrm{r}\mathrm{g}\left({\lambda }_{\mathrm{8,9}}\right)∣<{\displaystyle \frac{q\pi }{2}}.$$

Therefore, the equilibrium point does not satisfy the local asymptotic stability condition for fractional-order systems. It follows that, for the present model and the corresponding parameter set, the equilibrium point $E_D$ under DD conditions is locally unstable.

Further numerical simulations showed that, for the given positive initial condition, the core clock variables $R_1,P_1,R_2,P_2,R_3$and $P_3$ exhibited persistent bounded oscillations after a short transient and did not converge to the equilibrium point $E_D$. In contrast, $R_4$ and $P_4$ gradually decayed to zero, while the light-related variable $P_L$ approached its equilibrium value ${P_L}^{*}=1$. These results are consistent with the qualitative analysis under DD conditions and the preceding stability determination, indicating that the present model can generate intrinsic oscillations in the absence of periodic external light–dark inputs. The corresponding simulation results are shown in Figure 8.

6.2. Hopf Bifurcation for the Fractional Order $q$ Under DD Conditions

Based on the eigenvalue argument criterion applied in the previous section to the Jacobian matrix at the equilibrium $E_D$ under DD, the critical fractional order of the present model was determined as

$$q_c={\displaystyle \frac{2}{\pi }}\underset{1\le i\le 9}{\mathrm{m}\mathrm{i}\mathrm{n}}∣\mathrm{a}\mathrm{r}\mathrm{g}\left({\lambda }_i\right)∣=0.8831.$$

Therefore, when $q<q_c$, the equilibrium $E_D$ satisfies the local asymptotic stability condition for fractional-order systems. When $q>q_c$, this stability condition is no longer satisfied, and the equilibrium loses its local asymptotic stability. It should be noted that violation of the stability condition only indicates the loss of stability of the equilibrium and does not directly determine the long-term dynamical behavior of the system after destabilization. Further analysis based on numerical simulations is therefore required.

To investigate the effect of the fractional order $q$ on the long-term dynamical behavior of the system under DD, $q$ was scanned over the interval $q\in \left[\mathrm{0.70,1.00}\right]$ with a step size of 0.001. For the region $q<q_c$, the stable equilibrium branch corresponding to the $R_1$ component of $E_D$ was plotted. For the region $q>q_c$, long-term numerical simulations were performed from the same nonequilibrium initial condition, and the maximum and minimum values of $R_1$ over the interval $t\in \left[\mathrm{5000,6000}\right]$ were extracted. Meanwhile, the unstable equilibrium branch was retained and represented by a dashed line.

Numerical simulation results show that, when $q<q_c$, the equilibrium remains locally asymptotically stable, whereas when $q>q_c$, the equilibrium loses stability, which is consistent with the theoretical analysis presented above. For $q>q_c$, the maximum and minimum values of $R_1$ obtained from long-term simulations lie on opposite sides of the equilibrium value, forming separated upper and lower envelopes. As $q$ increases further, the distance between the upper and lower envelopes becomes larger, indicating that the amplitude of the long-term oscillations increases. These results show that variation in the fractional order $q$ can significantly affect both the stability of the equilibrium under DD and the long-term oscillatory behavior after the loss of stability. The specific results are shown in Figure 9.

To further verify the stability results described above and the differences in dynamical behavior shown in Figure 9, two representative fractional orders on either side of the critical value, $q=0.85$ and $q=0.91$, were selected. The time series of $R_1$ and the phase trajectories in the $(R_1,P_1)$ plane were then compared.

The numerical simulation results show that, when $q=0.85<q_c$, $R_1$ exhibits damped oscillations during the initial stage. The oscillation amplitude gradually decreases over time, and $R_1$ eventually converges to the equilibrium. The corresponding phase trajectory in the $(R_1,P_1)$ plane forms an inward-contracting spiral and ultimately converges to the equilibrium. When $q=0.91>q_c$, $R_1$ continues to exhibit pronounced non-decaying oscillations after a short transient period. The corresponding phase trajectory no longer contracts toward the equilibrium but instead forms a bounded, approximately closed orbit around it. These numerical results are consistent with the eigenvalue argument stability criterion and the changes in long-term dynamical behavior shown in Figure 9. The specific results are presented in Figure 10.

It should be noted that previous studies have shown that autonomous Caputo fractional-order systems do not admit nontrivial periodic solutions [36,37]. Accordingly, the observed dynamical behavior is not interpreted as a limit cycle in the sense of ordinary differential equations but is instead described as long-term bounded oscillations.

6.3. Effect of the Fractional Order on the Phases of Core Clock Variables Under Different Photoperiods

The previous section examined the effect of the fractional order $q$ on the long-term dynamical behavior of the system under DD conditions and showed that $q$ can significantly affect system stability and oscillation amplitude. Under periodic light–dark inputs, however, the circadian system must not only maintain stable oscillations but also preserve appropriate temporal relationships between the expression phases of the core clock modules and the external photoperiod. Therefore, to further characterize the regulatory effect of the fractional order on system dynamics, this section investigates the effect of $q$ on the phases of the core state variables under different photoperiods.

Three photoperiods, namely, 12L:12D, 16L:8D, and 8L:16D, were considered. The fractional order was scanned over the interval $q\in \left[\mathrm{0.8,1}\right]$, while all other model parameters were kept unchanged. For each combination of $q$ and photoperiod, the model was simulated continuously for 384 h from the same initial condition. The phases of $R_1$, $R_2$, $R_3$, and $R_4$ during the final five complete light–dark cycles were then extracted, and their mean phases were calculated. The onset of light on each day was defined as ZT0. It should be noted that ZT is defined over a complete 24 h light–dark cycle; therefore, ZT24 is equivalent to ZT0 of the subsequent cycle. The calculated results are shown in Figure 11.

As shown in Figure 11, the phase responses of the different state variables to $q$ differ markedly. Under all three photoperiods, the phase of $R_1$ remains close to ZT0, indicating that it is relatively insensitive to changes in $q$. The peak of $R_2$ occurs mainly between ZT4 and ZT6, and its phase exhibits only small fluctuations as $q$ varies. In contrast, $R_3$ and $R_4$ are more sensitive to changes in $q$. Under the 12L:12D and 16L:8D photoperiods, their phases are delayed as $q$ decreases. Under the 8L:16D photoperiod, the phase variations of $R_3$ and $R_4$ are relatively small, and the phase of $R_4$ remains approximately stable near ZT8.

Overall, the effect of the fractional order on the phases of the core state variables does not exhibit a monotonic pattern but instead shows clear variable dependence and photoperiod dependence. Together with the results of the previous section, these findings indicate that, when all other parameters are kept unchanged, the fractional order $q$ not only affects equilibrium stability and oscillation amplitude under continuous darkness but also modulates the phases of selected core variables under periodic light–dark inputs. Therefore, $q$ can serve as a dynamical regulatory parameter that influences the phases of selected core clock modules. It should be noted, however, that the specific biological interpretation of this phase-modulation effect requires further validation using experimental data.

6.4. Stability Transitions Induced by the CL Protein Degradation Rate $d_1$

Protein degradation controls protein lifetime and accumulation and thereby influences the strength and timing of circadian feedback regulation. In the present model, $d_1$ denotes the degradation rate of the CCA1/LHY protein. Because this protein regulates several other core clock modules, changes in $d_1$ may alter both its steady-state level and the overall stability of the feedback network. Therefore, $d_1$ was selected as a representative parameter for examining degradation-rate-dependent stability transitions.

To investigate its effect, the fractional order $q$ and all other parameters were fixed under DD conditions, while $d_1$ was varied over the prescribed interval. For each value of $d_1$, the equilibrium point was calculated numerically, and its local stability was determined from the Jacobian eigenvalues using the Matignon criterion. Stable and unstable equilibrium branches were then identified. In the unstable region, long-term numerical simulations were performed, and the extrema of the CL mRNA variable $R_1$ were extracted after removal of the transient response.

As shown in Figure 12, two critical values, denoted by $d_1^{\left(1\right)}$ and $d_1^{\left(2\right)}$, were identified. When $d_1<d_1^{\left(1\right)},$ the equilibrium is locally asymptotically stable, and the state variables converge to it after damped transient oscillations. When $d_1^{\left(1\right)}<d_1<d_1^{\left(2\right)},$ the equilibrium becomes locally unstable, and the maximum and minimum values of $R_1$ separate from the equilibrium branch, indicating sustained bounded oscillatory behavior. When $d_1>d_1^{\left(2\right)},$ the extrema merge again with the equilibrium branch, and the system returns to a locally stable equilibrium.

Related parameter-dependent oscillatory behavior has also been investigated in integer-order biological feedback models. Protein degradation rates have been shown to influence the period of circadian oscillations in minimal clock models [38], while relative degradation rates and feedback-loop architecture can modify the conditions for oscillations in gene-regulatory networks [39]. Parameter-dependent dynamical transitions have also been investigated through bifurcation analysis in an interlocked mammalian circadian-clock model [40]. Although these studies involve model structures and stability mechanisms different from those of the present Caputo fractional-order system, they support the broader conclusion that degradation processes and feedback architecture can shape oscillatory regimes in biological networks.

Overall, the influence of $d_1$ is nonmonotonic: both relatively low and relatively high degradation rates produce locally stable equilibria, whereas an intermediate interval leads to equilibrium instability and long-term bounded oscillations. These results indicate that the CCA1/LHY protein degradation rate affects not only its equilibrium abundance but also the stability regime of the fractional-order core circadian clock.

6.5. Hopf Bifurcation for the Fractional Order $q$ and the CL Protein Degradation Rate $d_1$ Under DD Conditions

To investigate the combined effects of the fractional order $q$ and the CL protein degradation rate $d_1$, the local stability regions in the $(q,d_1)$ parameter plane were determined under DD conditions, while all other model parameters were fixed. For each value of $d_1$, the equilibrium point and the eigenvalues of the corresponding Jacobian matrix were calculated. Because the equilibrium equations are independent of the fractional order, variation in $q$ affects the stability determination through the critical angle ${\displaystyle \frac{q\pi }{2}}$ in the Matignon criterion. The resulting stability regions are shown in Figure 13A.

When $q$ is relatively small, the equilibrium remains locally asymptotically stable throughout the investigated range of $d_1$. Once $q$ exceeds approximately $0.86$, an unstable region emerges in the parameter plane. For a fixed value of $q$ above this threshold, increasing $d_1$ may lead to successive transitions from stability to instability and then back to stability. The unstable interval is bounded by two critical degradation rates, and its width generally increases with $q$. This behavior can be explained by the increase in the critical angle ${\displaystyle \frac{q\pi }{2}}$, which makes the Matignon stability condition more restrictive as $q$ increases.

Representative parameter combinations from the stable and unstable regions were selected for numerical verification. As shown in Figure 13B, when $q=0.95$ and $d_1=4.5$, the variable $R_1$ converges to the equilibrium after damped transient oscillations. By contrast, when $q=0.95$ and $d_1=1.0$, the equilibrium is locally unstable, and the system exhibits sustained bounded oscillations, as shown in Figure 13C.

Overall, $d_1$ and $q$ influence equilibrium stability through different mechanisms. Variation in $d_1$ changes the equilibrium point and the eigenvalues of the Jacobian matrix, whereas variation in $q$ changes the critical angle in the fractional-order stability criterion. Their combined effects determine the stability boundary and the extent of the unstable region in the $\left(q,d_1\right)$ parameter plane.

6.6. Hopf Bifurcation for the CL Protein Degradation Rate $d_1$ and the P51 Protein Degradation Rate $d_{3D}$ Under DD Conditions

To investigate the combined effects of protein degradation in different core clock modules, the fractional order was fixed at $q=0.95$ under DD conditions, while the CL protein degradation rate $d_1$ and the P51 protein degradation rate $d_{3D}$ were varied jointly. All other model parameters were kept unchanged. For each parameter combination, the equilibrium point and the eigenvalues of the corresponding Jacobian matrix were calculated, and local stability was classified using the Matignon criterion. The resulting stability regions are shown in Figure 14A.

Within the investigated parameter ranges, the unstable region is enclosed by an approximately closed stability boundary. Parameter combinations outside this boundary correspond to locally asymptotically stable equilibria, whereas those inside the boundary correspond to unstable equilibria. Thus, equilibrium stability is determined by the combined values of $d_1$ and $d_{3D}$, rather than by either degradation rate alone. For a fixed value of one parameter, variation in the other may lead to successive transitions from stability to instability and then back to stability if the corresponding parameter path crosses the unstable region.

Representative parameter combinations were selected for numerical verification. As shown in Figure 14B, when $d_1=3.2$ and $d_{3D}=2.6$, the variable $R_1$ undergoes damped transient oscillations and eventually converges to the equilibrium, consistent with local asymptotic stability. By contrast, when $d_1=0.8$ and $d_{3D}=0.5$, the equilibrium is locally unstable, and $R_1$ exhibits sustained bounded oscillations, as shown in Figure 14C.

Overall, the two protein degradation rates interact to determine the stability region of the fractional-order core circadian clock. The approximately closed unstable region indicates that system stability depends on coordinated changes in degradation processes across different clock modules. Therefore, analysis of a single degradation rate alone may not fully characterize the parameter-dependent stability of the system.

7. Discussion

In this study, we developed a Caputo fractional-order model of the Arabidopsis core circadian clock based on the established compact model. The original regulatory structure was retained, and the fractional order $q$ was introduced to investigate its influence on system dynamics. We established the local existence and uniqueness, non-negativity, and boundedness of the model solutions. The model parameters were then estimated using wild-type mRNA expression data under 12L:12D, and the fitted model was further validated under other photoperiods and in several clock-gene mutants.

Compared with the original ODE model, the FDE model better described the post-peak decline and low-expression stage of the PRR5/TOC1 module. Under the evaluation method adopted in this study, the FDE model yielded lower MSE values under all four wild-type light conditions and under two of the four mutant conditions. In addition, without increasing the number of core state variables, the model reproduced the main expression order and rhythmic patterns of the four core clock modules under different photoperiods. These results indicate that the FDE model can better describe the temporal expression of core circadian clock genes.

The dynamical analysis showed that the fractional order affects not only the fitting results but also the behavior of the system. Under DD, the fractional order $q$ affected the stability of the equilibrium. Changes in $q$ also affected the oscillation amplitude and the phases of some core variables. In addition, the joint analyses of $q$, $d_1$, and $d_{3D}$ showed that protein degradation rates and the fractional order together determine the stable and unstable regions. These results are broadly consistent with previous studies showing that protein degradation rates can influence circadian periods [38], that relative degradation rates and feedback-loop architecture affect the conditions for oscillations in gene-regulatory networks [39], and that interlocked feedback loops shape parameter-dependent dynamics in mammalian circadian-clock models [40].

This study still has several limitations. The parameter fitting and model validation were mainly based on mRNA data extracted from published figures, and suitable protein time-series data were not available. The fractional order is used here as an effective mathematical parameter representing the combined historical effects of multiple regulatory processes, rather than as a quantity uniquely associated with a specific molecular mechanism. In addition, the improvement in mutant fitting was not uniform across all mutant conditions. Future studies could use more detailed mRNA and protein time-series data, include additional mutant and free-running experiments, investigate the use of different fractional orders for different core clock modules, and compare the present fractional-order formulation with alternative descriptions of historical dependence, such as time-delay models.

8. Conclusions

In this study, a Caputo fractional-order model of the Arabidopsis core circadian clock was developed based on the established compact model. The local existence and uniqueness, non-negativity, and boundedness of the model solutions were established. The model parameters were estimated using wild-type mRNA expression data under 12L:12D, and the fitted model was then validated under other photoperiods and in several clock-gene mutants. Compared with the original ODE model, the FDE model better described the post-peak decline and low-expression stage of the evening-phased PRR5/TOC1 module and yielded lower MSE values under all four wild-type light conditions and under two of the four mutant conditions.

The model results showed that the order $q$ affects equilibrium stability, oscillation amplitude, and the phases of some core clock variables. Under DD, changes in $q$ altered the local stability of the equilibrium, and long-term bounded oscillations were observed when the equilibrium became unstable. The joint stability analyses further showed that the order and protein degradation rates together determine the stable and unstable regions of the system. Overall, the FDE model can better describe the temporal expression and dynamical behavior of the Arabidopsis core circadian clock and provides a basis for further studying the effects of fractional order and protein degradation on circadian regulation.