In PabRod/sleepR: Sleep-wake Cycle Models

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

library(sleepR)

Model

The function philrob simulates Phillips and Robinson's model for sleep-wake dynamics.

Dynamics

The dynamics of the Phillips and Robinson model are given by the following system of ordinary differential equations:

$$\left[ \begin{array}{c} \tau_v \dot V_v + V_v \ \tau_m \dot V_m + V_m \ \chi \dot H + H \end{array} \right] = \left[ \begin{array}{ccc} 0 & -\nu_{vm} & \nu_{vh} \ -\nu_{mv} & 0 & 0 \ 0 & \mu & 0 \end{array} \right] \left[ \begin{array}{c} S(V_v) \ S(V_m) \ H \end{array} \right] + \left[ \begin{array}{c} -\nu_{vc} C(t) \ \nu_{ma} S(V_{a0}) \ 0 \end{array} \right]$$

where $S(V)$ is the saturation function:

$$S(V) = \frac{Q_{max}}{1 + e^{-\frac{V-\theta}{\sigma}}}$$

parms <- philrob_default_parms()
S <- function(V) {
  parms['Qmax']/(1 + exp((parms['theta']-V)/parms['sigma']))
}

xs <- seq(-10, 35, length.out = 100)
plot(xs, S(xs), type = 'l', xlab = 'V', ylab = 'S(V)', col = 'red', axes = FALSE)

# Asymptote
lines(c(-10, 35), c(parms['Qmax'], parms['Qmax']), lty = 'dashed')

# Half life
lines(c(-10, parms['theta']), c(parms['Qmax']/2, parms['Qmax']/2), lty = 'dashed')
lines(c(parms['theta'], parms['theta']), c(0, parms['Qmax']/2), lty = 'dashed')

axis(2, at = c(0, parms['Qmax']/2, parms['Qmax']), labels = c('0', 'Qmax/2', 'Qmax'))
axis(1, at = c(-10, 0, parms['theta'], 35), labels = c('', '0', 'theta', ''))

and the external forcing is typically given by $C(t)$:

$$C(t) = \frac{1}{2} \left( 1 + \cos(\omega t + \alpha) \right)$$

\newpage

Parameters

The default values for the parameters are listed in the table below:

| Symbol | Value | Units | |:--------------------|:----------------:|:-------------------| | $\tau_m$ | $10/3600$ | $h$ | | $\tau_v$ | $10/3600$ | $h$ | | $\chi$ | $10.8$ | $h$ | | $\nu_{vm}$ | $1.9/3600$ | $mV \cdot h$ | | $\nu_{mv}$ | $1.9/3600$ | $mV \cdot h$ | | $\nu_{vh}$ | $0.19$ | $mV \cdot nM^{-1}$ | | $\mu$ | $10^{-3}$ | $nM \cdot h$ | | $\nu_{vc}$ | $6.3$ | $mV$ | | $\nu_{ma}S(V_{a0})$ | $1$ | $mV$ | | $Q_{max}$ | $100 \cdot 3600$ | $h^{-1}$ | | $\theta$ | $10$ | $mV$ | | $\sigma$ | $3$ | $mV$ | | $\omega$ | $2 \pi / 24$ | $h^{-1}$ | | $\alpha$ | $0$ | $1$ |

State variables

The state variables are defined as:

| State variable | Units | Physiological interpetation | Informal interpretation | |:---------------|:-----:|:----------------------------|:------------------------| | $V_v$ | $mV$ | Activity of the VLPO | Stay asleep system | | $V_m$ | $mV$ | Activity of the MA | Stay awake system | | $H$ | $1$ | Homeostatic pressure | Somnogen level |

Diagram

\begin{center} \begin{tikzpicture}[node distance = 2cm, auto]

% Definitions
\tikzstyle{state} = [rectangle,draw,fill=blue!20,text width=5em,text centered,rounded corners,minimum height=4em]
\tikzstyle{line} = [draw, -latex']
\tikzstyle{source} = [draw, ellipse,fill=red!20, node distance=3cm,minimum height=2em]

% Place nodes
% States
\node [state] (Vv) {$V_v$};
\node [state, below of=Vv] (Vm) {$V_m$};
\node [state, below of=Vm] (H) {$H$};

% Sources
\node [source, right of=Vv] (C) {$C(t)$};
\node [source, right of=Vm] (A) {$A$};

% Draw edges
% Sources
\path [blue, line] (C) -> node{$\nu_{vc}$} (Vv);
\path [line] (A) -> node{$\nu_{ma}$} (Vm);

% Coupling
\path [line] (Vm) -> node{$\mu$} (H);
\draw[red, ->] (Vv) -- node[midway, right] {$\nu_{mv}$} (Vm);
\draw[red, ->] (Vm) -- node[midway, left] {$\nu_{vm}$} (Vv);
\draw[->] (H.west) .. controls +(left:20mm) and +(left:20mm) .. node{$\nu_{vh}$} (Vv.west);

% Decay
\draw[red, ->] (Vv.north) .. controls +(up:10mm) and +(up:10mm) .. (Vv.north west);
\draw[red, ->] (Vm.north west) .. controls +(left:10mm) and +(left:10mm) .. (Vm.south west);
\draw[red, ->] (H.south) .. controls +(down:10mm) and +(down:10mm) .. (H.south west);

\end{tikzpicture} \end{center}

Schematic summary of the dynamics. The blue nodes represent the system's states ($V_v$ the activity of the ventrolateral preoptic area, $V_m$ the activity of the mono aminergic group and $H$ the homeostatic pressure). The red nodes represent the external sources ($C(t)$, the astronomical light/dark forcing, and $A$, the acetylcholine group constant influence). The positive effects are coded as black arrows. Negative ones as red arrows. Blue arrows represent oscillating effects.

Reference

Phillips AJK, Robinson PA. A Quantitative Model of Sleep-Wake Dynamics Based on the Physiology of the Brainstem Ascending Arousal System. J Biol Rhythms. 2007 Apr 29;22(2):167–79. Available from: http://journals.sagepub.com/doi/10.1177/0748730406297512

\newpage

Examples of usage

Getting the time series

With default parameters:

## Problem setting
y0 <- c(Vv = -13, Vm = 1, H = 10) # Initial conditions

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

# Simulate
sol <- philrob(ts, y0)

With custom parameters:

## Problem setting
y0 <- c(Vv = -13, Vm = 1, H = 10) # Initial conditions

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

parms <- philrob_default_parms() # Load default parameters...
parms['vvc'] <- 6 # .. and modify one

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

With custom forcing:

## Problem setting
y0 <- c(Vv = -13, Vm = 1, H = 10) # Initial conditions

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

C <- function(t) { 0 }

# Simulate
sol <- philrob(ts, y0, parms = philrob_default_parms(), forcing = C)

With stabilization run of three days:

## Problem setting
y0 <- c(Vv = -13, Vm = 1, H = 10) # Initial conditions

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

# Simulate
sol <- philrob(ts, y0, tStabil = 3*24)

The output looks like:

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

where:

• time: the time (in h),
• Vv: activity of the ventrolateral preoptic area (in mV)
• Vm: activity of the monoaminergic group (in mV)
• H: homeostatic pressure / somnogen
• asleep: the asleep/awake status (TRUE if asleep, FALSE if awake)

Plotting results

Raster / somnogram plot

philrobPlot(sol)
rasterPlot(sol)

PabRod/sleepR documentation built on Dec. 21, 2020, 3:27 p.m.