# -*- coding: utf-8 -*-

"""
This is the channel module for Averaged Neuron (AN) model. 
"""

__author__ = 'Fumiya Tatsuki, Kensuke Yoshida, Tetsuya Yamada, \
              Takahiro Katsumata, Shoi Shi, Hiroki R. Ueda'
__status__ = 'Published'
__version__ = '1.0.0'
__date__ = '15 May 2020'


import os
import sys
"""
LIMIT THE NUMBER OF THREADS!
change local env variables BEFORE importing numpy
"""
os.environ['OMP_NUM_THREADS'] = '1'
os.environ['NUMEXPR_NUM_THREADS'] = '1'
os.environ['MKL_NUM_THREADS'] = '1'

import numpy as np
from typing import Optional

import params


params = params.Constants()


class Base:
    """ Keep basic attributes and helper functions for each channel.

    Parameters
    ----------
    g : float
        channel condactance
    e : float
        equilibrium (reversal) potential for a channel
    
    Attributes
    ----------
    g : float
        channel conductance
    e : float
        equilibrium (reversal) potential for a channel
    """

    def __init__(self, g: float, e: float) -> None:
        self.g = g
        self.e = e
    
    def set_g(self, new_g: float) -> None:
        """ Set a new conductance for a channel.

        Parameters
        ----------
        new_g : float
            new conductance set for a channel
        """
        self.g = new_g

    def get_g(self) -> float:
        ''' Get current channel conductance value.

        Returns
        ----------
        float
            current conductance
        '''
        return self.g

    def set_e(self, new_e: float) -> None:
        """ Set a new equiribrium potential for a channel

        Parameters
        ----------
            new equiribrium potential for a channel
        """
        self.e = new_e

    def get_e(self) -> float:
        ''' Get current equilibrium potential.

        Returns:
        ----------
        float
            current equilibrium potential.
        '''
        return self.e

class Leak(Base):
    """ Leak channel (sodium / potassium).

    Leak channel can be divided into leak sodium channel and leak
    potassium channel. Usually it doesn't have to be divided. If 
    it is separated, you can conduct more detailed analysis.

    Parameters
    ----------
    g : float
    channel conductance
    e : float
        equiribrium potential for the channel

    Attributes
    ----------
    g : float
        the channel conductance
    e : float
        equiribrium potential for the channel
    """
    def __init__(self, g: Optional[int]=None, e: float=params.vL) -> None:
        super().__init__(g, e)

    def i(self, v:float) -> float:
        """ Calculate current that flows through the channel.

        I = g * (v - e).
        v : membrane potential
        e : equiribrium potential (for sodium ion)

        Parameters
        ----------
        v : float
            membrane potential

        Returns
        ----------
        float
            current that flows through the channel
        """
        return self.g * (v - self.e)

    def set_div(self, vnal: float=params.vNaL, vkl: float=params.vK) -> None:
        """ Setting about deviding leak channel into Na leak and K leak.

        Conductances of leak potassium channel and leak sodium channel are 
        defined as:
        gkl = gleak * (vleak - vnal) / (vkl - vnal)
        gnal = gleak * (vleak - vkl) / (vnal - vkl).
        These definition sattisfies gleak=gkl+gnal.

        Parameters
        ----------
        vk : float
            equilibrium potential for leak potassium channel
        vnal : float
            equilibrium potential for leak sodium channel
        """
        self.vnal = vnal
        self.vkl = vkl
        self.gnal = self.g * (self.e - self.vkl) / (self.vnal - self.vkl)
        self.gkl = self.g * (self.e - self.vnal) / (self.vkl - self.vnal)

    def set_gna(self, new_gnal: float) -> None:
        """ Set a new conductance for a leak sodium channel.

        Parameters
        ----------
        new_gna : float
            new conductance set for a leak sodium channel
        """
        self.gnal = new_gnal
    
    def set_gk(self, new_gkl: float) -> None:
        """ Set a new conductance for a leak potassium channel.

        Parameters
        ----------
        new_gk : float
            new conductance set for a leak potassium channel
        """
        self.gkl = new_gkl

    def ikl(self, v: float) -> float:
        """ Calculate current that flows through the channel.

        I = g * (v - e).
        v : membrane potential
        e : equiribrium potential (for sodium ion)

        Parameters
        ----------
        v : float
            membrane potential

        Returns
        ----------
        float
            current that flows through the channel
        """
        return self.gkl * (v - self.vkl)

    def inal(self, v: float) -> float:
        """ Calculate current that flows through the channel.

        I = g * (v - e).
        v : membrane potential
        e : equiribrium potential (for sodium ion)

        Parameters
        ----------
        v : float
            membrane potential

        Returns
        ----------
        float
            current that flows through the channel
        """
        return self.gnal * (v - self.vnal)
    
    def i_div(self, v: float) -> float:
        return self.inal(v) + self.ikl(v)


