Electrical Circuit Model of the Plasma Membrane

For large-scale neural network simulations. Download and compile locally with gfortran.

! Hodgkin-Huxley Model in Fortran 90
! Simulates action potential propagation in squid giant axon
! Author: Cell Physiology Course
! Compile: gfortran -O3 -o hh hodgkin_huxley.f90

program hodgkin_huxley
    implicit none

    ! Physical constants (double precision for accuracy)
    real(8), parameter :: C_m = 1.0d0       ! Membrane capacitance (µF/cm²)
    real(8), parameter :: g_Na_max = 120.0d0 ! Max Na+ conductance (mS/cm²)
    real(8), parameter :: g_K_max = 36.0d0   ! Max K+ conductance (mS/cm²)
    real(8), parameter :: g_L = 0.3d0        ! Leak conductance (mS/cm²)
    real(8), parameter :: E_Na = 50.0d0      ! Na+ reversal potential (mV)
    real(8), parameter :: E_K = -77.0d0      ! K+ reversal potential (mV)
    real(8), parameter :: E_L = -54.4d0      ! Leak reversal potential (mV)

    ! Simulation parameters
    real(8), parameter :: dt = 0.01d0        ! Time step (ms)
    real(8), parameter :: T_max = 50.0d0     ! Total simulation time (ms)
    integer, parameter :: N_steps = nint(T_max / dt)

    ! State variables
    real(8) :: V, m, h, n
    real(8) :: I_Na, I_K, I_L, I_ext
    real(8) :: g_Na, g_K
    real(8) :: t

    ! RK4 intermediate variables
    real(8) :: k1_V, k2_V, k3_V, k4_V
    real(8) :: k1_m, k2_m, k3_m, k4_m
    real(8) :: k1_h, k2_h, k3_h, k4_h
    real(8) :: k1_n, k2_n, k3_n, k4_n
    real(8) :: V_tmp, m_tmp, h_tmp, n_tmp

    integer :: i, output_unit

    ! Initialize at resting state
    V = -65.0d0
    m = alpha_m(V) / (alpha_m(V) + beta_m(V))
    h = alpha_h(V) / (alpha_h(V) + beta_h(V))
    n = alpha_n(V) / (alpha_n(V) + beta_n(V))

    ! Open output file
    output_unit = 10
    open(unit=output_unit, file='hh_output.dat', status='replace')
    write(output_unit, '(A)') '# t(ms)  V(mV)  m  h  n  g_Na  g_K  I_Na  I_K'

    ! Main simulation loop using 4th-order Runge-Kutta
    do i = 0, N_steps
        t = i * dt
        if (t > 10.0d0 .and. t < 40.0d0) then
            I_ext = 10.0d0
        else
            I_ext = 0.0d0
        end if

        g_Na = g_Na_max * m**3 * h
        g_K = g_K_max * n**4
        I_Na = g_Na * (V - E_Na)
        I_K = g_K * (V - E_K)
        I_L = g_L * (V - E_L)

        if (mod(i, 10) == 0) then
            write(output_unit, '(9F12.4)') t, V, m, h, n, g_Na, g_K, I_Na, I_K
        end if

        ! RK4 integration
        k1_V = dVdt(V, m, h, n, I_ext)
        k1_m = dmdt(V, m); k1_h = dhdt(V, h); k1_n = dndt(V, n)

        V_tmp = V + 0.5d0*dt*k1_V
        m_tmp = m + 0.5d0*dt*k1_m
        h_tmp = h + 0.5d0*dt*k1_h
        n_tmp = n + 0.5d0*dt*k1_n
        k2_V = dVdt(V_tmp, m_tmp, h_tmp, n_tmp, I_ext)
        k2_m = dmdt(V_tmp, m_tmp); k2_h = dhdt(V_tmp, h_tmp); k2_n = dndt(V_tmp, n_tmp)

        V_tmp = V + 0.5d0*dt*k2_V
        m_tmp = m + 0.5d0*dt*k2_m
        h_tmp = h + 0.5d0*dt*k2_h
        n_tmp = n + 0.5d0*dt*k2_n
        k3_V = dVdt(V_tmp, m_tmp, h_tmp, n_tmp, I_ext)
        k3_m = dmdt(V_tmp, m_tmp); k3_h = dhdt(V_tmp, h_tmp); k3_n = dndt(V_tmp, n_tmp)

        V_tmp = V + dt*k3_V
        m_tmp = m + dt*k3_m
        h_tmp = h + dt*k3_h
        n_tmp = n + dt*k3_n
        k4_V = dVdt(V_tmp, m_tmp, h_tmp, n_tmp, I_ext)
        k4_m = dmdt(V_tmp, m_tmp); k4_h = dhdt(V_tmp, h_tmp); k4_n = dndt(V_tmp, n_tmp)

        V = V + (dt/6.0d0) * (k1_V + 2.0d0*k2_V + 2.0d0*k3_V + k4_V)
        m = m + (dt/6.0d0) * (k1_m + 2.0d0*k2_m + 2.0d0*k3_m + k4_m)
        h = h + (dt/6.0d0) * (k1_h + 2.0d0*k2_h + 2.0d0*k3_h + k4_h)
        n = n + (dt/6.0d0) * (k1_n + 2.0d0*k2_n + 2.0d0*k3_n + k4_n)

        m = max(0.0d0, min(1.0d0, m))
        h = max(0.0d0, min(1.0d0, h))
        n = max(0.0d0, min(1.0d0, n))
    end do

    close(output_unit)
    print *, 'Simulation complete. Output saved to hh_output.dat'

