Synchronization of the Circadian Rhythms with Memory: A Simple Fractional-Order Dynamical Model Based on Tow Coupled Oscillators

ABSTRACT Disruptions of the circadian rhythm are associated with internal desynchronization. It affects some internal functions of our body and behavior that are important to our health. Our modern lifestyle has contributed to millions of people developing some circadian rhythm disruptions, making the subject very important clinically as well as economically. Motivated by studying a simple mathematical model that can reveal some features of internal synchronization/desynchronization, in this contribution, we extend the coupling oscillator phase model proposed by Strogatz14 in the sense that memory is considered in the modeling. Such memory is a result of the introduction of Caputo-type fractional derivatives in the coupling oscillators’ phase model dynamics, resulting in a fractional phase model. We show that the proposed fractional coupling oscillator phase model is well-posed. Furthermore, we analyze the synchronization phenomena. We obtain the synchronized solutions explicitly when the memory is equally distributed between the oscillators. In contrast, when distinct levels of memory are imposed in the modeling, we obtain lower and upper bounds because any existing synchronized solution must be confined in between. We present numerical realizations that support the theoretical findings in great detail.


INTRODUCTION
Synchronization of oscillators have profound implications in many natural and technological phenomena, ranging for neuroscience, where the contributions are in the development of demand-controlled brain pacemakers for the therapy of neurological and psychiatric diseases, cardiac pacemaker e.g [11,16], circadian rhythms e.g.[15,22], chemistry e.g [9], to laser physics and electronics e.g [11,15,23], among others.
From the mathematical point of view, synchronization of coupled oscillators has been studied since the seminal works of Winfree [22] and Kuramoto [9].The general picture for the synchronization of coupled oscillators depends on the topology of the coupling and the number of oscillators.However, there are still many unsolved problems on this subject.[1,2,12,15] and references provide a comprehensive overview of the development and results of this setting.
Typically, circadian and sleep-wake rhythms are phase locked (synchronously) within and with the light-darkness environment.However, the modern human lifestyle often disrupts the usual schedule as a result of rotating shift work, jet lag, insomnia, or even today's online social networks.Such disruptions are so important clinically and economically that they are reported to affect millions of people every year [10,17].As a major consequence of the circadian disruption in the internal desynchronization of day-to-day functionality that regulates important functions such as behavior, metabolism, hormone levels, sleep and body temperature, [10,17].
In this contribution, we focus on a much simpler model 1 with only two nonlinear in-phase coupled oscillators, θ 1 and θ 2 , with the dynamics driven by a Caputo type fractional derivative 2 [5].The main motivation is to test the memory-like hypothesis of the circadian rhythms.New findings and models that incorporated memory in the circadian rhythms might help to propose and design treatments for people wich circadian disruption [10,17].For that fate, we chose to extend the analyses of the dynamical circadian rhythms system of in-phase coupled oscillators proposed by Strogatz [14], since it is a simple model for the sleep-wake θ 1 and body temperature θ 2 circadian rhythms.The mentioned extension deals with the memory imputed into the dynamics by a fractional order derivative (see Section 2 for details).In [3], the authors have contributed to the proposed models numerically, but, as far as we are aware, this contribution is the first to analyze the synchronization phenomena for the fractional phase model mathematically.
The main results of this contribution can be summarized as: • The proposed Caputo fractional order dynamics for in-phase coupled oscillators advance our understanding of memory's influences on circadian rhythms and their synchronization and dessynchronization phenomena; • We show the well-posedness of the proposed dynamical system; • We show that the synchronization phenomena for a pair of in-phase coupled oscillators proposed by Strogatz [14] can be understood if the memory is equally distributed between the oscillators.Hence, extending the conclusions in Strogatz [14]; 1 More circadian rhythm models can be found in [4,8,14] and the references therein. 2 We have chosen the Caputo-type fractional derivatives in these contributions since they implies in biological interpetrable initial condition for the dynamics [5].The analysis of the model in termos of ψ-Hilfer pseudo-fractional operator [13,19,20] will be investigate in future contributions.
Trends Comput.Appl.Math., 24, N. 2 (2023) • The same phenomenon becomes extremely complex in a general picture, if each oscillator has a distinct level of memory.However, given the necessary conditions for the existence of a synchronized solution, lower and upper bounds are derived based solely on the parameters, such that any existing synchronized solution must be confined in between.
• We present numerical realizations supporting the theoretical findings.
Outline: In Section 2, we present the proposed fractional phase model and corresponding coupling topology that will be analyzed during this contribution.In Section 3, we prove the wellposedness of the proposed model.In Section 4 we discuss the synchronization phenomena.In particular, we rapdly review the synchronization analysis of the phase model proposed by Strogatz in Subsection 4.1.We present explicitly the synchronized solutions for the fractional phase model in Subsection 4.2 and the bounds for any existing synchronized solution for the multiorder fractional phase model in Subsection 4.3.The distinct scenarios for synchronization of the proposed models are presented in Section 5. Finally, some conclusions and future directions are summarized in Section 6.