class NavHH(Base):
    """ Hodgkin-Huxley type volatage gated sodium channel.

    Note
    ----------
    This formulation is from Compute et al., 2003

    Parameters
    ----------
    g : float or None
        channel conductance, default None
    e : float
        equiribrium potential for the channel, 
        default anmodel.params.Constatns.VNa.

    Attributes
    ----------
    g : float
        HH type sodium channel conductance
    e : float
        equiribrium potential for the channel
        (in most cases, sodium equiribrium potential)
    """

    def __init__(self, g: Optional[float]=None, e: float=params.vNa) -> None:
        super().__init__(g, e)

    def am(self, v: float) -> float:
        """ Calculate voltage-dependent transition rate for activation states.

        In the two state model, activation variable m can be discribed as,
            dm/dt = a(1-m) + bm.
        In this method, 'a' for m can be calculated.

        Parameters
        ----------
        v : float
            membrane potential

        Returns
        ----------
        float
            transition rate for activation states
        """
        if v == -33.:
            return 1.
        else:
            return 0.1 * (v+33.0) / (1.0-np.exp(-(v+33.0)/10.0))

    def bm(self, v: float) -> float:
        """ Calculate voltage-dependent transition rate for activation states.

        'b' in the two state model for m can be calculated.

        Parameters
        ----------
        v : float
            membrane potential

        Returns
        ----------
        float
            transition rate for activation states
        """
        return 4.0 * np.exp(-(v+53.7)/12.0)

    def m_inf(self, v: float) -> float:
        """ Calculate activation variable for steady state.

        In the steady state two state model for activation variable, dm/dt = 0.
        Now, m can be calculated from 0 = a(1-m) + bm. 

        Parameters
        ----------
        v : float
            membrane potential

        Results
        ----------
        float
            activation variable for the channel
        """
        return self.am(v) / (self.am(v) + self.bm(v))

    def ah(self, v: float) -> float:
        """ Calculate voltage-dependent transition rate for inactivation states.

        'a' in the two state model for inactivation variable h can be calculated.

        Parameters
        ----------
        v : float
            membrane potential

        Results
        ----------
        float
            transition rate for inactivation states
        """
        return 0.07 * np.exp(-(v+50.0)/10.0)

    def bh(self, v: float) -> float:
        """ Calculate voltage-dependent transition rate for inactivation states.

        'b' in the two state model for h can be calculated.

        Parameters
        ----------
        v : float
            membrane potential
        
        Results
        ----------
        float
            transition rate for inactivation states
        """
        return 1.0 / (1.0 + np.exp(-(v+20.0)/10.0))

    def h_inf(self, v: float) -> float:
        """ Calculate inactivation rate in steady state.

        In the steady state two state model for activation variable, dh/dt = 0.
        Now, m can be calculated from 0 = a(1-h) + bh. 

        Note
        ----------
        This isn't used in the AN model (not steady state, in this case).

        Parameters
        ----------
        v : float
            membrane potential

        Returns
        ----------
        float
            inactivation variable for the channel
        """
        return self.ah(v) / (self.ah(v) + self.bh(v))

    def h_tau(self, v: float) -> float:
        return 1 / 4 * (self.ah(v)+self.bh(v))

    def dhdt(self, v: float, h: float) -> float:
        """ Differential equation for inactiavtion variable (not in steady state).

        In the two state model, inactivation variable can be formulated as
        dh/dt = a(1-h) + bh.
        In this case, right side is multiplied by constant.

        Parameters
        ----------
        v : float
            membrane potential
        h : float
            inactivation variable

        Returns
        ----------
        float
            dh/dt
        """
        return 4.0 * (self.ah(v)*(1-h) - self.bh(v)*h)

    def i(self, v: float, h: float) -> float:
        """ Calculate current that flows through the channel.

        I = g * m^3 * h * (v - e).
        m : activation variable (steady state)
        h : inactiavtion variable
        v : membrane potential
        e : equiribrium potential (for sodium ion)

        Parameters
        ----------
        v : float
            membrane potential
        h : float
            inactivation variable

        Returns
        ----------
        float
            current that flows through the channel
        """
        return self.g * (self.m_inf(v)**3) * h * (v-self.e)


