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)
# Specify HH model parameters
## Ena, Ek, Eca, El, ionic reversal potential for Na+, K+, Ca+, and leak
## C_m, the membrane capacitance
## gna, gk, gca, gl, maximum conductance of sodium, potassium, calcium, and leak
## I applied current to the membrane (set to 10 to have more than one spike)
parameters <- c(Vna = 115, Vk = -11, Vca = 115, Vl = 10.3, gna = 110,
gk = 36, gca = 50, gl= 0.25, C = 1, I = 10)
#
###### 1. Dopaminergic neurons
yini_Dopamin <- c(V = -14, m = 0.25, h = 0.55, n = 0.45, p=0.5) # incl. calcium
###### 2. Muscle neurons
yini_Muscle <- c(V = -14, m = 0.053, h = 0.68, n = 0.318, p=0.0) # excl. calcium
# Model equations
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) }
# ####################################
# Calcium
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
######################################################
## Derivatives
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)))
})
}
# Model app
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
## 1.971 0.026 2.018
print(system.time(
out_modelHH_Muscle <-
ode(yini_Muscle, func=modelHH, times=time, parms=parameters)))
## user system elapsed
## 1.807 0.028 1.869
datatable(summary(out_modelHH_Muscle), caption="Muscle Neurons")
# Specify HH model parameters
parameters <- c(Vna = 115, Vk = -11, Vca = 115, Vl = 10.3, gna = 110,
gk = 36, gca = 5, gl = 0.25, C = 1, I = 10)
# initial state values for Dopaminergic neurons and Muscle neurons
yini_Dopamin <- c(V = -14, m = 0.25, h = 0.55, n = 0.45, p = 0.5) # Dopaminergic neurons
yini_Muscle <- c(V = -14, m = 0.053, h = 0.68, n = 0.318, p = 0.0) # Muscle neurons
# Model equations
modelHH <- function(time, y, parms) {
with(as.list(c(y, parms)), {
# Ion channels functions/Gating functions for each ion channel
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
# Derivatives
dV <- (I - gna*h*(V-Vna)*m^3 - gk*(V-Vk)*n^4 - gca*(V-Vca)*m*h - 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)))
})
}
# Set integration times
time <- seq(from = 0, to = 100, by = 0.01)
# Model integration for Dopaminergic neurons
out_modelHH_Dopamin <- ode(yini_Dopamin, func = modelHH, times = time, parms = parameters)
# Model integration for Muscle neurons
out_modelHH_Muscle <- ode(yini_Muscle, func = modelHH, times = time, parms = parameters)
# probabilities of ion channel activation for dopaminergic neurons
plot(out_modelHH_Dopamin[, "time"], out_modelHH_Dopamin[, "m"], type = "l",
col = "blue", xlab = "Time", ylab = "Probability",
main = "Probabilities of Ion Channel Activation - Dopaminergic Neurons")
lines(out_modelHH_Dopamin[, "time"], out_modelHH_Dopamin[, "n"], col = "red")
lines(out_modelHH_Dopamin[, "time"], out_modelHH_Dopamin[, "p"], col = "green")
legend("topright", legend = c("m (Sodium)", "n (Potassium)", "p (Calcium)"),
col = c("blue", "red", "green"), lty = 1)
# ion channel activation over time for muscle neurons
plot(out_modelHH_Muscle[, "time"], out_modelHH_Muscle[, "m"], type = "l",
col = "blue", xlab = "Time", ylab = "Probability",
main = "Probabilities of Ion Channel Activation - Muscle Neurons")
lines(out_modelHH_Muscle[, "time"], out_modelHH_Muscle[, "n"], col = "red")
lines(out_modelHH_Muscle[, "time"], out_modelHH_Muscle[, "p"], col = "green")
legend("topright", legend = c("m (Sodium)", "n (Potassium)", "p (Calcium)"),
col = c("blue", "red", "green"), lty = 1)
## The spike distribution reflects the frequency and timing of action potentials generated by a neuron. Variations in ion channel probabilities directly influence the excitability and firing properties of neurons. Higher probabilities of ion channel activation lead to increased neuronal excitability and more frequent action potentials, resulting in a higher spike distribution. Conversely, lower probabilities result in decreased excitability and fewer action potentials, leading to a lower spike distribution. Therefore, changes in ion channel probabilities can profoundly impact the behavior of living cells by modulating their electrical activity and signaling patterns, ultimately influencing physiological processes and neural function.
# Define the probabilities for muscle and dopaminergic neurons
prob_Muscle <- c(m = 0.1, h = 0.7, n = 0.4, p = 0)
prob_Dopamin <- c(m = 0.3, h = 0.6, n = 0.5, p = 0.5)
# Model equations for dopaminergic neurons with calcium channels included
modelHH_Dopamin <- function(time, y, parms) {
with(as.list(c(y, parms)), {
alpha_p <- function(v) { 0.01*(15-v)/(exp((15-v)/15)-1) }
beta_p <- function(v) { 0.1*exp(-v/10) }
dV <- (I - gna*h*(V-Vna)*m^3 - gk*(V-Vk)*n^4 - gca*(V-Vca)*m*h - 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
dp <- alpha_p(V)*(1-p) - beta_p(V)*p
return(list(c(dV, dm, dh, dn, dp)))
})
}
# Extract membrane potentials for dopaminergic and muscle neurons
V_Dopamin <- out_modelHH_Dopamin[, "V"]
V_Muscle <- out_modelHH_Muscle[, "V"]
# Calculate spike distributions for dopaminergic and muscle neurons
spike_dist_Dopamin <- rowSums(diff(matrix(V_Dopamin, ncol = 1)) > 0) / length(time)
spike_dist_Muscle <- rowSums(diff(matrix(V_Muscle, ncol = 1)) > 0) / length(time)
# Plot spike distributions for both neuron types on a single plot
plot(time[-1], spike_dist_Muscle, type = 'l', col = 'blue',
xlab = 'Time', ylab = 'Spike Distribution (V/t)',
main = 'Spike Distribution Comparison')
lines(time[-1], spike_dist_Dopamin, type = 'l', col = 'red')
legend("topright", legend = c("Muscle Neurons", "Dopaminergic Neurons"),
col = c("blue", "red"), lty = 1)
## The spike distribution for dopaminergic neurons shows a similar pattern compared to Muscle neurons.Dopaminergic neurons, with their higher probabilities of ion channel activation, exhibit a simillar and frequent spike distribution compared to muscle neurons. This heightened excitability suggests both to generate simillar propensity for dopaminergic neurons to generate action potentials rapidly and sustainably, potentially supporting their role in modulating reward processing and motor control. Muscle neurons, with simillar probabilities of ion channel activation, exhibit a similllar frequent spike distribution, indicative of a more controlled and regulated firing pattern, which aligns with their function in mediating muscle contraction and movement coordination.
# HH model parameters
parameters <- c(Vna = 115, Vk = -11, Vca = 115, Vl = 10.3, gna = 110,
gk = 36, gca = 5, gl = 0.25, C = 1, I = 10)
# Set the initial for Dopaminergic with and without calcium channels
yini_Dopamin_with_Ca <- c(V = -14, m = 0.25, h = 0.55, n = 0.45, p = 0.5) # Dopaminergic neurons with calcium channels
yini_Dopamin_without_Ca <- c(V = -14, m = 0.25, h = 0.55, n = 0.45) # Dopaminergic neurons without calcium channels
# equations for Dopaminergic with calcium channels
modelHH_Dopamin_with_Ca <- 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
dV <- (I - gna*h*(V-Vna)*m^3 - gk*(V-Vk)*n^4 - gca*(V-Vca)*m*h - 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)))
})
}
# equations for Dopaminergic without calcium channels
modelHH_Dopamin_without_Ca <- 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) }
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)))
})
}
# Model with calcium channels
out_modelHH_Dopamin_with_Ca <- ode(yini_Dopamin_with_Ca, func = modelHH_Dopamin_with_Ca, times = time, parms = parameters)
# Model without calcium channels
out_modelHH_Dopamin_without_Ca <- ode(yini_Dopamin_without_Ca, func = modelHH_Dopamin_without_Ca, times = time, parms = parameters)
# spike distributions
spike_distributions_Dopamin_with_Ca <- rowSums(diff(out_modelHH_Dopamin_with_Ca[, "V", drop = FALSE]) > 0) / length(time)
spike_distributions_Dopamin_without_Ca <- rowSums(diff(out_modelHH_Dopamin_without_Ca[, "V", drop = FALSE]) > 0) / length(time)
# Plot spike distributions for Dopaminergic neurons
plot(time[-1], spike_distributions_Dopamin_with_Ca, type = 'l', col = 'blue',
xlab = 'Time', ylab = 'Spike Distribution (V/t)',
main = 'Spike Distribution of Dopaminergic Neurons with and without Ca channels')
lines(time[-1], spike_distributions_Dopamin_without_Ca, type = 'l', col = 'red')
legend("topright", legend = c("With Ca Channels", "Without Ca Channels"),
col = c("blue", "red"), lty = 1)
## The spike distribution for dopaminergic neurons with calcium channels exhibits a higher frequency of spikes compared to those without calcium channels. This suggests that the presence of calcium channels enhances neuronal excitability, potentially leading to more frequent action potentials. The difference in spike distributions highlights the significant impact of calcium channels on the firing behavior of dopaminergic neurons, indicating their crucial role in regulating cellular activity and neurotransmitter release.
# Load required libraries
library(deSolve)
library(ggplot2)
# Specify HH model parameters
parameters <- c(Vna = 50, Vk = -77, Vca = 110, Vl = -54.4, gna = 120,
gk = 36, gca = 5, gl = 0.3, C = 1, I = 10)
# Set the initial state values for Muscle neurons
yini_Muscle_custom <- c(V = -65, m = 0.2, h = 0.7, n = 0.3, p = 0)
# Model equations
modelHH <- function(time, y, parms) {
with(as.list(c(y, parms)), {
# Ion channels functions/Gating functions for each ion channel
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
# Derivatives
dV <- (I - gna*h*(V-Vna)*m^3 - gk*(V-Vk)*n^4 - gca*(V-Vca)*m*h - 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)))
})
}
# Set integration times
time <- seq(from = 0, to = 100, by = 0.01)
# Model integration for custom Muscle neurons
out_modelHH_Muscle_custom <- ode(yini_Muscle_custom, func = modelHH, times = time, parms = parameters)
# Plotting spike profiles of muscle neurons
plot(out_modelHH_Muscle_custom[, "time"], out_modelHH_Muscle_custom[, "V"],
type = "l", col = "blue", xlab = "Time", ylab = "Membrane Potential (mV)",
main = "Spike Profile of Custom Muscle Neurons")
lines(out_modelHH_Muscle[, "time"], out_modelHH_Muscle[, "V"], col = "red")
legend("topright", legend = c("New Initial Conditions", "Original Initial Conditions"),
col = c("blue", "red"), lty = 1)
##By comparing the spike profiles of custom muscle neurons with the original muscle neurons from part a, we observe differences in their membrane potential dynamics. The new muscle neurons exhibit a distinct pattern of spiking, which could suggest altered excitability or ion channel kinetics compared to the original conditions. These differences underscore the sensitivity of neuronal behavior to changes in initial conditions and channel probabilities, reflecting the complex interplay between ion channels and membrane properties in regulating neuronal activity.