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 Strogatz¹⁴ 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.

KEYWORDS:
circadian rhythms; synchronization; memory

1 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⁹, 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²² and Kuramoto⁹. 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¹ with only two nonlinear in-phase coupled oscillators, θ ₁ and θ ₂, 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¹⁴, since it is a simple model for the sleep-wake θ ₁ and body temperature θ ₂ circadian rhythms. The mentioned extension deals with the memory imputed into the dynamics by a fractional order derivative (see Section 2 for details). In³, the authors have contributed to the proposed models numerically, but, as far as we are aware, this contribution is the first to analyze th e 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 ¹⁴ can be understood if the memory is equally distributed between the oscillators. Hence, extending the conclusions in Strogatz¹⁴;

• 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 well-posedness 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 multi-order 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.

2 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¹⁴, in this contribution we shall consider the coupling topology presented in the Figure 1.

Figure 1
Schematic coupling topology of the oscillators system.

By observing the coupling topology in Figure 1, the fractional-order coupled oscillators, θ ₁ (the temperature) and θ ₂ (the sleep-wake), follow the dynamic

$D^{{\alpha }_1}{\theta }_1\left(t\right)={\omega }_1^{{\alpha }_1}-B^{{\alpha }_1}\cos \left(2\pi \left({\theta }_1\left(t\right)\right)\right)$ (2.1)

${\theta }_2\left(t\right)=0,$ (2.2)

for t∈[0, f ₂] and

$D^{{\alpha }_1}{\theta }_1\left(t\right)={\omega }_1^{{\alpha }_1}-B^{{\alpha }_1}\cos \left(2\pi \left({\theta }_2\left(t\right)-{\theta }_1\left(t\right)\right)\right)$ (2.3)

$D^{{\alpha }_2}{\theta }_2\left(t\right)={\omega }_2^{{\alpha }_2}+A^{{\alpha }_2}\cos \left(2\pi \left({\theta }_1\left(t\right)-{\theta }_2\left(t\right)\right)\right)$ (2.4)

for t≥f ₂.

In (2.1)-(2.3), ω _j are strictly positive constant that represents the intrinsic frequency given by ${\omega }_j=\frac{1}{{\tau }_j}$ for the period (in hours) τ _j , A, B are the coupling strength and D ^αj (·) is the Caputo fractional derivative operator⁵ 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⁶ 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 ₂=0 the system (2.3)-(2.4) recover the phase model proposed by Strogatz¹⁴. Hence, we are analyzing a generalization of the model proposed by Strogatz¹⁴, 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 [⁵, 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 θ ₂ to θ ₁ 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 ₂ of the phase θ ₂(t). Hence, for t ₀=0 we will have the initial condition

${\theta }_2\left(t_0=0\right)=0,$ (2.5)

or no activities (sleeping). As the Figure 1 suggest, we have

${\theta }_2\left(t\right)=0fort\in \left[0,f_2\right[and{\theta }_2\left(t=f_2\right)=F_2,$ (2.6)

In this setting, F ₂>0 means abrupt changes in the asleep stage. On the other hand, F ₂=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 ₂=0.

Furthermore, we assume that the temperature before the activity is lower. Therefore, after normalization, we assume that

${\theta }_1\left(t_0=0\right)=0.$ (2.7)

The system (2.3)-(2.4) is a multi-order fractional differential system with initial conditions (2.7)-(2.5) e.g.,⁵. In particular, for α ₁=α ₂=1 and f ₂=0, we have the phase coupled oscillators proposed by Strogatz¹⁴. We will refer the model (2.1)-(2.3) as the Strogatz phase model when α ₁=α ₂=1, fractional phase model when α ₁=α ₂ ou multi-order fractional phase model when α ₁≠α ₂.

3 WELL-POSEDNESS

In this section, we prove the well-posedness of the fractional system of nonlinear coupled oscillators (2.3)-(2.4) with initial conditions given by (2.7)-(2.5) with α _j ∈]0, 1], for j=1, 2.

Lemma 3.1.The fractional order initial value problem(2.1)-(2.7)has a unique continue solution on [0, f ₂]. It solution depends continuously on the parameters, the fractional order α ₁∈]0, 1] and on the initial condition-(2.7).