class KvHH(Base):
    """ Hodgkin-Huxley type voltage gated potassium channel (delayed rectifier).

    Note
    ----------
    This formulation is from Compute et al., 2003

    Parameters
    ----------
    g : float or None
        channel conductance, default None
    e : float
        equiribrium potential for the channel, 
        default anmodel.params.Constants.VK.

    Attributes
    ----------
    g : float
        HH type potassium channel conductance
    e : float
        equiribrium potential for the channel
        (in most cases, potassium equiribrium potential)
    """

    def __init__(self, g: Optional[float]=None, e: float=params.vK) -> None:
        super().__init__(g, e)

    def an(self, v: float) -> float:
        """ Calculate voltage-dependent transition rate for activation states.

        'a' in the two state model for activation state n can be calculated.

        Parameters
        ----------
        v : float
            membrane potential

        Results
        ----------
        float
            transition rate for activation states
        """
        if v == -34.:
            return 0.1
        else:
            return 0.01 * (v+34.0) / (1.0-np.exp(-(v+34.0)/10.0))

    def bn(self, v: float) -> float:
        """ Calculate voltage-dependent transition rate for activation states.

        'b' in the two state model for activation state n can be calculated.

        Parameters
        ----------
        v : float
            membrane potential

        Results
        ----------
        float
            transition rate for activation states
        """
        return 0.125 * np.exp(-(v+44.0)/25.0)

    def n_inf(self, v: float) -> float:
        """ Calculate activation variable in steady state.

        In the steady state two state model for activation state, dn/dt = 0.
        Now, n can be calculated from 0 = a(1-n) + bn. 

        Note
        ----------
        This isn't used in the AN model (not steady state, in this case).

        Parameters
        ----------
        v : float
            membrane potential

        Results
        ----------
        float
            activation variable for the channel
        """
        return self.an(v) / (self.an(v)+self.bn(v))

    def n_tau(self, v: float) -> float:
        return 1 / (4 * (self.an(v) + self.bn(v)))

    def dndt(self, v: float, n: float) -> float:
        """ Differential equation for actiavtion variable (not in steady state).

        In the two state model, activation variable can be formulated as
        dn/dt = a(1-n) + bn.
        In this case, right side is multiplied by constant.

        Parameters
        ----------
        v : float
            membrane potential
        n : float
            activation variable

        Returns
        ----------
        float
            dn/dt
        """
        return 4.0 * (self.an(v)*(1-n)-self.bn(v)*n)

    def i(self, v: float, n: float) -> float:
        """ Calculate current that flows through the channel.

        I = g * n^4 * (v - e).
        n : activation variable (steady state)
        v : membrane potential
        e : equiribrium potential (for potassium ion)

        Parameters
        ----------
        v : float
            membrane potential
        h : float
            inactivation variable

        Returns
        ----------
        float
            current that flows through the channel
        """
        return self.g * n**4 * (v-self.e)


