---
title: "Hodgkin-Huxley"
author: "Daniel Kaschek, Mirjam Fehling-Kaschek"
date: "May 18, 2017"
output: html_document
runtime: shiny
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


Load some required packages

```{r}
library(deSolve)
library(tidyr)
library(ggplot2)
library(cowplot)
theme_set(ggthemes::theme_few())
```

Define the Hodgkin-Huxley differential equations

```{r}


alpha_n <- function(v) 0.01*(v + 55)/(1 - exp(-(v + 55)/10))
beta_n <- function(v)  0.125*exp(-(v + 65)/80)
alpha_m <- function(v) 0.1*(v+40)/( 1 - exp(-(v+40)/10) )
beta_m <- function(v) 4*exp(-(v+65)/18)
alpha_h <- function(v) 0.07*exp(-(v+65)/20)
beta_h <- function(v) 1/( 1 + exp(-(v+35)/10) )

hodge <- function(t, x, pars) {
  with(as.list(c(x, pars)), {
    dIapp <- 0
    dv <- (-g_K*n^4*(v-v_K) - g_Na*m^3*h*(v-v_Na) - g_L*(v-v_L) + Iapp)/C_m
    dn <- alpha_n(v)*(1-n) - beta_n(v)*n
    dm <- alpha_m(v)*(1-m) - beta_m(v)*m
    dh <- alpha_h(v)*(1-h) - beta_h(v)*h
    
    list(c(Iapp = dIapp, v = dv, n = dn, m = dm, h = dh))
    
  })
}

pars <- c(
  C_m = 1.0, #µF/cm²
  g_L = 0.3, #m.mho/cm²
  g_K = 36, #m.mho/cm²
  g_Na = 120, #m.mho/cm²
  v_K = -77, #mV
  v_Na = 50, #mV
  v_L = -54.4 #mV
)

ini <- c(
  Iapp = 0, 
  v = -60, # mV
  n = 0, 
  m = 0,
  h = 1
)


```

Plot function to create a plot for a specific value of Iapp.

```{r}

plotHH <- function(tmin = 40, tmax = 70, Iapp = 3) {
  
  times <- seq(0, 200, .1)
  events <- data.frame(
    var = "Iapp",
    time = c(50, 140),
    value = c(Iapp, 0),
    method = "rep"
  )
  
  out <- as.data.frame(ode(ini, times, hodge, pars, events = list(data = events)))
  out <- gather(out, "gate", "value", -time)
  out <- subset(out, time >= tmin & time <= tmax)
  
  P1 <- ggplot(subset(out, gate %in% c("Iapp", "v")), aes(x = time, y = value)) + facet_wrap(~gate, scales = "free", nrow = 2) + geom_line()
  
  P2 <- ggplot(subset(out, gate %in% c("n", "m", "h")), aes(x = time, y = value, group = gate, color = gate)) + geom_line() + theme(legend.position = "top")
  
  plot_grid(P1, P2, align = "v")
  
}


```


Implement the plot as interactive shiny plot

```{r eruptions, echo=FALSE}
inputPanel(
  sliderInput("tminmax", label = "Time interval:",
              min = 0, max = 200, value = c(45, 70), step = 1),
  sliderInput("Iapp", label = "Applied current Iapp:",
              min= -10, max = 50, value = 1, step = 1)
)

renderPlot({
  plotHH(input$tminmax[1], input$tminmax[2], input$Iapp)
})
```