---
title: "The Parameters of a \n Dark Adaptation Model Explained"
author: "Dr Jeremiah MF Kelly"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Parameters Explained}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

The following plots explain the parameters of the model used to predict thresholds ($Thrs$) following a photo-bleach which transforms more than 30% of the retinal photopigment. 
These data were extracted from Figure 1 in Pugh 1975. 

There are typically three phases to the recovery of visual sensitivity in dark adaptation data; an early cone phase, the second rod phase sometimes called S2, and a final rod phase or S3 and are shown in the figure below. 

The initial phase depicts the cone response ($Thrs_c$), and can be represented by three parameters; the absolute cone threshold ($\theta_1$, log$_{10}$(cd.m$^{\text{-}2}$) scotopic), the cone coefficient  ($\theta_2$, log$_{10}$(cd.m$^{\text{-}2}$) scotopic), and the time constant of cone recovery, sometimes known as tau ($\theta_3$, min).

$$
Thrs_c = \theta_1 + \theta_2.\exp\left(\frac{\text{-}t}{\theta_3}\right) .
$$

The cone coefficient can be thought of as the cone threshold at the instant following the flash, but as this variable can take non-physiological values it is best perhaps thought of a free variable that aids the fitting process.


```{r, fig.align='center', fig.width=8, fig.height=8, echo=FALSE}
load("dark.rda")
tmp<-dark
  par(las = 1, bty = "n", 
      mfrow = c(1, 1), 
      font = 2, font.lab = 2, 
      cex.lab = 1.5, 
      oma = c(1, 1, 1, 1), 
      mar = c(6, 6, 2, 3)
      )
  XL <- expression(bold(Time ~ (min)))
  YL <- expression(bold(Threshold ~ (log ~ units)))
  
  plot(tmp$time, tmp$thrs, xlab = XL, ylab = YL, pch = 1, cex = 1)
  axis(1, lwd = 2)
  axis(2, lwd = 2)
 
  lines(tmp$time, tmp$fit, col = 2)
  CT <- tmp$opt[1]
  CC <- tmp$opt[2]
  Tau <- tmp$opt[3]
  S2 <- tmp$opt[4]
  Al <- tmp$opt[5]
  S3 <- tmp$opt[6]
  Be <- tmp$opt[7]
  
  
 
  rect(0, CT, Al, max(tmp$thrs), border = NA, col = rgb(1, 1, 0, 0.2))
idx<-tmp$time<Al
points(tmp$time[idx], tmp$thrs[idx], pch=16, cex=0.7)

Off <- c(1.5, 3) * -0.2
text(Al, max(tmp$thrs), "Cone", adj = Off)

Lty = 4

Sym <- expression(theta[1] ~ Cone ~ Threshold)
arrows(Be + 10, CT, -2, CT, length = 0.125, angle = 20, lty = Lty)
text(Be + 10, CT, Sym, adj = c(0, 0), cex = 1.2)



Sym <- expression(theta[2] ~ Cone ~ Coefficient)
arrows(Be + 10, CT + CC, -2, CC + CT, length = 0.125, angle = 20, lty = Lty)
text(Be + 10, CC + CT, Sym, adj = c(0, 0), cex = 1.2)


Sym <- expression(theta[3] ~ Cone ~ Time ~ Constant ~ (tau))
# arrows(Be+10,CT+CC/2, 3, CT+CC/2, length = 0.125, angle = 20, lty=Lty)
text(3, CT + CC/2, Sym, adj = c(0, 0), cex = 1.2)
```

After the cone mediated phase there follows the rod phase, S2, which represents the rate limited restoration of photopigment and is defined by its slope, ($\theta_4$, log$_{10}$(cd.m$^{\text{-}2}$).min$^{\text{-}1}$) and the rod-cone transition time or alpha point ($\theta_5$, min). It is clear from the figure that the alpha point is defined by the slope and the intercept of the S2 phase.

The linear (on a log-lin plot) rod phase S2 is described in detail in Lamb and Pugh 2004.


```{r, fig.align='center', fig.width=8, fig.height=8, echo=FALSE}
source('../R/H.R')
source('../R/P7c.R')
tmp<-dark
  par(las = 1, bty = "n", 
      mfrow = c(1, 1), 
      font = 2, font.lab = 2, 
      cex.lab = 1.5, 
      oma = c(1, 1, 1, 1), 
      mar = c(6, 6, 2, 3)
      )
  XL <- expression(bold(Time ~ (min)))
  YL <- expression(bold(Threshold ~ (log ~ units)))
  
  plot(tmp$time, tmp$thrs, xlab = XL, ylab = YL, pch = 1, cex = 1)
  axis(1, lwd = 2)
  axis(2, lwd = 2)
 
  lines(tmp$time, tmp$fit, col = 2)
  CT <- tmp$opt[1]
  CC <- tmp$opt[2]
  Tau <- tmp$opt[3]
  S2 <- tmp$opt[4]
  Al <- tmp$opt[5]
  S3 <- tmp$opt[6]
  Be <- tmp$opt[7]
  
  
 
rect(Al, P7c(tmp$opt, Be), Be, P7c(tmp$opt, Al), border = NA, col = rgb(1, 0, 0, 0.2))

idx<-tmp$time>Al & tmp$time<Be

points(tmp$time[idx], tmp$thrs[idx], pch=16, cex=0.7)

mod<-lm(tmp$thrs[idx]~tmp$time[idx])
abline(mod, col=2, lty=2)
Off <- c(1, 1) * 3 * -0.08
text(Be, P7c(tmp$opt, Al), "S2", adj = Off)

Lty = 4
Sym <- expression(atop(theta[4] ~ Second ~ Rod ~ Recovery, ~Rate ~ (S2)))
# arrows(Be+10,CT+CC/2, 3, CT+CC/2, length = 0.125, angle = 20, lty=Lty)
text(Be - 3, mean(tmp$thrs) - 0.2, Sym, adj = c(0, 0), cex = 1.2)

Sym <- expression(atop(theta[5] ~ Cone ~ Rod, ~Breakpoint ~ (alpha)))
arrows(Al - 3, P7c(tmp$opt, Al) - 1.4, Al, P7c(tmp$opt, Al), length = 0.125, angle = 20, 
       lty = Lty)
text(Al - 5, P7c(tmp$opt, Al) - 1.8, Sym, adj = c(0.5, 0), cex = 1.2)
```