THE FRACTIONAL PHASE MODEL
One essential point in the modeling of the proposed model is the coupling topology.Many coupling strategies have been studied in the literature, for example [9,15] and references therein.Motivated by the circadian rhythm modeling proposed by Strogatz [14], in this contribution we shall consider the coupling topology presented in the Figure 1.By observing the coupling topology in Figure 1, the fractional-order coupled oscillators, θ 1 (the temperature) and θ 2 (the sleep-wake), follow the dynamic for t ∈ [0, f 2 ] and ) In (2.1)-(2.3),ω j are strictly positive constant that represents the intrinsic frequency given by ω j = 1 τ j for the period (in hours) τ j , A, B are the coupling strength and D α j (•) is the Caputo fractional derivative operator [5] of order α j ∈]0, 1], for j = 1, 2, respectively.
The left hand side of equations (2.3)-(2.1)has dimension of (time) −α j , for j = 1, 2, respectively.Therefore, in order of time dimension match (be consistent), the parameters ω j , B, A on the right hand side of (2.3)-(2.4)have powers α j , for j = 1, 2, respectively.See [6] for further details on the dimension consistence of the proposed system.It is also worth noting that, in the limit case of α j → 1, with f 2 = 0 the system (2.3)-(2.4)recover the phase model proposed by Strogatz [14].Hence, we are analyzing a generalization of the model proposed by Strogatz [14], with the premise that our circadian rhythms have memory.Here, memory is recovered thanks to the properties of the fractional Caputo derivative as motivated by [5,Remark 6.4].
Furthermore, we assume that oscillators' movements are counter-clockwise.The signals on the coupling constants are determined to attend to this hypothesis.Therefore, it follows from the schematic system in Figure 1 that the coupling from θ 2 to θ 1 is contrary to the movement.As a result, the signal representing the coupling forces B is negative, whereas the signal representing the coupling A is positive.
The system (2.3)-(2.4) is completed by initial conditions as follows: we assume that the period that we are at sleep is a fraction f 2 of the phase θ 2 (t).Hence, for t 0 = 0 we will have the initial condition or no activities (sleeping).As the Figure 1 suggest, we have In this setting, F 2 > 0 means abrupt changes in the asleep stage.On the other hand, F 2 = 0 means that you have a smooth change between the sleep-asleep stages and can be interpreted as a normally asleep situation.Hence, whatever is necessary in the forthcoming analysis, we assume that F 2 = 0.
Proof.It follows from Lemma 3.1 that item i) and item ii) are true on the interval [0, f 2 ], for with Therefore, is enough to prove the claims for the multi-order fractional differential system (2.3)- It is done similarly to the arguments in Lemma 3.1.
Since The next corollary is an essential part of the synchronization analysis presented in the next section.
Corollary 3.1.Let the assumption on Theorem 3.1 holds true and let f 2 = 0. Then the unique solution of the multi-order fractional system(2.3)-(2.4),with initial conditions (2.7)-(2.5),are given by where In this setting, we will consider synchronization as the well-known emerging phenomenon arising in the analysis of oscillators for which, independently of their differences, there is a spontaneous lock of one to another such that, they will oscillate exactly at the same frequency.In other words, by defining the phase difference we will say that the phase oscillators θ 1 and θ 2 are synchronized when ψ(t) = c, for some constant c.Therefore, this synchronization of two oscillators only make sense for θ 2 ̸ = 0, otherwise it does not influence in θ 1 .Hence, we only analyze the synchronization for with t ≥ f 2 .
Below we will present a complete analysis of analytic synchronization for the case of . Furthermore, we will persuade the reader that spontaneous synchronization for multifractional fractional phase oscillators with ) is an extremely complex phenomenon, as the presence of synchronization is difficult to demonstrate analytically.
On the other hand, we show lower and upper bounds for any existing synchronized solutions.

