Problem #1

library(deSolve)
library(ggplot2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(see)
library(DT)

# Parameters for Muscle neurons
parameters <- c(Vna = 115, Vk = -11, Vca = 115, Vl = 10.3, gna = 110, 
                gk = 36, gca = 50, gl= 0.25, C = 1, I = 10)

# Initial conditions for Muscle neurons
yini_Muscle <- c(V = -14, m = 0.053, h = 0.68, n = 0.318, p=0.0)

# HH model equations for Muscle neurons
modelHH <- function(time, y, parms) {
  with(as.list(c(y, parms)), {
    alpha_m <- function(v) { 0.1*(25-v)/(exp((25-v)/10)-1) }
    beta_m  <- function(v) { 4*exp(-v/18) }
    alpha_h <- function(v) { 0.07*exp(-v/20) }
    beta_h  <- function(v) { 1/(exp((30-v)/10)+1) }
    alpha_n <- function(v) { 0.01*(10-v)/(exp((10-v)/10)-1) }
    beta_n  <- function(v) { 0.125*exp(-v/80) }
    
    alpha_p <- function(v) { 0.01*(15-v)/(exp((15-v)/15)-1) }
    beta_p  <- function(v) { 0.1*exp(-v/10) }

    dp <- alpha_p(V)*(1-p)-beta_p(V)*p 
    dp <- alpha_p(V)*(1-p)-beta_p(V)*p 
    
    dV <- (I - gna*h*(V-Vna)*m^3-gk*(V-Vk)*n^4-gl*(V-Vl))/C
    dm <- alpha_m(V)*(1-m)-beta_m(V)*m
    dh <- alpha_h(V)*(1-h)-beta_h(V)*h
    dn <- alpha_n(V)*(1-n)-beta_n(V)*n 
    return(list(c(dV, dm, dh, dn, dp)))
  })
}

time = seq(from =0, to=100, by = 0.01)

print(system.time(
  out_modelHH_Muscle <- 
    ode(yini_Muscle, func=modelHH, times=time, parms=parameters)))
##    user  system elapsed 
##   0.262   0.004   0.267
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
modelHH_Muscle_df <- as.data.frame(out_modelHH_Muscle)

plot_ly(data = modelHH_Muscle_df, type="scatter", mode="lines") %>%
  add_trace(x=~time, y=~V/100, name="V/100 (Membrane Action Potential)") %>%
  layout(title="Muscle Neurons: ODE Solutions of the Hodgkin-Huxley Model",
         legend = list(orientation = 'h'))

Variations in ion channel probabilities are key determinants of the electrophysiological properties of living cells. For instance, increased activation of sodium channels facilitates rapid depolarization, potentially accelerating action potential initiation and promoting rapid signaling in neurons and muscle cells. On the other hand, altered probabilities of potassium channel activation can modify the rate of repolarization and the neuron’s readiness for subsequent firing, influencing rhythmic firing and overall cellular excitability. These intricacies in ion channel behavior are essential for the proper functioning of cells and can impact critical activities, from muscle contraction coordination to the modulation of neural circuits involved in cognitive functions.

Problem #2

library(deSolve)
library(ggplot2)
library(dplyr)
library(see)
library(DT)

# Parameters for Muscle neurons
parameters <- c(Vna = 115, Vk = -11, Vca = 115, Vl = 10.3, gna = 110, 
                gk = 36, gca = 50, gl= 0.25, C = 1, I = 10)

# Initial conditions
yini_Muscle <- c(V = -14, m = 0.1, h = 0.7, n = 0.4, p=0.0)
yini_Dopamin <- c(V = -14, m = 0.3, h = 0.6, n = 0.5, p=0.5)

# HH model equations for Muscle neurons
modelHH <- function(time, y, parms) {
  with(as.list(c(y, parms)), {
    alpha_m <- function(v) { 0.1*(25-v)/(exp((25-v)/10)-1) }
    beta_m  <- function(v) { 4*exp(-v/18) }
    alpha_h <- function(v) { 0.07*exp(-v/20) }
    beta_h  <- function(v) { 1/(exp((30-v)/10)+1) }
    alpha_n <- function(v) { 0.01*(10-v)/(exp((10-v)/10)-1) }
    beta_n  <- function(v) { 0.125*exp(-v/80) }
    
    alpha_p <- function(v) { 0.01*(15-v)/(exp((15-v)/15)-1) }
    beta_p  <- function(v) { 0.1*exp(-v/10) }

    dp <- alpha_p(V)*(1-p)-beta_p(V)*p 
    dp <- alpha_p(V)*(1-p)-beta_p(V)*p 
    
    dV <- (I - gna*h*(V-Vna)*m^3-gk*(V-Vk)*n^4-gl*(V-Vl))/C
    dm <- alpha_m(V)*(1-m)-beta_m(V)*m
    dh <- alpha_h(V)*(1-h)-beta_h(V)*h
    dn <- alpha_n(V)*(1-n)-beta_n(V)*n 
    return(list(c(dV, dm, dh, dn, dp)))
  })
}

time = seq(from =0, to=100, by = 0.01)

print(system.time(
  out_modelHH_Muscle <- 
    ode(yini_Muscle, func=modelHH, times=time, parms=parameters)))
##    user  system elapsed 
##   0.248   0.007   0.254
print(system.time(
  out_modelHH_Dopamin <- 
    ode(yini_Dopamin, func=modelHH, times=time, parms=parameters)))
##    user  system elapsed 
##   0.238   0.004   0.243
library(plotly)
modelHH_Muscle_df <- as.data.frame(out_modelHH_Muscle)

modelHH_Dopamin_df <- as.data.frame(out_modelHH_Dopamin)

ggplot() +
  geom_line(data = modelHH_Muscle_df, aes(x = time, y = V), color = "blue",) +
  geom_line(data = modelHH_Dopamin_df, aes(x = time, y = V), color = "red",) +
  labs(title = "Action Potential Comparison For Dopaminergic and Muscle Neurons", x = "Time (ms)", y = "Membrane Potential (mV)") +
  theme_minimal() +
  scale_color_manual(values = c("Muscle" = "blue", "Dopaminergic" = "red"))

Spike distribution differences between dopaminergic neurons (red) and muscle neurons (blue) are largely attributable to the presence of calcium (Ca²⁺) channels in dopaminergic neurons and the resulting distinct electrophysiological characteristics. Dopaminergic neurons may demonstrate a longer-duration action potential with a more pronounced plateau phase due to the Ca²⁺ influx, which is crucial in neurotransmitter release and neuronal signaling. Muscle neurons, which typically lack this plateau phase in their action potentials, are designed for rapid and short bursts of activation, essential for quick and forceful muscle contractions. Consequently, the inclusion of Ca²⁺ channels in dopaminergic neurons contributes to their unique firing patterns compared to muscle neurons, reflecting their specialized roles within the nervous system and the neuromuscular junctions, respectively.

Problem #3

library(deSolve)
library(ggplot2)
library(dplyr)
library(see)
library(DT)

# Parameters for Muscle neurons
parameters <- c(Vna = 115, Vk = -11, Vca = 115, Vl = 10.3, gna = 110, 
                gk = 36, gca = 50, gl= 0.25, C = 1, I = 10)

# Initial conditions
yini_Dopamin <- c(V = -14, m = 0.3, h = 0.6, n = 0.5, p=0.0)
yini_DopaminCa <- c(V = -14, m = 0.3, h = 0.6, n = 0.5, p=0.5)

# HH model equations for Muscle neurons
modelHH <- function(time, y, parms) {
  with(as.list(c(y, parms)), {
    alpha_m <- function(v) { 0.1*(25-v)/(exp((25-v)/10)-1) }
    beta_m  <- function(v) { 4*exp(-v/18) }
    alpha_h <- function(v) { 0.07*exp(-v/20) }
    beta_h  <- function(v) { 1/(exp((30-v)/10)+1) }
    alpha_n <- function(v) { 0.01*(10-v)/(exp((10-v)/10)-1) }
    beta_n  <- function(v) { 0.125*exp(-v/80) }
    
    alpha_p <- function(v) { 0.01*(15-v)/(exp((15-v)/15)-1) }
    beta_p  <- function(v) { 0.1*exp(-v/10) }

    dp <- alpha_p(V)*(1-p)-beta_p(V)*p 
    dp <- alpha_p(V)*(1-p)-beta_p(V)*p 
    
    dV <- (I - gna*h*(V-Vna)*m^3-gk*(V-Vk)*n^4-gl*(V-Vl))/C
    dm <- alpha_m(V)*(1-m)-beta_m(V)*m
    dh <- alpha_h(V)*(1-h)-beta_h(V)*h
    dn <- alpha_n(V)*(1-n)-beta_n(V)*n 
    return(list(c(dV, dm, dh, dn, dp)))
  })
}

time = seq(from =0, to=100, by = 0.01)

print(system.time(
  out_modelHH_Dopamin <- 
    ode(yini_Dopamin, func=modelHH, times=time, parms=parameters)))
##    user  system elapsed 
##   0.254   0.002   0.255
print(system.time(
  out_modelHH_DopaminCa <- 
    ode(yini_DopaminCa, func=modelHH, times=time, parms=parameters)))
##    user  system elapsed 
##   0.238   0.002   0.240
library(plotly)
modelHH_Dopamin_df <- as.data.frame(out_modelHH_Dopamin)

modelHH_DopaminCa_df <- as.data.frame(out_modelHH_DopaminCa)

ggplot() +
  geom_line(data = modelHH_Dopamin_df, aes(x = time, y = V), color = "blue",) +
  geom_line(data = modelHH_DopaminCa_df, aes(x = time, y = V), color = "red",) +
  labs(title = "Dopaminergic Neuron Spike Distribution With and Without Ca Channels", x = "Time (ms)", y = "Membrane Potential (mV)") +
  theme_minimal() +
  scale_color_manual(values = c("Dopaminergenic without Ca" = "blue", "Dopaminergic with Ca" = "red"))

The presence of calcium (Ca₂⁺) channels typically has a marked effect on the action potential profile of dopaminergic neurons. When calcium channels are considered (red), the spike distribution often shows a plateau phase due to the extended Ca₂⁺ influx, leading to prolonged depolarization. This contrasts with the distribution observed without calcium channels (green), where the action potential has a sharper peak and quicker repolarization, resulting in a more traditional, brief spike shape. These disparities highlight the crucial role of Ca₂⁺ in modulating the rate and duration of action potentials, which can affect the neuron’s synaptic transmission and overall firing patterns.

Problem #4

library(deSolve)
library(ggplot2)
library(dplyr)
library(see)
library(DT)

# Parameters for Muscle neurons
params_a <- c(Vna = 115, Vk = -11, Vl = 10.3, gna = 110, gk = 36, gl = 0.25, C = 1, I = 10)
params_new <- c(C = 1, gl = 0.3, Vl = -54.4, gk = 36, Vk = -77, gna = 120, Vna = 50, I = 10, gca = 5, Vca = 110)

# Initial conditions
yini_a <- c(V = -14, m = 0.1, h = 0.7, n = 0.4, p = 0)
yini_new <- c(V = -65, m = 0.2, h = 0.7, n = 0.3, p = 0)

# HH model equations for Muscle neurons
modelHH <- function(time, y, parms) {
  with(as.list(c(y, parms)), {
    alpha_m <- function(v) { 0.1*(25-v)/(exp((25-v)/10)-1) }
    beta_m  <- function(v) { 4*exp(-v/18) }
    alpha_h <- function(v) { 0.07*exp(-v/20) }
    beta_h  <- function(v) { 1/(exp((30-v)/10)+1) }
    alpha_n <- function(v) { 0.01*(10-v)/(exp((10-v)/10)-1) }
    beta_n  <- function(v) { 0.125*exp(-v/80) }
    
    alpha_p <- function(v) { 0.01*(15-v)/(exp((15-v)/15)-1) }
    beta_p  <- function(v) { 0.1*exp(-v/10) }

    dp <- alpha_p(V)*(1-p)-beta_p(V)*p 
    dp <- alpha_p(V)*(1-p)-beta_p(V)*p 
    
    dV <- (I - gna*h*(V-Vna)*m^3-gk*(V-Vk)*n^4-gl*(V-Vl))/C
    dm <- alpha_m(V)*(1-m)-beta_m(V)*m
    dh <- alpha_h(V)*(1-h)-beta_h(V)*h
    dn <- alpha_n(V)*(1-n)-beta_n(V)*n 
    return(list(c(dV, dm, dh, dn, dp)))
  })
}

time = seq(from =0, to=100, by = 0.01)

print(system.time(
  out_modelHH_a <- 
    ode(yini_a, func=modelHH, times=time, parms=params_a)))
##    user  system elapsed 
##   0.244   0.002   0.246
print(system.time(
  out_modelHH_new <- 
    ode(yini_new, func=modelHH, times=time, parms=params_new)))
##    user  system elapsed 
##   0.222   0.003   0.225
library(plotly)
modelHH_a_df <- as.data.frame(out_modelHH_a)

modelHH_new_df <- as.data.frame(out_modelHH_new)

ggplot() +
  geom_line(data = modelHH_a_df, aes(x = time, y = V), color = "blue",) +
  geom_line(data = modelHH_new_df, aes(x = time, y = V), color = "red",) +
  labs(title = "Simulation of Muscle Neuron Spike Profiles With Varying Initial Conditions", x = "Time (ms)", y = "Membrane Potential (mV)") +
  theme_minimal() +
  scale_color_manual(values = c("Part a" = "blue", "New" = "red"))

The spike profile from part A (blue) accurately represents the typical behavior of muscle neurons, generating repeated action potentials necessary for muscle contractions. Conversely, the new set of initial conditions (red) fails to elicit spiking, likely indicating a membrane potential below the necessary threshold, which would result in muscle inactivity in a biological context. The observed differences underscore the crucial role of correct ion channel function and initial membrane potential in regulating muscle neuron excitability and, by extension, muscle function. An inability to fire action potentials means muscle fibers would not receive the signals required for contraction, highlighting the essential nature of these electrophysiological properties.