class KvA(Base):
    """ Fast A-type potassium channel.

    Note
    ----------
    This formulation is from Compute et al., 2003

    Parameters
    ----------
    g : float or None
        channel conductance, default None
    e : float
        equiribrium potential for the channel, 
        default anmodel.params.Constants.VK.
    tau : float
        time constant for the differential equation dh/dt, 
        default anmodel.params.Constants.tauA

    Attribute
    ----------
    g : float
        the channel conductance
    e : float
        equiribrium potential for the channel
        (in most cases, potassium equiribrium potential)
    tau : float
        time constant for the differential equation dh/dt
    """
    def __init__(self, g: Optional[float]=None, e: float=params.vK, 
                 tau: float=params.tau_a) -> None:
        super().__init__(g, e)
        self.tau = tau

    def m_inf(self, v: float) -> float:
        """ Calculate activation variable in steady state.

        In this case, we use fitted formulation, not two state model.

        Parameters
        ----------
        v : float
            membrane potential

        Results
        ----------
        float
            activation variable for the channel
        """
        return 1.0 / (1.0 + np.exp(-(v+50.0)/20.0))

    def h_inf(self, v: float) -> float:
        """ Calculate inactivation variable in steady state.

        We use fitted formulation, not two state model.

        Note
        ----------
        This isn't used in the AN model (not steady state, in this case).

        Parameters
        ----------
        v : float
            membrane potential

        Results
        ----------
        float
            inactivation variable for the channel
        """
        return 1.0 / (1.0 + np.exp((v+80.0)/6.0))

    def dhdt(self, v: float, h: float) -> float:
        """ Differential equation for inactiavtion variable (not in steady state).

        dh/dt = (h_inf - h) / tau

        Parameters
        ----------
        v : float
            membrane potential
        h : float
            inactivation variable

        Returns
        ----------
        float
            dh/dt
        """
        return (self.h_inf(v)-h) / self.tau

    def i(self, v: float, h: float) -> float:
        """ Calculate current that flows through the channel.

        I = g * m^3 * h * (v - e).
        m : activation variable (steady state)
        h : inactivation variable
        v : membrane potential
        e : equiribrium potential (for potassium ion)

        Parameters
        ----------
        v : float
            membrane potential
        h : float
            inactivation variable

        Returns
        ----------
        float
            current that flows through the channel
        """
        return self.g * (self.m_inf(v)**3) * h * (v-self.e)


class KvSI(Base):
    """ Slowly inactivating potassium channel (a kind of delayed rectifier).

    Slowly inactivating potassium channel doesn't have inactivation variable.

    Note
    ----------
    This formulation is from Compute et al., 2003

    Parameters
    ----------
    g : float
    channel conductance
    e : float
        equiribrium potential for the channel

    Attributes
    ----------
    g : float
        the channel conductance
    e : float
        equiribrium potential for the channel
        (in most cases, potassium equiribrium potential)
    """
    def __init__(self, g: Optional[float]=None, e: float=params.vK) -> None:
        super().__init__(g, e)

    def m_inf(self, v: float) -> float:
        """ Calculate activation variable in steady state.

        In this case, we use fitted formulation, not two state model.

        Parameters
        ----------
        v : float
            membrane potential

        Results
        ----------
        float
            activation variable for the channel
        """
        return 1.0 / (1.0 + np.exp(-(v+34.0)/6.5))

    def m_tau(self, v: float) -> float:
        """ Calculate time constant for the differential equation dm/dt

        Parameters
        ----------
        v : float
            membrane potential

        Results
        ----------
        float
            time constant for the differential equation dm/dt
        """
        return 8.0 / (np.exp(-(v+55.0)/30.0) + np.exp((v+55.0)/30.0))

    def dmdt(self, v: float, m: float) -> float:
        """ Differential equation for actiavtion variable (not in steady state).

        dm/dt = (m_inf - m) / tau

        Parameters
        ----------
        v : float
            membrane potential
        m : float
            activation variable

        Returns
        ----------
        float
            dm/dt
        """
        return (self.m_inf(v)-m) / self.m_tau(v)

    def i(self, v: float, m: float) -> float:
        """ Calculate current that flows through the channel.

        I = g * m  * (v - e).
        m : activation variable (steady state)
        v : membrane potential
        e : equiribrium potential (for potassium ion)

        Parameters
        ----------
        v : float
            membrane potential
        h : float
            inactivation variable

        Returns
        ----------
        float
            current that flows through the channel
        """
        return self.g * m * (v-self.e)


