O objetivo deste projeto é de resolver numericamente as equações do modelo de Hodgkin- Huxley, simulando a resposta do modelo a uma injeção de corrente do tipo degrau. As equações do modelo são as seguintes e elas representam as funçoes de abertura (alfa) e fechamento (beta) dos parâmetros n, m e h:

\[\alpha_{n}(V) = 0.01*\frac{10-V}{\exp(\frac{10-V}{10}) - 1}\]

\[\beta_{n}(V) = 0.125*exp(\frac{-V}{80})\]

\[\alpha_{m}(V) = 0.1*\frac{25-V}{\exp(\frac{25-V}{10}) - 1}\]

\[\beta_{m}(V) = 4*exp(\frac{-V}{18})\]

\[\alpha_{h}(V) = 0.07*exp(\frac{-V}{20})\]

\[\beta_{h}(V) = \frac{1}{\exp(\frac{30-V}{10}) + 1}\]

alfaN <- function(v){ # can result inf if v = 10
    val <- 0.01*((10 - v)/(exp((10 - v)/10) - 1))
}

betaN <- function(v){
    val <- 0.125 * exp(-v/80)
};

alfaM <- function(v){   # can result inf if v = 10
    val <- 0.1 * ((25 - v)/(exp((25 - v)/10) - 1))
}

betaM <- function(v){ 
    val <- 4 * exp(-v/18)
}

alfaH <- function(v){ 
    val <- 0.07 * exp(-v/20)
}

betaH <- function(v){
    val <- 1/(exp((30 - v)/10) + 1)
}

A voltagem V que aparece nestas equações é o potencial de membrana medido em relação ao potencial de repouso, V = Vm - Vrep, em mV. A densidade de corrente injetada Jinj é medida em \(\mu\)A/cm2, as condutâncias são medidas em mS/cm2 e a capacitância é medida em \(\mu\)F/cm2.

Alfa e Beta são funções lineares mostrando a probabilidade de abertura e fechamento de cada parte do canal iônico.

O significado das variáveis n, m e h é o seguinte. Um canal de Na+ é formado por três portões do tipo m e por um portão do tipo h, cada um dos quais pode estar aberto ou fechado. Supondo que os portões funcionam independentemente, a fração de canais de Na+ abertos a cada instante de tempo é igual a m3h. Igualmente, um canal de K+ é composto por quatro portões do tipo n e a fração de canais de K+ abertos a cada instante é n4.

expInic <- 0; expEnd <- 30; deltat <- 0.025;
expLength <- length(seq(0,expEnd,deltat));
vo <- seq(expInic,expEnd, deltat);
valuesAlfaN <- alfaN(vo);
valuesBetaN = betaN(vo);

par(mfrow = c(1,2));
plot.ts(valuesAlfaN);
plot.ts(valuesBetaN);

plot of chunk unnamed-chunk-2

Criaçao das equações diferenciais, usando o método de euler

steps <- c(length = expLength);
n <- vector(length = expLength);
m <- vector(length = expLength);
h <- vector(length = expLength);
v <- vector(length = steps);
v[1] <- 0;                            # stablish initial conditions
n[1] <- alfaN(0)/(alfaN(0)+betaN(0)); #
m[1] <- alfaM(0)/(alfaM(0)+betaM(0)); #
h[1] <- alfaH(0)/(alfaH(0)+betaH(0)); #

Parâmetros iniciais

gNa <- 120;  ENa <- 115; #sodium conductance and equilibrium
gK <- 36;    EK <- -12;  #potassium conductance and equilibrium
gVaz <- 0.3; EVaz <- 10.163;        #leakage conductance and equilibrium
Cm <- 1;                 #capacitance
Jdens <- 100; Jini <- 5; Jend <- 7; # current density and time interval of injection

Simulaçao do programa principal para a solução do sistema de equações (1)-(4). Voltagem do modelo inicializado em V = Vrep = 0 antes de começar o loop principal da simulação. Implementado solução usando o método de Euler com um passo de integração igual a 0,025 ms.