Synchronization for the Strogatz phase model -
For the sake of completeness, in this subsection, we revisit the Strogatz analysis on synchronization of the phase model (2.3)-(2.4).That means, we are looking for the system of phase oscillators (2.3)-(2.4)with α 1 = α 2 = 1.
Consider the phase difference ψ(t) as in (4.1).It follows from (2.3)-(2.4)(with α 1 = α 2 = 1) that the derivative of ψ(t) w.r.t.t satisfies the equation where Ω = ω 1 − ω 2 is the intrinsic phase difference of the oscillators and D = A + B is the total coupling strength of the system.
As a result, the synchronization of the phase oscillators θ 1 and θ 2 coincide with the stable stationary point of the dynamic (4.2), i.e., ψ ′ (t) = 0. Using this in (4.2), we obtain Notice that the spontaneous synchronization of the phase oscillators implies the following relation between the intrinsic phase difference Ω and the total coupling strength D, given by The relation (4.4) is known as the Winfree constraints for spontaneous synchronization, e.g.[15,22].
Notice that, from (4.5), the frequency during synchronization ω * = θ ′ (1) = θ ′ 2 is given by It differs from the natural frequency ω 1 and ω 2 , by respectively.Hence, we can conclude that, during synchronization, the frequencies of the oscillators are shifted from their intrinsic values in proportion to the coupling strengths, i.e., .9)

Synchronization for the fractional phase model
In this subsection, we will show that there exist an analytic synchronized solution for the fractional system (2.3)-(2.4),with α 1 = α 2 = α ∈]0, 1[.The idea is to follow the same steps on Subsection (4.1).
According to (2.3)-(2.4)we have that where Ω α = ω α 1 − ω α 2 and D α = A α + B α are the α-power intrinsic frequency differences and the total coupling strength of the system, respectively.Because Caputo's derivatives of constants are zero (see the definition of Caputo derivative in [5], synchronization of the oscillators θ 1 and θ 2 corresponds to the stable stationary points of (4.10), i.e.D α ψ 1 (t) = 0.In other words, the stationary phase difference satisfies on the synchronization stage.Using the fractional integral of order α (see [5]) on both side of (4.12) and the initial conditions, we have Integrating the above equations implies that, during synchronization It follows from (4.12) that the frequency during synchronization ω * α = D α θ (1) = D α θ 2 is given by that differs from the natural frequency ω α 1 and ω α 2 , by respectively.Similarly to the analysis on Subsection 4.1, the oscillator frequencies are shifted from their intrinsic values in proportion to the fractional power coupling strengths, i.e., during synchronization.
From the analysis carried out before, we can argue about the influence of the fractional order dynamics on the synchronization of the phase oscillator model (2.3)-(2.4).Indeed, from (4.9), the oscillator frequencies are shifted from their intrinsic values in a linear proportion to the coupling strengths for the case of α = 1 (the Strogatz model) and as a fractional power low for the case of α 1 = α 2 = α.It follows that the frequency during synchronization is a monotone decreasing (increasing) function of α ∈]0, 1] towards B A , if A > B (respectively, B > A).

