In PabRod/sleepR: Sleep-wake Cycle Models

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(sleepR)

Model

The function strogatz simulates Strogatz's model for sleep-wake dynamics. The dynamics have the following structure:

$$ \dot \theta_1 = \omega_1 - C_1 \cdot \cos(2 \pi (\theta_2 - \theta_1)) \ \dot \theta_2 = \omega_2 + C_2 \cdot \cos(2 \pi (\theta_2 - \theta_1)) $$

Analysis

It is useful to note that, using the relative phase $\psi \equiv \theta_1 - \theta_2$, the system can be simplified to:

$$ \dot \psi = \Omega - C \cos(2 \pi \psi) $$ where $\Omega \equiv \omega_1 - \omega_2$ and $C \equiv C_1 + C_2$. The relative phase can be stabilized if and only if $\dot \psi = 0$, and this is only possible if:

$$ \lvert \frac{\Omega}{C} \rvert \leq 1 $$

Reference

Strogatz, S. H. (1987). Human sleep and circadian rhythms: a simple model based on two coupled oscillators. Journal of Mathematical Biology, 25(3), 327–347. http://doi.org/10.1007/BF00276440

\newpage

Examples of usage

Getting the time series

With default parameters:

## Problem setting
y0 <- c(th1 = 0.1, th2 = 0.05) # Initial conditions

nDays <- 5
ts <- seq(0, nDays*24, length.out=nDays*24*20) # Times to simulate

# Simulate
sol <- strogatz(ts, y0)

With custom parameters:

## Problem setting
y0 <- c(th1 = 0.1, th2 = 0.05) # Initial conditions

nDays <- 5
ts <- seq(0, nDays*24, length.out=nDays*24*20) # Times to simulate
parms <- c(w1=1/24, w2=0.85/24, C1=0/24, C2=0.16/24) # Parameters

# Simulate
sol <- strogatz(ts, y0, parms)

The output looks like:

knitr::kable(head(sol, 5))

where:

time: the time (in h),
th1: the phase of the circadian oscillator
th2: the phase of the sleep-wake cycle
asleep: the asleep/awake status (TRUE if asleep, FALSE if awake)

Plotting results

Raster / somnogram plot

plot(sol$time/24, sol$asleep, type = 'line', xlab = 'Time (d)', ylab = 'Asleep')
rasterPlot(sol)

Results

Entrained vs. non entrained case

## Problem setting
y0 <- c(th1 = 0.1, th2 = 0.05) # Initial conditions

nDays <- 60
ts <- seq(0, nDays*24, length.out=nDays*24*20) # Times to simulate

# Parameters
# With this settings, the bifurcation happens at w2 = 0.84/24 h^-1
parms_entrained <- c(w1=1/24, w2=0.845/24, C1=0/24, C2=0.16/24)
parms_non_entrained <- c(w1=1/24, w2=0.835/24, C1=0/24, C2=0.16/24)

# Simulate
sol_entrained <- strogatz(ts, y0, parms_entrained)
sol_non_entrained <- strogatz(ts, y0, parms_non_entrained)

# Plot
rasterPlot(sol_entrained)
title('Entrained')
rasterPlot(sol_non_entrained)
title('Not entrained')