contains
    real(8) function alpha_m(V)
        real(8), intent(in) :: V
        real(8) :: x
        x = V + 40.0d0
        if (abs(x) < 1.0d-6) then
            alpha_m = 1.0d0
        else
            alpha_m = 0.1d0 * x / (1.0d0 - exp(-x/10.0d0))
        end if
    end function

    real(8) function beta_m(V)
        real(8), intent(in) :: V
        beta_m = 4.0d0 * exp(-(V + 65.0d0) / 18.0d0)
    end function

    real(8) function alpha_h(V)
        real(8), intent(in) :: V
        alpha_h = 0.07d0 * exp(-(V + 65.0d0) / 20.0d0)
    end function

    real(8) function beta_h(V)
        real(8), intent(in) :: V
        beta_h = 1.0d0 / (1.0d0 + exp(-(V + 35.0d0) / 10.0d0))
    end function

    real(8) function alpha_n(V)
        real(8), intent(in) :: V
        real(8) :: x
        x = V + 55.0d0
        if (abs(x) < 1.0d-6) then
            alpha_n = 0.1d0
        else
            alpha_n = 0.01d0 * x / (1.0d0 - exp(-x/10.0d0))
        end if
    end function

    real(8) function beta_n(V)
        real(8), intent(in) :: V
        beta_n = 0.125d0 * exp(-(V + 65.0d0) / 80.0d0)
    end function

    real(8) function dVdt(V, m, h, n, I_ext)
        real(8), intent(in) :: V, m, h, n, I_ext
        real(8) :: I_Na, I_K, I_L
        I_Na = g_Na_max * m**3 * h * (V - E_Na)
        I_K = g_K_max * n**4 * (V - E_K)
        I_L = g_L * (V - E_L)
        dVdt = (I_ext - I_Na - I_K - I_L) / C_m
    end function

    real(8) function dmdt(V, m)
        real(8), intent(in) :: V, m
        dmdt = alpha_m(V) * (1.0d0 - m) - beta_m(V) * m
    end function

    real(8) function dhdt(V, h)
        real(8), intent(in) :: V, h
        dhdt = alpha_h(V) * (1.0d0 - h) - beta_h(V) * h
    end function

    real(8) function dndt(V, n)
        real(8), intent(in) :: V, n
        dndt = alpha_n(V) * (1.0d0 - n) - beta_n(V) * n
    end function
end program hodgkin_huxley