Richness in Synchronization for the multi-order fractional phase model
In this subsection, we take a look on the richness of the dynamics of the multi-order fractional phase model (2.3)-(2.4)with α 1 ̸ = α 2 .We will provide pieces of evidence that, if the dynamics of our body is governed by such a model, then, the complexity of our circadian rhythm synchronization is hidden behind a simple explanation as pointed out in Subsections (4.1)-(4.2).However, some glimpses of the relationship between the coupling strengths and internal frequencies show up in the calculations.
For ease of presentation, we assume that f 2 = 0, otherwise, we shall analyze the synchronization for t > f 2 .Considering the difference of (3.1) and (3.2) given in Corollary 3.1) and using a standard calculation, we have For the next steps we shall assume the following: Assumption 4.1.g 1 (t) − g 2 (t) and g 2 (t) does not change signal, during synchronization.
Before we continue with the analysis, it is worth noting some facts in relation to the Assumption (4.1).
Remark 4.2.Let g 2 (t) > g 1 (t) in Assumption (4.1).Then, it follows from the definition of g 1 (t) and g 2 (t) that during the synchronization.Consequently, we have that is a generalization of Winfree constraints for synchronization, e.g., [22].
On the other hand, the requirement of g 2 (t) does not change the signal during synchronization in Assumption 4.1, can be interpreted as the weakly coupled ideas necessary in the analysis of Kuramoto's and Winfree's approaches [9,15,22], since A α 2 shall be small with relation to ω α 2 2 .If the Assumption 4.1 holds true, we can apply the Mean Value Theorem [6, Theorem 2.1], from with, we conclude that there exists τ ∈]0,t] such that, equation (4.16) is rewritten as where J α 2 0 denotes the fractional integral operator of order α 2 ∈]0, 1], [5].Defining the phase difference ψ 2 (t) := (θ 1 − θ 2 )(t), we have Next, we will apply the chain rule for fractional derivative [5] to the equation (4.17).Since the chain rule for fractional derivative involves an infinity series, we assume that such a series is convergent.Furthermore, to have some glimpses into the relationship between the coupling strengths and the natural frequencies for the synchronization of the oscillators, we present the lower order terms of the series explicitly.It means that we have Trends Comput.Appl.Math., 24, N. 2 (2023) Using that D α 2 −1 = J 1−α 2 for α 2 ∈]0, 1], see [5] and developing the terms for k = 1 above, we arrive in As before, the synchronized solution satisfies D α 2 ψ 2 (t) = 0.However, differently from the cases in Subsections 4.1-4.2,there is a richness of factors around the coupling strengths and natural frequencies, since the expression is dependent on t.Therefore, it is not clear whether to have an explicit synchronized solution.
Because of the aforementioned analysis, we conjecture that the synchronization of two oscillators with different memories depends on other external forces known in the chronobiology and medicine literature as zeitgebers, e.g.[1,2,8,14].It will be investigated in the future.
In the remaining part of this section, we will use the comparison principle [21, Theorem 2] together with the necessary conditions of Assumption 4.1 to derive bounds that encapsulate the synchronized solutions for the oscillators.It will be described in the following proposition.