Proof. Since |cos(x)|≤1≤1+|x|, it follows that the right hand side of the system (2.1) is continuous w.r.t.t and Lipschitz continuous w.r.t. the second argument. Therefore, [⁵, Theorem 6.1 - Theorem 6.5] implies the existence of a unique solution on the interval $\left[0,{\hat{K}}_{*}\right]$, for some ${\hat{K}}_{*}>0$, that depends only on the parameters, the Lipschitz constant and on α ₁∈]0, 1].

Moreover, since cos(x) is uniformly bounded, it follows that cos(x)≤1+|x| ^µ , for some µ∈[0,1]. Hence, the right hand side of (2.1) easily satisfies [¹⁸, Theorem 3.2], for with the continuous extension of a solution of (2.1)-(2.7) to the interval [0, f ₂] follows.

Now, we can formulate the following result

Theorem 3.1.Let α_j ∈]0, 1], for j=1, 2, respectively. Assume F ₂=0. Then:

1. There exist a unique continuous solution (θ ₁(t), θ ₂(t))^T for multi-order fractional differential system(2.3)-(2.4)with initial conditions(2.7)-(2.5)on the interval [0, K _*], for some K _* >0.

2. The solution (θ ₁(t), θ ₂(t))^T of the multi-order fractional differential system(2.3)-(2.4)depends continuously on the initial conditions(2.7)-(2.5), on the parameters of the system and on the order of the derivatives α _j ∈]0, 1], for j=1, 2.

3. The solution (θ ₁(t), θ ₂(t))^T of the multi-order fractional differential system(2.3)-(2.4)can be continuously extended to the interval [0, ∞].

Proof. It follows from Lemma 3.1 that item i) and item ii) are true on the interval [0, f ₂], for with θ ₂(t) is identically zero on [0, f ₂]. Furthermore, θ ₁(f ₂) is well defined.

Therefore, is enough to prove the claims for the multi-order fractional differential system (2.3)-(2.4) with initial conditions θ ₁(t ₀)=θ ₁(f ₂) and θ ₂(t ₀)=F ₂:=0.

It is done similarly to the arguments in Lemma 3.1.