class Cav(Base):
    """ Voltage-gated calcium channel.

    Parameters
    ----------
    g : float
    channel conductance
    e : float
        equiribrium potential for the channel

    Attributes
    ----------
    g : float
        the channel conductance
    e : float
        equiribrium potential for the channel
        (in most cases, calcium equiribrium potential)
    """
    def __init__(self, g: Optional[float]=None, e: float=params.vCa) -> None:
        super().__init__(g, e)

    def m_inf(self, v: float) -> float:
        """ Calculate activation variable in steady state.

        In this case, we use fitted formulation, not two state model.

        Parameters
        ----------
        v : float
            membrane potential

        Results
        ----------
        float
            activation variable for the channel
        """
        return 1.0 / (1.0 + np.exp(-(v+20.0)/9.0))

    def i(self, v: float) -> float:
        """ Calculate current that flows through the channel.

        I = g * m^2  * (v - e).
        m : activation variable (steady state)
        v : membrane potential
        e : equiribrium potential (for potassium ion)

        Parameters
        ----------
        v : float
            membrane potential

        Returns
        ----------
        float
            current that flows through the channel
        """
        return self.g * self.m_inf(v)**2 * (v-self.e)


class NaP(Base):
    """ Persistent sodium channel.

    Note
    ----------
    This formulation is from Compute et al., 2003

    Parameters
    ----------
    g : float
    channel conductance
    e : float
        equiribrium potential for the channel

    Attributes
    ----------
    g : float
        the channel conductance
    e : float
        equiribrium potential for the channel
        (in most cases, sodium equiribrium potential)
    """
    def __init__(self, g: Optional[float]=None, e: float=params.vNa) -> None:
        super().__init__(g, e)

    def m_inf(self, v: float) -> float:
        """ Calculate activation variable in steady state.

        In this case, we use fitted formulation, not two state model.

        Parameters
        ----------
        v : float
            membrane potential

        Results
        ----------
        float
            activation variable for the channel
        """
        return 1.0 / (1.0 + np.exp(-(v+55.7)/7.7))

    def i(self, v: float) -> float:
        """ Calculate current that flows through the channel.

        I = g * m^3  * (v - e).
        m : activation variable (steady state)
        v : membrane potential
        e : equiribrium potential (for potassium ion)

        Parameters
        ----------
        v : float
            membrane potential

        Returns
        ----------
        float
            current that flows through the channel
        """
        return self.g * self.m_inf(v)**3 * (v-self.e)


class KCa(Base):
    """ Calcium-dependent potassium channel.

    Parameters
    ----------
    g : float
    channel conductance
    e : float
        equiribrium potential for the channel
    kd_ca : float
        dissociation constant of calcium-dependent pottasium channels 

    Attributes
    ----------
    g : float
        the channel conductance
    e : float
        equiribrium potential for the channel
        (in most cases, calcium equiribrium potential)
    kd_ca : float
        dissociation constant of calcium-dependent pottasium channels 
    """
    def __init__(self, g: Optional[float]=None, e: float=params.vK, 
                 kd_ca: float=params.kd_ca) -> None:
        super().__init__(g, e)
        self.kd_ca = kd_ca

    def m_inf(self, ca: float) -> float:
        """ Calculate activation variable in steady state.

        In this case, we use fitted formulation, not two state model.

        Parameters
        ----------
        ca : float
            intracellular calcium concentration

        Results
        ----------
        float
            activation variable for the channel
        """
        return 1.0 / (1.0 + (self.kd_ca/ca)**(3.5))

    def i(self, v: float, ca: float) -> float:
        """ Calculate current that flows through the channel.

        I = g * m^3  * (v - e).
        m : activation variable (steady state)
        v : membrane potential
        e : equiribrium potential (for potassium ion)

        Parameters
        ----------
        v : float
            membrane potential
        ca : float
            intracellular calcium concentration

        Returns
        ----------
        float
            current that flows through the channel
        """
        return self.g * self.m_inf(ca) * (v-self.e)