NUMERICAL EXPERIMENTS
In this section, we present some numerical interpretation of the synchronization analysis carried out during the writing of this manuscript.The parameters of the model are selected to satisfy the Winfree restrictions where available (i.e., for the phase model and the fractional phase model).
In any of the cases, the parameters are calibrated from real data.Calibration (identification) will be investigated in future contributions.For the multi-fractional phase model, where explicit synchronized solutions are not available, we solve the multi-fractional differential equations (2.3)-(2.4)numerically, using the numerical scheme proposed in [7] for solving multi-term fractional differential equations.
It is worth mentioning that, since we are looking for synchronized solutions, in all the numerical experiments we assume that f 2 = F 2 = 0. Furthermore, the intrinsic frequencies ω i of each oscillator θ i for i = 1, 2 is related to the period τ i , in hours, for with each oscillator complete their on cycle.In other words, ω i = 1 τ i .In the next table, we present the parameter value choices for each numerical example 3 .Figure 4 correspond to synchronized solutions of the fractional phase model, for the fractional order of Caputo derivatives given by α 1 = α 2 = 0, 9 and the same intrinsic frequencies and coupling forces used in the solution presented in Figure 3 (see the second and third line of Table (1).This example collaborates with the thesis that memory induces synchronization.It is  On the other hand, Figure 5 shows a desynchronized solution of the fractional phase model (see the fourth line of Table 1 for the parameter choices).It is also important to mention that the simulated situation does not satisfy the generalized Winfree restriction for synchronization in Remark 4.1.Given the numerical examples presented above, it seems that the generalized Winfree restriction (see Remark 4.1) is a bifurcation point for the synchronization.However, the conclusion of such an observation will not be analyzed in this contribution.
In Figures 6 and 7, we present a synchronized and a desynchronized solution for the multi-order fractional phase model studied in Subsection 4.3, respectively.Since for this approach, we were unable to derive analytic synchronized solutions, the parameter choices (given by the fifty and sixth lines of Table 1, respectively) were considered satisfying a generalized Winfree's kind of restrictions, given by A , respectively.In particular, Figure 6 numerically verifies the conclusions on Proposition 4.1.Figure 6: Solutions θ i , θ iL and θ iU (as pointed in Proposition 4.1) of the multi-order fractional phase model with the parameter choice in the fifth line on Table 1.The phase difference is constant, indicating that the oscillators are synchronized.The generalized Winfree restriction is unavailable but it satisfies | what we conjecture that indicates the synchronization.
Figure 7: Solutions θ 1 and θ 2 of the multi-order fractional phase model with the parameter choice in the sixth line on Table 1.The phase difference is not constant, indicating that the oscillators are desynchronized.The generalized Winfree restriction is unavailable, but it satisfies | what we conjecture that indicates the desynchronization.

CONCLUSIONS AND FUTURE WORKS
In this contribution, we have extended the coupling oscillators phase model proposed by Strogatz [14] employing a Caputo fractional derivative dynamics (fractional coupling oscillators phase model).Such an approach has the potential to reveal some features related to the internal synchronization/desynchronization of the circadian rhythm with memory.
We have shown the well-posedness of the proposed model as well as provided explicitly the synchronized solutions for the fractional phase model.We also derived the Winfree restrictions, which allow attesting synchronization using only the system's parameters.We obtain lower and upper bounds for the multi-order fractional phase model, for any existing synchronized solution must be confined in between.We also provide numerical simulations supporting the theoretical findings.
Many questions regarding the present approach remain open and will be investigated in future contributions.
• The calibration (identification) of the model parameters from real data.
• Analyse the Winfree restriction (see Remark 4.1) as bifurcation point for the synchronization.
• Improve the understanding of synchronization of multi-order fractional phase models.
• Investigate the synchronization of the multi-order fractional phase model and its relation with external forces, known in the chronobiology and medicine literature as zeitgebers, e.g.[1,2,8,14].

Figure 1 :
Figure 1: Schematic coupling topology of the oscillators system.

Figures 2 and 3
Figures 2 and 3 represents synchronized and desynchronized solutions of the Strogatz phase model discussed in Subsection 4.1, respectively.The parameters for the model correspond to the Winfree restriction (4.4), respectively.

Figure 2 :
Figure 2: Solutions θ 1 and θ 2 of the phase model with parameter choice corresponding the first line on Table (1).The phase difference is constant, indicating the synchronization.The Winfree restriction (4.4) is satisfied.

Figure 3 :
Figure 3: Solutions θ 1 and θ 2 of the phase model with the parameter choice in the second line onTable(1).The phase difference is not constant, indicating that the oscillators are desynchronized.The Winfree restriction (4.4) is not satisfied, which also indicates the desynchronization.

Figure 4 :
Figure 4: Solutions θ 1 and θ 2 of the fractional phase model with the parameter choice in the third line on Table (1).The phase difference is constant indicating that the oscillators are synchronized.The generalized Winfree restriction (see Remark 4.1) is satisfied, which indicates the synchronization.

Figure 5 :
Figure 5: Solutions θ 1 and θ 2 of the fractional phase model with the parameter choice in the fourth line on Table 1.The phase difference is not constant indicating that the oscillators are desynchronized.The generalized Winfree restriction (see Remark (4.1) is not satisfied, which indicates the desynchronization.

Figure 8 :
Figure 8: Asymptotic behaviour of the solutions θ i , θ iL and θ iU of the multi-order fractional phase model as α 1 → α 2 .
[5,os(x)| ≤ |x| + 1 it follows that the right hand side of the system (2.3)-(2.4) is continuous w.r.t.t and Lipschitz continuous w.r.t. the second argument.It follows from[5, Theorem 8.11]the existence of a unique solution as claimed in item i).Following the steps in [5, Theorem 8.9] one can see that the system (2.3)-(2.4)canbetransformed in a system with the same fractional derivative order in each equation, for with, if the right-hand side is continuous and Lipschitz continuous the transformed system also is continuous and Lipshitz continuous.Hence, one can apply now[5, Theorem 8.11]to conclude the existence of a unique continue solution on [ f 2 , K * ] that depends continuously on the initial data, the parameters, and the derivative orders as claimed in i) ii).
The extension of the solution of (2.3)-(2.4)to the interval [0, f 2 ] was already proved in Lemma 3.1.Again, since cos(x) is uniformly bounded, it follows that the right hand side of the system (2.3)-(2.4),with initial conditions on f 2 as above, easily satisfies [18, Theorem 3.], for with the continuously extension to the interval [ f 2 , ∞[ follows.It concludes the claim of item iii).□ It is worth noting that if F 2 > 0, then we we cannot prove the continuity of θ 2 in [0, ∞[, becouse θ 2 (t) = 0 for t ∈ [0, f 2 [.However, we could consider a piecewise continue solution (θ 1 (t), 0) T on [0, f 2 [ given by Lemma 3.1 and (θ 1 (t), θ 2 (t)) T on [ f 2 , ∞[, which follows the line of proof of Theorem 3.1.It will appear explicitly in the next sections.

Table 1 :
Table with the parameter values of the corresponding numerical examples.