Since |cos(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 [⁵, Theorem 8.11] the existence of a unique solution as claimed in item i). Following the steps in [⁵, Theorem 8.9] one can see that the system (2.3)-(2.4) can be transformed 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 [⁵, Theorem 8.11] to conclude the existence of a unique continue solution on [f ₂ , 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 ₂] 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 ₂ as above, easily satisfies [¹⁸, Theorem 3.], for with the continuously extension to the interval [f ₂ , ∞[ follows. It concludes the claim of item iii).

It is worth noting that if F ₂>0, then we we cannot prove the continuity of θ ₂ in [0, ∞[, because θ ₂(t)=0 for t∈[0, f ₂[. However, we could consider a piecewise continue solution (θ ₁(t), 0)^T on [0, f ₂[ given by Lemma 3.1 and (θ ₁(t), θ ₂(t))^T on [f ₂ , ∞[, which follows the line of proof of Theorem 3.1. It will appear explicitly in the next sections.

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₂=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

${\theta }_1\left(t\right)=\frac{1}{\Gamma \left({\alpha }_1\right)}{\int }_0^t{\left(t-s\right)}^{{\alpha }_1-1}{g}_1\left(s\right)ds$ (3.1)

${\theta }_2\left(t\right)=\frac{1}{\Gamma \left({\alpha }_2\right)}{\int }_0^t{\left(t-s\right)}^{{\alpha }_2-1}{g}_2\left(s\right)ds$ (3.2)

where

$g_1\left(t\right)={\omega }_1^{{\alpha }_1}-B^{{\alpha }_1}\cos \left(2\pi \left({\theta }_2\left(t\right)-{\theta }_1\left(t\right)\right)\right)$ (3.3)

$g_2\left(t\right)={\omega }_2^{{\alpha }_2}+A^{{\alpha }_2}\cos \left(2\pi \left({\theta }_1\left(t\right)-{\theta }_2\left(t\right)\right)\right)$ (3.4)

Proof. The existence and uniqueness follows from Theorem 3.1. Then, the claim of this corollary follows from a component-wise application of [⁵, Lemma 6.2] in (2.3)-(2.4) with the null initial conditions (2.7)-(2.5).

4 ON THE SYNCHRONIZATION ANALYSIS

In this section, we investigate the synchronization of the fractional phase model (2.3)-(2.4). More precisely, we will investigate the macroscopic synchronization phenomena associated with the model (2.3)-(2.4) for different choices of the fractional derivatives α ₁ , α ₂∈]0, 1].

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

$\psi \left(t\right)={\theta }_1\left(t\right)-{\theta }_2\left(t\right),$ (4.1)

we will say that the phase oscillators θ ₁ and θ ₂ are synchronized when ψ(t)=c, for some constant c. Therefore, this synchronization of two oscillators only make sense for θ ₂≠0, otherwise it does not influence in θ ₁. Hence, we only analyze the synchronization for with t≥f ₂.

Below we will present a complete analysis of analytic synchronization for the case of α ₁=α ₂∈]0, 1]. Furthermore, we will persuade the reader that spontaneous synchronization for multifractional fractional phase oscillators with α ₁≠α ₂ (although α ₁ , α ₂∈]0, 1]) 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.

4.1 Synchronization for the Strogatz phase model - α₁=α₂=1

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.

Consider the phase difference ψ(t) as in (4.1). It follows from (2.3)-(2.4) (with α ₁=α ₂=1) that the derivative of ψ(t)w.r.t.t satisfies the equation

$\psi \text{'}\left(t\right)=\Omega -D\cos \left(2\pi \psi \left(t\right)\right),$ (4.2)

where Ω=ω ₁ −ω ₂ 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 θ ₁ and θ ₂ coincide with the stable stationary point of the dynamic (4.2), i.e., ψ′(t)=0. Using this in (4.2), we obtain

$\psi \left(t\right)=\frac{1}{2\pi }\arccos \left(\frac{\Omega }{D}\right).$ (4.3)

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

$K:=\left|\frac{D}{\Omega }\right|>1,$ (4.4)

The relation (4.4) is known as the Winfree constraints for spontaneous synchronization, e.g.^{15), (22}.

Using (4.3) in the system (2.3)-(2.4) (with α ₁=α ₂=1), we have that, during synchronization, the dynamics satisfies

${\theta }_1^{\text{'}}\left(t\right)=\frac{A{\omega }_1+B{\omega }_2}{A+B}={\theta }_2^{\text{'}}\left(t\right).$ (4.5)

Integrating both sides of (4.5) w.r.t.t and using the initial conditions (2.7)-(2.5), as well as the fact that θ ₂(t)=0 on the interval [0, f ₂], we find that the synchronized analytic solution of (2.3)-(2.4) is given by

${\theta }_1\left(t\right)=\left(\frac{A{\omega }_1+B{\omega }_2}{A+B}\right)t$ (4.6)

${\theta }_2\left(t\right)=\left(\frac{A{\omega }_1+B{\omega }_2}{A+B}\right)\left(t-f_2\right)+F_2,$ (4.7)

It follows from (4.1) and (4.6)-(4.7) that, in the synchronization stage, the phase differences, $\psi \left(t\right)=f_2\left(\frac{A{\omega }_1+B{\omega }_2}{A+B}\right)-F_2$ is constant. This influence of f ₂ and θ ₂(f ₂) was not capture in the original paper by Strogatz¹⁴.

Notice that, from (4.5), the frequency during synchronization $\omega *=\theta \text{'}\left(1\right)={\theta }_2^{\text{'}}$ is given by

${\omega }^{*}=\frac{A{\omega }_1+B{\omega }_2}{A+B}.$

It differs from the natural frequency ω ₁ and ω ₂, by

$\begin{array}{c}\delta {\omega }_1={\omega }^{*}-{\omega }_1=-\frac{B\Omega }{D},\\ \delta {\omega }_2={\omega }^{*}-{\omega }_2=\frac{A\Omega }{D},\end{array}$ (4.8)

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.,

$\left|\frac{\delta {\omega }_1}{\delta {\omega }_2}\right|=\frac{B}{A}.$ (4.9)

4.2 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 α ₁=α ₂=α∈]0, 1[. The idea is to follow the same steps on Subsection (4.1).

Lets take the Caputo fractional derivative of order α∈]0, 1[, of the phase difference ψ(t). According to (2.3)-(2.4) we have that

$D^{\alpha }{\psi }_1\left(t\right)={\Omega }_{\alpha }-D_{\alpha }\cos \left(2\pi {\psi }_1\left(t\right)\right)$ (4.10)

where ${\Omega }_{\alpha }={\omega }_1^{\alpha }-{\omega }_2^{\alpha }$ 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⁵, synchronization of the oscillators θ ₁ and θ ₂ corresponds to the stable stationary points of (4.10), i.e. D ^α ψ ₁(t)=0. In other words, the stationary phase difference satisfies

${\psi }_1\left(t\right)=\frac{1}{2\pi }\arccos \left(\frac{{\Omega }_{\alpha }}{D_{\alpha }}\right).$ (4.11)

Remark 4.1.Notice that, since the domain of arccos∈[−1, 1], it follows from(4.11), that a necessary condition for synchronization is that |Ω_α |<D _α . This condition is a generalization of Winfree’s restriction, e.g.^{15), (22}, to the fractional phase model.

Therefore, substituting (4.11) into the system (2.3)-(2.4), we have that

$D^{\alpha }{\theta }_1\left(t\right)=\frac{A^{\alpha }{\omega }_1^{\alpha }+B^{\alpha }{\omega }_2^{\alpha }}{A^{\alpha }+B^{\alpha }}=D^{\alpha }{\theta }_2\left(t\right),$ (4.12)

on the synchronization stage. Using the fractional integral of order α (see⁵) on both side of (4.12) and the initial conditions, we have

${\theta }_1\left(t\right)=\left(\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{\left(t-x\right)}^{\alpha -1}\left(\frac{A^{\alpha }{\omega }_1^{\alpha }+B^{\alpha }{\omega }_2^{\alpha }}{A^{\alpha }+B^{\alpha }}\right)dx\right)$

and

${\theta }_2\left(t\right)=\left({\theta }_2\left(f_2\right)+\frac{1}{\Gamma \left(\alpha \right)}{\int }_{f_2}^t{\left(t-x\right)}^{\alpha -1}\left(\frac{A^{\alpha }{\omega }_1^{\alpha }+B^{\alpha }{\omega }_2^{\alpha }}{A^{\alpha }+B^{\alpha }}\right)dx\right),$

since θ ₂(t)=0 for t∈[0, f ₂[.

Integrating the above equations implies that, during synchronization

$\begin{array}{c}{\theta }_1\left(t\right)=\frac{\left(A^{\alpha }{\omega }_1^{\alpha }+B^{\alpha }{\omega }_2^{\alpha }\right)t^{\alpha }}{\left(A^{\alpha }+B^{\alpha }\right)\alpha \Gamma \left(\alpha \right)}\\ {\theta }_2\left(t\right)=\frac{\left(A^{\alpha }{\omega }_1^{\alpha }+B^{\alpha }{\omega }_2^{\alpha }\right){\left(t-f_2\right)}^{\alpha }}{\left(A^{\alpha }+B^{\alpha }\right)\alpha \Gamma \left(\alpha \right)}+F_2.\end{array}$ (4.13)

It follows from (4.12) that the frequency during synchronization ${\omega }_{\alpha }^{*}=D^{\alpha }\theta \left(1\right)=D^{\alpha }{\theta }_2$ is given by

${\omega }_{\alpha }^{*}=\frac{A^{\alpha }{\omega }_1^{\alpha }+B^{\alpha }{\omega }_2^{\alpha }}{A^{\alpha }+B^{\alpha }},$

that differs from the natural frequency ${\omega }_1^{\alpha }$ and ${\omega }_2^{\alpha }$, by

${\delta }_{\alpha }{\omega }_1={\omega }_{\alpha }^{*}-{\omega }_1^{\alpha }=-\frac{B^{\alpha }{\Omega }_{\alpha }}{D_{\alpha }},{\delta }_{\alpha }{\omega }_2={\omega }_{\alpha }^{*}-{\omega }_2^{\alpha }=\frac{A^{\alpha }{\Omega }_{\alpha }}{D_{\alpha }},$ (4.14)

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.,

$\left|\frac{{\delta }_{\alpha }{\omega }_1}{{\delta }_{\alpha }{\omega }_2}\right|={\left(\frac{B}{A}\right)}^{\alpha },$ (4.15)

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 α ₁=α ₂=α. It follows that the frequency during synchronization is a monotone decreasing (increasing) function of α∈]0, 1] towards $\frac{B}{A}$, if A>B (respectively, B>A).

4.3 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 α ₁≠α ₂. 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 ₂=0, otherwise, we shall analyze the synchronization for t>f ₂. Considering the difference of (3.1) and (3.2) given in Corollary 3.1) and using a standard calculation, we have

$\left({\theta }_1-{\theta }_2\right)\left(t\right)={\int }_0^t\frac{\left(t-s\right){}^{{\alpha }_2-1}}{\Gamma \left({\alpha }_2\right)}\left[g_1\left(s\right)-g_2\left(s\right)\right]\left[\frac{\Gamma \left({\alpha }_2\right)}{\Gamma \left({\alpha }_1\right)}\left(t-s\right){}^{{\alpha }_1-{\alpha }_2}\right]ds+{\int }_0^t\frac{{\left(t-s\right)}^{{\alpha }_2-1}}{\Gamma \left({\alpha }_2\right)}{g}_2\left(s\right)\left[\frac{\Gamma \left({\alpha }_2\right)}{\Gamma \left({\alpha }_1\right)}{\left(t-s\right)}^{{\alpha }_1-{\alpha }_2}-1\right]ds.$ (4.16)

For the next steps we shall assume the following:

Assumption 4.1.g₁(t)−g₂(t) and g₂(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 ₂ (t)>g ₁ (t) in Assumption (4.1). Then, it follows from the definition of g ₁ (t) and g ₂ (t) that

$A^{{\alpha }_2}+B^{{\alpha }_1}\cos \left(2\pi \left({\theta }_2\left(t\right)-{\theta }_1\left(t\right)\right)\right)\ge {\omega }_1^{{\alpha }_1}-{\omega }_2^{{\alpha }_2}$

during the synchronization. Consequently, we have

$A^{{\alpha }_2}+B^{{\alpha }_1}\ge \left|{\omega }_1^{{\alpha }_1}-{\omega }_2^{{\alpha }_2}\right|$

that is a generalization of Winfree constraints for synchronization, e.g., ²² .

On the other hand, the requirement of g ₂ (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 ${\omega }_2^{{\alpha }_2}$ .

If the Assumption 4.1 holds true, we can apply the Mean Value Theorem [⁶, Theorem 2.1], from with, we conclude that there exists τ∈]0,t] such that, equation (4.16) is rewritten as

$\left({\theta }_1-{\theta }_2\right)\left(t\right)=\frac{\Gamma \left({\alpha }_2\right)}{\Gamma \left({\alpha }_1\right)}{\left(t-\tau \right)}^{{\alpha }_1-{\alpha }_2}{J}_0^{{\alpha }_2}\left(\left(g_1-g_2\right)\left(t\right)\right)+\left(\frac{\Gamma \left({\alpha }_2\right)}{\Gamma \left({\alpha }_1\right)}{\left(t-\tau \right)}^{{\alpha }_1}-1\right)J_0^{{\alpha }_2}\left(g_2\left(t\right)\right),$

where $J_0^{{\alpha }_2}$ denotes the fractional integral operator of order α₂∈]0, 1],⁵. Defining the phase difference ψ ₂(t):=(θ ₁ −θ ₂)(t), we have

