knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) library(sleepR)
The function philrob
simulates Phillips and Robinson's model for sleep-wake 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
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$ |
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 |
\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.
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
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 / somnogenasleep
: the asleep/awake status (TRUE
if asleep, FALSE
if awake)philrobPlot(sol) rasterPlot(sol)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.