class KIR(Base):
    """ Inwardly rectifying potassium channel.

    Note
    ----------
    This formulation is from Compute et al., 2003

    Parameters
    ----------
    g : float or None
    channel conductance, default None
    e : float
        equiribrium potential for the channel

    Attributes
    ----------
    g : float
        the channel conductance
    e : float
        equiribrium potential for the channel
        (in most cases, calcium equiribrium potential)
    """
    def __init__(self, g: Optional[float]=None, e: float=params.vK) -> None:
        super().__init__(g, e)

    def h_inf(self, v: float) -> float:
        """ Calculate inactivation variable in steady state.

        We use fitted formulation, not two state model.

        Parameters
        ----------
        v : float
            membrane potential

        Results
        ----------
        float
            inactivation variable for the channel
        """
        return 1.0/(1.0 + np.exp((v + 75.0)/4.0))

    def i(self, v: float) -> float:
        """ Calculate current that flows through the channel.

        I = g * h  * (v - e).
        h : inactivation variable (steady state)
        v : membrane potential
        e : equiribrium potential (for potassium ion)

        Parameters
        ----------
        v : float
            membrane potential

        Returns
        ----------
        float
            current that flows through the channel
        """
        return self.g * self.h_inf(v) * (v-self.e)


class AMPAR(Base):
    """ AMPA receptor.

    Note
    ----------
    This formulation is from Compute et al., 2003

    Parameters
    ----------
    g : float or None
    receptor conductance, default None
    e : float
        equiribrium potential for the receptor
    s_a : float
        coefficient of f(V)
    tau_a : float
        time constant for differential equation of gating variable s

    Attributes
    ----------
    g : float
        the receptor conductance
    e : float
        equiribrium potential for the receptor
    s_a : float
        coefficient of f(V)
    tau_a : float
        time constant for differential equation of gating variable s
    """
    def __init__(self, g: Optional[float]=None, e: float=params.vAMPAR, 
                 s_a: float=params.s_a_ampar, s_tau: float=params.s_tau_ampar) -> None:
        super().__init__(g, e)
        self.s_a = s_a
        self.tau_a = s_tau

    def f(self, v: float) -> float:
        """ Function that converts membrane potential into firing rate.

        Parameters
        ----------
        v : float
            membrane potential

        Returns
        ----------
        float
            f(v) : firing rate
        """
        return 1.0 / (1.0 + np.exp(-(v-20.0)/2.0))

    def dsdt(self, v: float, s: float) -> float:
        """ Differential equation for gating variable s.

        ds/dt = af(Vpre) - s/tau

        Parameters
        ----------
        v : float
            membrane potential
        s : float
            gating variable

        Returns
        ----------
        float
            ds/dt
        """
        return self.s_a * self.f(v) - s/self.tau_a

    def i(self, v: float, s: float) -> float:
        """ Calculate current that flows by the receptor.

        I = g * s  * (v - e).
        g : gating variable
        v : membrane potential
        e : equiribrium potential (for potassium ion)

        Parameters
        ----------
        v : float
            membrane potential
        s : float
            gating varriable

        Returns
        ----------
        float
            current that flows by the receptor
        """
        return self.g * s * (v - self.e)