${\psi }_2\left(t\right)=\frac{\Gamma \left({\alpha }_2\right)}{\Gamma \left({\alpha }_1\right)}{\left(t-\tau \right)}^{{\alpha }_1-{\alpha }_2}{J}_0^{{\alpha }_2}\left(\left(\left({\omega }_1^{{\alpha }_1}-{\omega }_2^{{\alpha }_2}\right)-\left(A^{{\alpha }_1}+B^{{\alpha }_2}\right)\cos \left(2\pi {\psi }_2\left(t\right)\right)\right)\right)+\frac{\Gamma \left({\alpha }_2\right)}{\Gamma \left({\alpha }_1\right)}\left({\left(t-\tau \right)}^{{\alpha }_1}-1\right)J_0^{{\alpha }_2}\left(\left({\omega }_2^{{\alpha }_2}+A^{{\alpha }_2}\cos \left(2\pi {\psi }_t\left(s\right)\right)\right)\right)$ (4.17)

Next, we will apply the chain rule for fractional derivative⁵ 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

$D^{{\alpha }_2}{\psi }_2\left(t\right)=k=0\frac{\Gamma \left({\alpha }_2\right)}{\Gamma \left({\alpha }_1\right)}\left({\left(t-\tau \right)}^{{\alpha }_1-{\alpha }_2}\left(g_1\left(t\right)-g_2\left(t\right)\right)+\left({\left(t-\tau \right)}^{{\alpha }_1}-1\right)g_2\left(t\right)\right)k=1+\frac{\Gamma \left({\alpha }_2+1\right)}{\Gamma \left({\alpha }_1\right)}\left[D^1{\left(t-\tau \right)}^{{\alpha }_1-{\alpha }_2}\right]\left(D^{{\alpha }_2-1}\left[J^{{\alpha }_2}\left(g_1\left(t\right)-g_2\left(t\right)\right)\right]\right)+\left[D^1\left({\left(t-\tau \right)}^{{\alpha }_1}-1\right)\right]\left[D^{{\alpha }_2-1}\left(J^{{\alpha }_2}{g}_2\left(t\right)\right)\right]+higherorderterms,k\ge 2$