Similarly the next rod phase, is described by the sum of the slopes of S2 and S3 ($\theta_6$, log$_{10}$(cd.m$^{\text{-}2}$).min$^{\text{-}1}$) and the rod-rod transition time, occasionally referred to as the beta point ($\theta_7$, min).

The final phase of rod recovery S3 may not be found in all subjects, e.g. a subject with reduced rod function may take an hour or more to reach this final phase. 


```{r, fig.align='center', fig.width=8, fig.height=8, echo=FALSE}
library(Dark)
data(dark)
tmp<-dark
  par(las = 1, bty = "n", 
      mfrow = c(1, 1), 
      font = 2, font.lab = 2, 
      cex.lab = 1.5, 
      oma = c(1, 1, 1, 1), 
      mar = c(6, 6, 2, 3)
      )
  XL <- expression(bold(Time ~ (min)))
  YL <- expression(bold(Threshold ~ (log ~ units)))
  
  plot(tmp$time, tmp$thrs, xlab = XL, ylab = YL, pch = 1, cex = 1)
  axis(1, lwd = 2)
  axis(2, lwd = 2)
 
  lines(tmp$time, tmp$fit, col = 2)
  CT <- tmp$opt[1]
  CC <- tmp$opt[2]
  Tau <- tmp$opt[3]
  S2 <- tmp$opt[4]
  Al <- tmp$opt[5]
  S3 <- tmp$opt[6]
  Be <- tmp$opt[7]
rect(Be, min(tmp$thrs), max(tmp$time), P7c(tmp$opt, Be), border = NA, col = rgb(1, 0, 
                                                                                1, 0.2))

Off <- c(1, 1) * 3 * -0.08
text(max(tmp$time), P7c(tmp$opt, Be), "S3", Off)  
idx<-tmp$time>Be

points(tmp$time[idx], tmp$thrs[idx], pch=16, cex=0.7)

mod<-lm(tmp$thrs[idx]~tmp$time[idx])
abline(mod, col=2, lty=2)

Lty = 4
Sym <- expression(atop(theta[6] ~ Third ~ Rod ~ Recovery, ~Rate ~ (-(S2 + S3))))
# arrows(Be+10,CT+CC/2, 3, CT+CC/2, length = 0.125, angle = 20, lty=Lty)
text(Be + 12, P7c(tmp$opt, Be) - 0.5, Sym, adj = c(0, 0), cex = 1.2)



Sym <- expression(atop(theta[7] ~ Rod ~ Rod, ~Breakpoint ~ (beta)))
arrows(Al, P7c(tmp$opt, Be) - 0.4, Be, P7c(tmp$opt, Be), length = 0.125, angle = 20, lty = Lty)
text(Al - 5, P7c(tmp$opt, Be) - 0.6, Sym, adj = c(0.5, 0), cex = 1.2)

```

The model uses the transition times to control switch functions and hence allow each phase to combine.

The response threshold $Thrs$ can therefore be summarized as follows:
$$
Thrs = Thrs_c +
\underbrace{ \theta_4.h(t,\theta_5) }_{\mbox{S2}}+\underbrace{ \theta_6.h(t,\theta_7) }_{\mbox{S3}}
$$
where $h(t,\theta_i)$ is a switch function:
$$
h(t,\theta_i)
\begin{cases} 0 & \text{if $t-\theta_i\leq0$,}
\\
t-\theta_i &\text{if $t- \theta_i> 0$.}
\end{cases}
\qquad i \in (5,7).
$$

The models in this package use a logistic function 
$$
h(t,\theta_i)\rightarrow \frac{1}{1+exp(-k(t-\theta_i))}.(t-\theta_i)\qquad i \in (5,7).
$$

This continuous switch function is computationally more efficient and allows an analytic approach to estimates of parameter variance. 

Notice that since the threshold data are in log units, when plotted against time the apparently linear phases S2 and S3 describe exponential processes.


# References
E. Pugh. Rushton's paradox: rod dark adaptation after flash photolysis. The Journal of Physiology, 248(2):413, 1975.

T. Lamb and E. Pugh. Dark adaptation and the retinoid cycle of vision. Progress in Retinal and Eye Research, 23(3):307–380, 2004.