class NMDAR(Base):
    """ NMDA receptor.

    The gating variable of NMDA receptor follows a two first-order 
    kinetic scheme so that EPSC has a slower rise phase and saturates
    at high presynaptic firing rates.

    Note
    ----------
    This formulation is from Compute et al., 2003 and Wang, 1999

    Parameters
    ----------
    g : float or None
    receptor conductance, default None
    e : float
        equiribrium potential for the receptor
    s_a : float
        coefficient of (1 - s)
    s_tau : float
        time constant for differential equation of gating variable s
    x_a : float
        coefficient of f(V)
    x_tau : float
        time constant for differential equation of second-order
        gating variable x

    Attributes
    ----------
    g : float
        the receptor conductance
    e : float
        equiribrium potential for the receptor
    s_a : float
        coefficient of f(V)
    tau_a : float
        time constant for differential equation of gating variable
    """
    def __init__(self, g: Optional[float]=None, e: float=params.vNMDAR, 
                 s_a: float=params.s_a_nmdar, s_tau: float=params.s_tau_nmdar, 
                 x_a: float=params.x_a_nmdar, x_tau: float=params.x_tau_nmdar,
                 ion: bool=False, ex_mg: Optional[float]=None) -> None:
        super().__init__(g, e)
        self.s_a = s_a
        self.s_tau = s_tau
        self.x_a = x_a
        self.x_tau = x_tau
        self.ion = ion
        self.ex_mg = ex_mg

    def f(self, v: float) -> float:
        """ Function that converts membrane potential into firing rate.

        Parameters
        ----------
        v : float
            membrane potential

        Returns
        ----------
        float
            f(v) : firing rate
        """
        return 1.0 / (1.0 + np.exp(-(v-20.0)/2.0))

    def dxdt(self, v: float, x: float) -> float:
        """ Differential equation for second-order gating variable x.

        dx/dt = af(Vpre) - x/tau

        Parameters
        ----------
        v : float
            membrane potential
        x : float
            second-order gating variable

        Returns
        ----------
        float
            dx/dt
        """
        return self.x_a * self.f(v) - x/self.x_tau

    def dsdt(self, v: float, s: float, x: float) -> float:
        """ Differential equation for gating variable s.

        ds/dt = ax(1 - s) - s/tau

        Parameters
        ----------
        v : float
            membrane potential
        s : float
            gating variable
        x : float
            second-order gating variable

        Returns
        ----------
        float
            ds/dt
        """
        return self.s_a * x * (1-s) - s/self.s_tau

    def i(self, v: float, s: float) -> float:
        """ Calculate current that flows by the receptor.

        I = a * g * s  * (v - e).
        a : scaling variable, 1 (not considering ion concentration) or
            1.1 / (1+Mg/8) (considering ion concentration)
        g : gating variable
        v : membrane potential
        e : equiribrium potential (for potassium ion)

        Parameters
        ----------
        v : float
            membrane potential
        s : float
            gating varriable

        Returns
        ----------
        float
            current that flows by the receptor
        """
        if self.ion:
            return 1.1 / (1.0+self.ex_mg/8.0) * self.g * s * (v-self.e)
        else:
            return self.g * s * (v-self.e)


class GABAR(Base):
    """ GABA receptor.

    Note
    ----------
    This formulation is from Compute et al., 2003

    Parameters
    ----------
    g : float or None
    receptor conductance, default None
    e : float
        equiribrium potential for the receptor
    s_a : float
        coefficient of f(V)
    tau_a : float
        time constant for differential equation of gating variable s

    Attributes
    ----------
    g : float
        the receptor conductance
    e : float
        equiribrium potential for the receptor
    s_a : float
        coefficient of f(V)
    tau_a : float
        time constant for differential equation of gating variable s
    """
    def __init__(self, g: Optional[float]=None, e: float=params.vGABAR, 
                 s_a: float=params.s_a_gabar, s_tau: float=params.s_tau_gabar) -> None:
        super().__init__(g, e)
        self.s_a = s_a
        self.s_tau = s_tau

    def f(self, v: float) -> float:
        """ Function that converts membrane potential into firing rate.

        Parameters
        ----------
        v : float
            membrane potential

        Returns
        ----------
        float
            f(v) : firing rate
        """
        return 1.0 / (1.0 + np.exp(-(v-20.0)/2.0))

    def dsdt(self, v: float, s: float) -> float:
        """ Differential equation for gating variable s.

        ds/dt = af(Vpre) - s/tau

        Parameters
        ----------
        v : float
            membrane potential
        s : float
            gating variable

        Returns
        ----------
        float
            ds/dt
        """
        return self.s_a * self.f(v) - s/self.s_tau

    def i(self, v: float, s: float) -> float:
        """ Calculate current that flows by the receptor.

        I = g * s  * (v - e).
        g : gating variable
        v : membrane potential
        e : equiribrium potential (for potassium ion)

        Parameters
        ----------
        v : float
            membrane potential
        s : float
            gating varriable

        Returns
        ----------
        float
            current that flows by the receptor
        """
        return self.g * s * (v-self.e)