Using that $D^{{\alpha }_2-1}=J^{1-{\alpha }_2}$ for α ₂∈]0, 1], see⁵ and developing the terms for k=1 above, we arrive in

$D^{{\alpha }_2}{\psi }_2\left(t\right)=\frac{\Gamma \left({\alpha }_2\right)}{\Gamma \left({\alpha }_1\right)}\left({\left(t-\tau \right)}^{{\alpha }_1-{\alpha }_2}\left(g_1\left(t\right)-g_2\left(t\right)\right)+\left({\left(t-\tau \right)}^{{\alpha }_1}-1\right)g_2\left(t\right)\right)+\frac{\Gamma \left({\alpha }_2+1\right)}{\Gamma \left({\alpha }_1\right)}{\left(t-\tau \right)}^{{\alpha }_1-\left(1+{\alpha }_2\right)}\left[\left({\alpha }_2-{\alpha }_1\right){\int }_0^t\left(g_1\left(s\right)-g_2\left(s\right)\right)ds+{\alpha }_1{\left(t-\tau \right)}^{-\left(1+{\alpha }_2\right)}{\int }_0^t{g}_2\left(s\right)ds\right]+higherorderterms,k\ge 2$

As before, the synchronized solution satisfies $D^{{\alpha }_2}{\psi }_2\left(t\right)=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 [²¹, 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.

Proposition 4.1.Let Assumption 4.1 holds true with g₁(t)≥g ₂(t) and g ₂(t)≥0. Then, the unique synchronized solution θ ₁(t), θ ₂(t) of(2.3)-(2.4)are enclosed by θ _jL (t)≤θ _j (t)≤θ _jU (t), for j=1, 2, given, respectively, by

${\theta }_{lL}\left(t\right)=\frac{A^{{\alpha }_2}{\omega }_1^{{\alpha }_1}+B^{{\alpha }_1}{\omega }_2^{{\alpha }_2}}{{\alpha }_1\Gamma \left({\alpha }_1\right)\left(A^{{\alpha }_2}+B^{{\alpha }_1}\right)}{t}^{{\alpha }_1},{\theta }_{2L}\left(t\right)=0$ (4.18)

and

${\theta }_{lU}\left(t\right)=\frac{A^{{\alpha }_2}{\omega }_1^{{\alpha }_1}+B^{{\alpha }_1}{\omega }_2^{{\alpha }_2}}{{\alpha }_1\Gamma \left({\alpha }_1\right)A^{{\alpha }_2}}{t}^{{\alpha }_1},{\theta }_{2U}\left(t\right)=\frac{A^{{\alpha }_2}{\omega }_1^{{\alpha }_1}+B^{{\alpha }_1}{\omega }_2^{{\alpha }_2}}{{\alpha }_2\Gamma \left({\alpha }_2\right)\left(A^{{\alpha }_2}+B^{{\alpha }_1}\right)}{\left(t-f_2\right)}^{{\alpha }_2}+F_2.$ (4.19)

Proof. The uniqueness follows directly from Theorem 3.1. Then, from Assumption 4.1, it follows that

$-\frac{{\omega }_2^{{\alpha }_2}}{A^{{\alpha }_2}}\le \cos \left(2\pi {\psi }_2\right)\le \frac{{\omega }_1^{{\alpha }_1}-{\omega }_2^{{\alpha }_2}}{A^{{\alpha }_2}+B^{{\alpha }_1}}.$ (4.20)

Hence, the claim follows by substitute (4.20) back in (2.3) and (2.4) and applying the comparison principle [²¹, Theorem 2].

An interesting conjecture about the synchronization for α ₁ ≠α ₂ is possible by compare the result in Proposition 4.1 with the synchronized solution for α ₁=α ₂ in (4.13): there exists an external mechanism that slow down the oscillator θ ₁ into the direction of θ _1L and speed-up the oscillator θ ₂ into the direction of θ _2U as |α ₁ −α ₂ |→0, asymptotically.

5 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⁷ 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 ₂=F ₂=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, ${\omega }_i=\frac{1}{{\tau }_i}$.

In the next table, we present the parameter value choices for each numerical example³.

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

Example	Parameter Value						Winfree restrictions
	A	B	τ1	τ2	α1	α2	θ1 and θ2
1	0,016	0,03	22	24	1	1	synchronized by (4.4)
2	0,008	0,009	16	23	1	1	desynchronized by (4.4)
3	0,008	0,009	16	23	0,9	0,9	synchronized by Remark 4.1
4	0,0005	0,009	18	26	0,6	0,6	desynchronized by Remark 4.1
5	0,003	0,008	20	26	0,5	0,4	unavailable
6	0,006	0,004	20	24	0,5	0,8	unavailable

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
Solutions θ ₁ and θ ₂ 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
Solutions θ ₁ and θ ₂ of the phase model with the parameter choice in the second line on Table (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 correspond to synchronized solutions of the fractional phase model, for the fractional order of Caputo derivatives given by α ₁=α ₂=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 worth mentioning that the parameter choices satisfy the generalized Winfree restrictions derived from Remark 4.1.

Figure 4
Solutions θ ₁ and θ ₂ 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.

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.

Figure 5
Solutions θ ₁ and θ ₂ 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.

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^{{\alpha }_1}+B^{{\alpha }_2}>\left|{\omega }_1^{{\alpha }_1}-{\omega }_2^{{\alpha }_2}\right|$ and $A^{{\alpha }_1}+B^{{\alpha }_2}<\left|{\omega }_1^{{\alpha }_1}-{\omega }_2^{{\alpha }_2}\right|$, respectively.

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 $A^{{\alpha }_1}+B^{{\alpha }_2}>\left|{\omega }_1^{{\alpha }_1}-{\omega }_2^{{\alpha }_2}\right|$ what we conjecture that indicates the synchronization.

Figure 7
Solutions θ ₁ and θ ₂ 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 $A^{{\alpha }_1}+B^{{\alpha }_2}<\left|{\omega }_1^{{\alpha }_1}-{\omega }_2^{{\alpha }_2}\right|$ what we conjecture that indicates the desynchronization.

In particular, Figure 6 numerically verifies the conclusions on Proposition 4.1.

In the remained numerical example presented in Figure 8a-d, we test the conjecture presented just after Proposition 4.1, regarding the asymptotic convergence of a possible synchronized solution for the multi-order fractional phase model towards the up and lower bound solutions, as α ₁→α ₂. For that fate, we choose the coupling forces and intrinsic frequencies corresponding the values in the fifth line on Table 1. Then, we plot the corresponding solutions θ _i , θ _iL and θ _iU , for i=1, 2, for α ₁=0, 48; 0, 45; 0, 42; and 0, 41, while α ₂=0, 4, respectively in Figure 8a-d. As it can be see in Figure 8, the conjecture is numerically verified.

Figure 8
Asymptotic behaviour of the solutions θ _i , θ _iL and θ _iU of the multi-order fractional phase model as α ₁→α ₂.

6 CONCLUSIONS AND FUTURE WORKS

In this contribution, we have extended the coupling oscillators phase model proposed by Strogatz¹⁴ 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.

• Analyse the proposed model using other versions of fractional derivatives, e.g.,^13),(19),(20.

• 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}.

REFERENCES

• ¹
  E. Bard, P. Youngmin & W. Dan. Recent advances in coupled oscillator theory. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 377(2160) (2019), 20190092. doi: 10.1098/rsta.2019.0092.
  » https://doi.org/10.1098/rsta.2019.0092
• ²
  C. Bick, G. Goodfellow, L.C. R. & E.A. Martens. Understanding the dynamics of biological and neural oscillator networks through mean-field reductions: a review (2019).
• ³
  N. Contessa & A. De Cezaro. Derivadas fracionárias na modelagem do ritmo circadiano. Scientia Plena, 13(4) (2017), 0499091-11.
• ⁴
  S. Da Silveira Glaeser, A.T. Dos Santos, A. De Cezaro, C.M.D.S. Machado & D.F. Adamatti. Modeling of Circadian Rhythm under influence of Pain: an approach based on Multi-agent Simulation. ADCAIJ: Advances in Distributed Computing and Artificial Intelligence Journal, 7(2) (2018), 17-25.
• ⁵
  K. Diethelm. “The analysis of fractional differential equations”. Lectutes Notes in Mathematics. Springer, Diethelm (2004).
• ⁶
  K. Diethelm. The mean value theorems and a Nagumo-type uniqueness theorem for Caputo’s fractional calculus. Fractional Calculus and Applied Analysis, 15(2) (2012), 304-313.
• ⁷
  R. Garrappa. Stability-preserving high-order methods for multiterm fractional differential equations. International Journal of Bifurcation and Chaos, 22(04) (2012), 1250073.
• ⁸
  E.B. Klerman & M.S. Hilaire. Review: On Mathematical Modeling of Circadian Rhythms, Performance, and Alertness. Journal of Biological Rhythms, 22(2) (2007), 91-102.
• ⁹
  Y. Kuramoto. “Chemical Oscillations, Waves, and Turbulence”. Springer, Berlin (1984).
• ¹⁰
  M.C. Moore-Ede, C.A. Czeisler & G.S. Richardson. Circadian Timekeeping in Health and Disease. New England Journal of Medicine, 309(9) (1983), 530-536.
• ¹¹
  A. Pikovisky, M. Rosemblum & J. Kurths. “Synchronization: A Universal Concept in Nonlinear Sciences”. Cambridge University Press, Cambridge (2001).
• ¹²
  F.A. Rodrigues, T.K.D. Peron, P. Ji & J. Kurths. The Kuramoto model in complex networks. Physics Reports, 610(Supplement C) (2016), 1-98.
• ¹³
  J. Sousa, G. Frederico & E. Oliveira. ψ-Hilfer pseudo-fractional operator: new results about fractional calculus. Comp. Appl. Math., 39(254) (2020), 020-01304-6.
• ¹⁴
  S.H. Strogatz. Human sleep and circadian rhythms: a simple model based on two coupled oscillators. J. Math. Biol, 25 (1987), 327-347.
• ¹⁵
  S.H. Strogatz. From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators. Physica D., 143 (2000), 1-20.
• ¹⁶
  P.A. Tass. “Phase Resetting in Medicine and Biology: Stochastic Modelling and Data Analysis”. Springer, NY (1999).
• ¹⁷
  M.J. Thorpy, E. Korman, A.J. Spielman & P.B. Glovinsky. Delayed sleep phase syndrome in adolescents. Journal of Adolescent Health Care, 9(1) (1988), 22-27.
• ¹⁸
  C. Tisdell. Solutions to Fractional Diffetential Equations that Extend. Journal of Classical Analysis, 5(2) (2014), 129-136.
• ¹⁹
  J. Vanterler da C. Sousa & E. Capelas de Oliveira. On the ψ-Hilfer fractional derivative. Communications in Nonlinear Science and Numerical Simulation, 60 (2018), 72-91.
• ²⁰
  J. Vanterler da C. Sousa & E. Capelas de Oliveira. Leibniz type rule: ψ-Hilfer fractional operator. Communications in Nonlinear Science and Numerical Simulation, 77 (2019), 305-311.
• ²¹
  Z. Wang, D. Yang & T.E.A. Ma. Stability analysis for nonlinear fractional-order systems based on comparison principle. Nonlinear Dyn., (75) (2014), 387-402.
• ²²
  A. Winfree. “The Geometry of Biological Time”. Springer-Verlag, NY (2001).
• ²³
  M. Zhang, G.S. Wiederhecker, S. Manipatruni, A. Barnard, P. McEuen & M. Lipson. Synchronization of micromechanical oscillators using light. Phys. Rev. Lett., (109) (2012), 233906.

• 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⁵. The analysis of the model in termos of ψ-Hilfer pseudo-fractional operator^{13), (19), (20} will be investigate in future contributions.
• 3
  The Winfree restrictions correspond to the derived Winfree restrictions derived for the parameter values in each of the phase models analysed. For the multi-order fractional phase model, those restrictions are unavailable.

Publication Dates

• Publication in this collection
  29 May 2023
• Date of issue
  Apr-Jun 2023

History

• Received
  07 Mar 2021
• Accepted
  18 Sept 2022