for (i in 1:(expLength-1)){
    if(i >= Jini & i <=Jend){
        j <- Jdens;
    }else{
        j <- 0;
    }
    #update the main function
    v[i+1] <- v[i] +  (-gNa * m[i]^3 * h[i] * (v[i] - ENa) - gK * n[i]^4 * (v[i] - EK) - gVaz * (v[i] - EVaz) + j) * deltat;
    #update n channel
    if(v[i+1] != 10){ #case when v = 10 (alphaN(v) is indefinite)
        n[i+1] <- n[i] + (alfaN(v[i+1]) * (1 - n[i]) - betaN(v[i+1]) * n[i]) * deltat;
    }else{
        n[i+1] <- n[i] + (0.1 * (1 - n[i]) - betaN(v[i+1]) * n[i]) * deltat;
    }
    #################
    #update m channel
    if(v[i+1] != 25){ #case when v = 25 (alphaM(25) is indefinite)
        m[i+1] <- m[i] + (alfaM(v[i+1]) * (1 - m[i]) - betaM(v[i+1]) * m[i]) * deltat;
    }else{
        m[i+1] <- m[i] + (1 * (1 - m[i]) - betaM(v[i+1]) * m[i]) * deltat;
    }
    #################
    #update h channel
    h[i+1] <- h[i] + (alfaH(v[i+1]) * (1 - h[i]) - betaH(v[i+1]) * h[i]) * deltat;
}

par(mfrow = c(2,1));
plot.ts(n,col = "red",
     main = paste("Pulso ligado em ",Jini," e desligado em ",Jend),
     xlab = "Tempo", ylab = "n, m, h adimensional",
     ylim = c(-0.1,1));
lines(m,col = "green");
lines(h,col = "pink");
legend("topright",legend = c('n','m','h'), lty=c(1,1,1), col=c("red","green","pink"))

plot.ts(v,col = 'blue', main = "Potencial da Membrana (mV)",
     xlab = "Tempo",
     ylab = "mV");

plot of chunk unnamed-chunk-5

Os valores de estado estacionário de n (para o K+) e m (para o Na+) crescem monotonamente à medida que V aumenta, enquanto que o valor de estado estacionário de h (para o Na+) diminui monotonamente à medida que V aumenta. É por isso que n e m são chamadas de variáveis de ativação (do potássio e do sódio, respectivamente) e h é chamada de variável de inativação (do sódio). Observe que o valor da constante de tempo de m é muito menor do que os valores das constantes de tempo de n e h.

Simulaçao da resposta do modelo para densidades de corrente despolarizantes constantes com os seguintes valores (todas em unidades de \(\mu\)A/cm2): 1,0; 5,0; 10,0 e 50,0. Suponha que essas injeções constantes de corrente durem por 100 ms (de t = 5 ms a t = 105 ms) e tempo total de simulação de 105 ms.

expLength <- 3500;
steps <- 105;
n = vector(length = expLength);
m = vector(length = expLength);
h = vector(length = expLength);
v = vector(length = steps);
v[1] = 0;                            # stablish initial conditions
n[1] = alfaN(0)/(alfaN(0)+betaN(0)); #
m[1] = alfaM(0)/(alfaM(0)+betaM(0)); #
h[1] = alfaH(0)/(alfaH(0)+betaH(0)); #

Jdens <- c(1, 5, 10, 170);
Jini <- 5; Jend <- expLength;

par(mfrow = c(2,2))
for (k in 1:length(Jdens)){
    for (i in 1:(expLength-1)){
        if(i >= Jini & i <=Jend){
            j = Jdens[k];
        }else{
            j = 0;
        }
        v[i+1] = v[i] +  (-gNa * m[i]^3 * h[i] * (v[i] - ENa) - gK * n[i]^4 * (v[i] - EK) - gVaz * (v[i] - EVaz) + j) * deltat;
        if(v[i+1] != 10){
            n[i+1] = n[i] + (alfaN(v[i+1]) * (1 - n[i]) - betaN(v[i+1]) * n[i]) * deltat;
        }else{
            n[i+1] = n[i] + (0.1 * (1 - n[i]) - betaN(v[i+1]) * n[i]) * deltat;
        }
        
        if(v[i+1] != 25){
            m[i+1] = m[i] + (alfaM(v[i+1]) * (1 - m[i]) - betaM(v[i+1]) * m[i]) * deltat;
        }else{
            m[i+1] = m[i] + (1 * (1 - m[i]) - betaM(v[i+1]) * m[i]) * deltat;
        }
        h[i+1] = h[i] + (alfaH(v[i+1]) * (1 - h[i]) - betaH(v[i+1]) * h[i]) * deltat;
    }
    plot.ts(v, ylab = "mV", ylim = c(-50,150),col="blue");
    legend("topright",legend = 'v', lty=1, col="blue")
}

plot of chunk unnamed-chunk-6

-Primeiro caso: Não ocorre disparo, apenas uma perturbação muito pequena no início do gráfico, melhor vizualisada com um zoom.

-Segundo caso: Há apenas 1 spike e uma pequena perturbação que desparece rápido.

-Terceiro caso: Ocorre um trem de disparos.

-Quarto caso: Ocorre apenas um spike inicial seguido de oscilações de baixa amplitude.