diff --git a/PySONIC/neurons/cortical.py b/PySONIC/neurons/cortical.py
index 575692f..7c9c2ce 100644
--- a/PySONIC/neurons/cortical.py
+++ b/PySONIC/neurons/cortical.py
@@ -1,393 +1,394 @@
 # -*- coding: utf-8 -*-
 # @Author: Theo Lemaire
 # @Email: theo.lemaire@epfl.ch
 # @Date:   2017-07-31 15:19:51
 # @Last Modified by:   Theo Lemaire
-# @Last Modified time: 2020-01-14 15:24:51
+# @Last Modified time: 2020-01-20 15:22:30
 
 import numpy as np
 from ..core import PointNeuron
 
 
 class Cortical(PointNeuron):
     ''' Generic cortical neuron
 
         Reference:
         *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac,
         Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for
         different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.*
     '''
 
     # ------------------------------ Biophysical parameters ------------------------------
 
     # Resting parameters
     Cm0 = 1e-2  # Membrane capacitance (F/m2)
 
     # Reversal potentials (mV)
     ENa = 50.0   # Sodium
     EK = -90.0   # Potassium
     ECa = 120.0  # Calcium
 
     # ------------------------------ Gating states kinetics ------------------------------
 
     @classmethod
     def alpham(cls, Vm):
         return 0.32 * cls.vtrap(13 - (Vm - cls.VT), 4) * 1e3  # s-1
 
     @classmethod
     def betam(cls, Vm):
         return 0.28 * cls.vtrap((Vm - cls.VT) - 40, 5) * 1e3  # s-1
 
     @classmethod
     def alphah(cls, Vm):
         return 0.128 * np.exp(-((Vm - cls.VT) - 17) / 18) * 1e3  # s-1
 
     @classmethod
     def betah(cls, Vm):
         return 4 / (1 + np.exp(-((Vm - cls.VT) - 40) / 5)) * 1e3  # s-1
 
     @classmethod
     def alphan(cls, Vm):
         return 0.032 * cls.vtrap(15 - (Vm - cls.VT), 5) * 1e3  # s-1
 
     @classmethod
     def betan(cls, Vm):
         return 0.5 * np.exp(-((Vm - cls.VT) - 10) / 40) * 1e3  # s-1
 
     @staticmethod
     def pinf(Vm):
         return 1.0 / (1 + np.exp(-(Vm + 35) / 10))
 
     @classmethod
     def taup(cls, Vm):
         return cls.TauMax / (3.3 * np.exp((Vm + 35) / 20) + np.exp(-(Vm + 35) / 20))  # s
 
     # ------------------------------ States derivatives ------------------------------
 
     @classmethod
     def derStates(cls):
         return {
             'm': lambda Vm, x: cls.alpham(Vm) * (1 - x['m']) - cls.betam(Vm) * x['m'],
             'h': lambda Vm, x: cls.alphah(Vm) * (1 - x['h']) - cls.betah(Vm) * x['h'],
             'n': lambda Vm, x: cls.alphan(Vm) * (1 - x['n']) - cls.betan(Vm) * x['n'],
             'p': lambda Vm, x: (cls.pinf(Vm) - x['p']) / cls.taup(Vm)
         }
 
     # ------------------------------ Steady states ------------------------------
 
     @classmethod
     def steadyStates(cls):
         return {
             'm': lambda Vm: cls.alpham(Vm) / (cls.alpham(Vm) + cls.betam(Vm)),
             'h': lambda Vm: cls.alphah(Vm) / (cls.alphah(Vm) + cls.betah(Vm)),
             'n': lambda Vm: cls.alphan(Vm) / (cls.alphan(Vm) + cls.betan(Vm)),
             'p': lambda Vm: cls.pinf(Vm)
         }
 
     # ------------------------------ Membrane currents ------------------------------
 
     @classmethod
     def iNa(cls, m, h, Vm):
         ''' Sodium current '''
         return cls.gNabar * m**3 * h * (Vm - cls.ENa)  # mA/m2
 
     @classmethod
     def iKd(cls, n, Vm):
         ''' delayed-rectifier Potassium current '''
         return cls.gKdbar * n**4 * (Vm - cls.EK)  # mA/m2
 
     @classmethod
     def iM(cls, p, Vm):
         ''' slow non-inactivating Potassium current '''
         return cls.gMbar * p * (Vm - cls.EK)  # mA/m2
 
     @classmethod
     def iLeak(cls, Vm):
         ''' non-specific leakage current '''
         return cls.gLeak * (Vm - cls.ELeak)  # mA/m2
 
     @classmethod
     def currents(cls):
         return {
             'iNa': lambda Vm, x: cls.iNa(x['m'], x['h'], Vm),
             'iKd': lambda Vm, x: cls.iKd(x['n'], Vm),
             'iM': lambda Vm, x: cls.iM(x['p'], Vm),
             'iLeak': lambda Vm, _: cls.iLeak(Vm)
         }
 
 
 class CorticalRS(Cortical):
     ''' Cortical regular spiking neuron
 
         Reference:
         *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac,
         Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for
         different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.*
     '''
 
     # Neuron name
     name = 'RS'
 
     # ------------------------------ Biophysical parameters ------------------------------
 
     # Resting parameters
     Vm0 = -71.9  # Membrane potential (mV)
 
     # Reversal potentials (mV)
     ELeak = -70.3  # Non-specific leakage
 
     # Maximal channel conductances (S/m2)
     gNabar = 560.0  # Sodium
     gKdbar = 60.0   # Delayed-rectifier Potassium
     gMbar = 0.75    # Slow non-inactivating Potassium
     gLeak = 0.205   # Non-specific leakage
 
     # Additional parameters
     VT = -56.2       # Spike threshold adjustment parameter (mV)
     TauMax = 0.608   # Max. adaptation decay of slow non-inactivating Potassium current (s)
-    area = 11.88e-9  # Cell membrane area (m2)
+    area = 11.84e-9  # Cell membrane area (m2)
 
     # ------------------------------ States names & descriptions ------------------------------
     states = {
         'm': 'iNa activation gate',
         'h': 'iNa inactivation gate',
         'n': 'iKd gate',
         'p': 'iM gate'
     }
 
 
 class CorticalFS(Cortical):
     ''' Cortical fast-spiking neuron
 
         Reference:
         *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac,
         Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for
         different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.*
     '''
 
     # Neuron name
     name = 'FS'
 
     # ------------------------------ Biophysical parameters ------------------------------
 
     # Resting parameters
     Vm0 = -71.4  # Membrane potential (mV)
 
     # Reversal potentials (mV)
     ELeak = -70.4  # Non-specific leakage
 
     # Maximal channel conductances (S/m2)
     gNabar = 580.0  # Sodium
     gKdbar = 39.0   # Delayed-rectifier Potassium
     gMbar = 0.787   # Slow non-inactivating Potassium
     gLeak = 0.38    # Non-specific leakage
 
     # Additional parameters
     VT = -57.9       # Spike threshold adjustment parameter (mV)
     TauMax = 0.502   # Max. adaptation decay of slow non-inactivating Potassium current (s)
     area = 10.17e-9  # Cell membrane area (m2)
 
     # ------------------------------ States names & descriptions ------------------------------
     states = {
         'm': 'iNa activation gate',
         'h': 'iNa inactivation gate',
         'n': 'iKd gate',
         'p': 'iM gate'
     }
 
 
 class CorticalLTS(Cortical):
     ''' Cortical low-threshold spiking neuron
 
         References:
         *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac,
         Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for
         different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.*
 
         *Huguenard, J.R., and McCormick, D.A. (1992). Simulation of the currents involved in
         rhythmic oscillations in thalamic relay neurons. J. Neurophysiol. 68, 1373–1383.*
 
     '''
 
     # Neuron name
     name = 'LTS'
 
     # ------------------------------ Biophysical parameters ------------------------------
 
     # Resting parameters
     Vm0 = -54.0  # Membrane potential (mV)
 
     # Reversal potentials (mV)
     ELeak = -50.0  # Non-specific leakage
 
     # Maximal channel conductances (S/m2)
     gNabar = 500.0  # Sodium
     gKdbar = 40.0   # Delayed-rectifier Potassium
     gMbar = 0.28    # Slow non-inactivating Potassium
     gCaTbar = 4.0   # Low-threshold Calcium
     gLeak = 0.19    # Non-specific leakage
 
     # Additional parameters
-    VT = -50.0    # Spike threshold adjustment parameter (mV)
-    TauMax = 4.0  # Max. adaptation decay of slow non-inactivating Potassium current (s)
-    Vx = -7.0     # Voltage-dependence uniform shift factor at 36°C (mV)
-    area = 25e-9  # Cell membrane area (m2)
+    VT = -50.0       # Spike threshold adjustment parameter (mV)
+    TauMax = 4.0     # Max. adaptation decay of slow non-inactivating Potassium current (s)
+    Vx = -7.0        # Voltage-dependence uniform shift factor at 36°C (mV)
+    area = 25.00e-9  # Cell membrane area (m2)
 
     # ------------------------------ States names & descriptions ------------------------------
     states = {
         'm': 'iNa activation gate',
         'h': 'iNa inactivation gate',
         'n': 'iKd gate',
         'p': 'iM gate',
         's': 'iCaT activation gate',
         'u': 'iCaT inactivation gate'
     }
 
     # ------------------------------ Gating states kinetics ------------------------------
 
     @classmethod
     def sinf(cls, Vm):
         return 1.0 / (1.0 + np.exp(-(Vm + cls.Vx + 57.0) / 6.2))
 
     @classmethod
     def taus(cls, Vm):
         x = np.exp(-(Vm + cls.Vx + 132.0) / 16.7) + np.exp((Vm + cls.Vx + 16.8) / 18.2)
         return 1.0 / 3.7 * (0.612 + 1.0 / x) * 1e-3  # s
 
     @classmethod
     def uinf(cls, Vm):
         return 1.0 / (1.0 + np.exp((Vm + cls.Vx + 81.0) / 4.0))
 
     @classmethod
     def tauu(cls, Vm):
         if Vm + cls.Vx < -80.0:
             return 1.0 / 3.7 * np.exp((Vm + cls.Vx + 467.0) / 66.6) * 1e-3  # s
         else:
             return 1.0 / 3.7 * (np.exp(-(Vm + cls.Vx + 22) / 10.5) + 28.0) * 1e-3  # s
 
     # ------------------------------ States derivatives ------------------------------
 
     @classmethod
     def derStates(cls):
         return {**super().derStates(), **{
             's': lambda Vm, x: (cls.sinf(Vm) - x['s']) / cls.taus(Vm),
             'u': lambda Vm, x: (cls.uinf(Vm) - x['u']) / cls.tauu(Vm)
         }}
 
     # ------------------------------ Steady states ------------------------------
 
     @classmethod
     def steadyStates(cls):
         return {**super().steadyStates(), **{
             's': lambda Vm: cls.sinf(Vm),
             'u': lambda Vm: cls.uinf(Vm)
         }}
 
     # ------------------------------ Membrane currents ------------------------------
 
     @classmethod
     def iCaT(cls, s, u, Vm):
         ''' low-threshold (T-type) Calcium current '''
         return cls.gCaTbar * s**2 * u * (Vm - cls.ECa)  # mA/m2
 
     @classmethod
     def currents(cls):
         return {**super().currents(), **{
             'iCaT': lambda Vm, x: cls.iCaT(x['s'], x['u'], Vm)
         }}
 
 
 class CorticalIB(Cortical):
     ''' Cortical intrinsically bursting neuron
 
         References:
         *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac,
         Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for
         different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.*
 
         *Reuveni, I., Friedman, A., Amitai, Y., and Gutnick, M.J. (1993). Stepwise
         repolarization from Ca2+ plateaus in neocortical pyramidal cells: evidence
         for nonhomogeneous distribution of HVA Ca2+ channels in dendrites.
         J. Neurosci. 13, 4609–4621.*
     '''
 
     # Neuron name
     name = 'IB'
 
     # ------------------------------ Biophysical parameters ------------------------------
 
     # Resting parameters
     Vm0 = -71.4  # Membrane potential (mV)
 
     # Reversal potentials (mV)
     ELeak = -70  # Non-specific leakage
 
     # Maximal channel conductances (S/m2)
     gNabar = 500   # Sodium
     gKdbar = 50    # Delayed-rectifier Potassium
     gMbar = 0.3    # Slow non-inactivating Potassium
     gCaLbar = 1.0  # High-threshold Calcium
     gLeak = 0.1    # Non-specific leakage
 
     # Additional parameters
-    VT = -56.2      # Spike threshold adjustment parameter (mV)
-    TauMax = 0.608  # Max. adaptation decay of slow non-inactivating Potassium current (s)
+    VT = -56.2       # Spike threshold adjustment parameter (mV)
+    TauMax = 0.608   # Max. adaptation decay of slow non-inactivating Potassium current (s)
+    area = 28.95e-9  # Cell membrane area (m2)
 
     # ------------------------------ States names & descriptions ------------------------------
     states = {
         'm': 'iNa activation gate',
         'h': 'iNa inactivation gate',
         'n': 'iKd gate',
         'p': 'iM gate',
         'q': 'iCaL activation gate',
         'r': 'iCaL inactivation gate'
     }
 
     # ------------------------------ Gating states kinetics ------------------------------
 
     @classmethod
     def alphaq(cls, Vm):
         return 0.055 * cls.vtrap(-(Vm + 27), 3.8) * 1e3  # s-1
 
     @staticmethod
     def betaq(Vm):
         return 0.94 * np.exp(-(Vm + 75) / 17) * 1e3  # s-1
 
     @staticmethod
     def alphar(Vm):
         return 0.000457 * np.exp(-(Vm + 13) / 50) * 1e3  # s-1
 
     @staticmethod
     def betar(Vm):
         return 0.0065 / (np.exp(-(Vm + 15) / 28) + 1) * 1e3  # s-1
 
     # ------------------------------ States derivatives ------------------------------
 
     @classmethod
     def derStates(cls):
         return {**super().derStates(), **{
             'q': lambda Vm, x: cls.alphaq(Vm) * (1 - x['q']) - cls.betaq(Vm) * x['q'],
             'r': lambda Vm, x: cls.alphar(Vm) * (1 - x['r']) - cls.betar(Vm) * x['r']
         }}
 
     # ------------------------------ Steady states ------------------------------
 
     @classmethod
     def steadyStates(cls):
         return {**super().steadyStates(), **{
             'q': lambda Vm: cls.alphaq(Vm) / (cls.alphaq(Vm) + cls.betaq(Vm)),
             'r': lambda Vm: cls.alphar(Vm) / (cls.alphar(Vm) + cls.betar(Vm))
         }}
 
     # ------------------------------ Membrane currents ------------------------------
 
     @classmethod
     def iCaL(cls, q, r, Vm):
         ''' high-threshold (L-type) Calcium current '''
         return cls.gCaLbar * q**2 * r * (Vm - cls.ECa)  # mA/m2
 
     @classmethod
     def currents(cls):
         return {**super().currents(), **{
             'iCaL': lambda Vm, x: cls.iCaL(x['q'], x['r'], Vm)
         }}
diff --git a/PySONIC/neurons/stn.py b/PySONIC/neurons/stn.py
index 6be3e2b..a219132 100644
--- a/PySONIC/neurons/stn.py
+++ b/PySONIC/neurons/stn.py
@@ -1,453 +1,455 @@
 # -*- coding: utf-8 -*-
 # @Author: Theo Lemaire
 # @Email: theo.lemaire@epfl.ch
 # @Date:   2018-11-29 16:56:45
 # @Last Modified by:   Theo Lemaire
-# @Last Modified time: 2019-12-09 16:05:13
+# @Last Modified time: 2020-01-20 15:35:44
 
 import numpy as np
 from ..core import PointNeuron
 from ..constants import FARADAY, Z_Ca
 from ..utils import findModifiedEq
 
 
 class OtsukaSTN(PointNeuron):
     ''' Sub-thalamic nucleus neuron
 
         References:
         *Otsuka, T., Abe, T., Tsukagawa, T., and Song, W.-J. (2004). Conductance-Based Model
         of the Voltage-Dependent Generation of a Plateau Potential in Subthalamic Neurons.
         Journal of Neurophysiology 92, 255–264.*
 
         *Tarnaud, T., Joseph, W., Martens, L., and Tanghe, E. (2018). Computational Modeling
         of Ultrasonic Subthalamic Nucleus Stimulation. IEEE Trans Biomed Eng.*
     '''
 
     # Neuron name
     name = 'STN'
 
     # ------------------------------ Biophysical parameters ------------------------------
 
     # Resting parameters
     Cm0 = 1e-2   # Membrane capacitance (F/m2)
     Vm0 = -58.0  # Membrane potential (mV)
     Cai0 = 5e-9  # Intracellular Calcium concentration (M)
 
     # Reversal potentials (mV)
     ENa = 60.0     # Sodium
     EK = -90.0     # Potassium
     ELeak = -60.0  # Non-specific leakage
 
     # Maximal channel conductances (S/m2)
     gNabar = 490.0   # Sodium
     gLeak = 3.5      # Non-specific leakage
     gKdbar = 570.0   # Delayed-rectifier Potassium
     gCaTbar = 50.0   # Low-threshold Calcium
     gCaLbar = 150.0  # High-threshold Calcium
     gAbar = 50.0     # A-type Potassium
     gKCabar = 10.0   # Calcium-dependent Potassium
 
     # Physical constants
     T = 306.15  # K (33°C)
 
     # Calcium dynamics
     Cao = 2e-3         # extracellular Calcium concentration (M)
     taur_Cai = 0.5e-3  # decay time constant for intracellular Ca2+ dissolution (s)
 
     # Fast Na current m-gate
     thetax_m = -40   # mV
     kx_m = -8        # mV
     tau0_m = 0.2e-3  # s
     tau1_m = 3e-3    # s
     thetaT_m = -53   # mV
     sigmaT_m = -0.7  # mV
 
     # Fast Na current h-gate
     thetax_h = -45.5  # mV
     kx_h = 6.4        # mV
     tau0_h = 0e-3     # s
     tau1_h = 24.5e-3  # s
     thetaT1_h = -50   # mV
     thetaT2_h = -50   # mV
     sigmaT1_h = -15   # mV
     sigmaT2_h = 16    # mV
 
     # Delayed rectifier K+ current n-gate
     thetax_n = -41   # mV
     kx_n = -14       # mV
     tau0_n = 0e-3    # s
     tau1_n = 11e-3   # s
     thetaT1_n = -40  # mV
     thetaT2_n = -40  # mV
     sigmaT1_n = -40  # mV
     sigmaT2_n = 50   # mV
 
     # T-type Ca2+ current p-gate
     thetax_p = -56    # mV
     kx_p = -6.7       # mV
     tau0_p = 5e-3     # s
     tau1_p = 0.33e-3  # s
     thetaT1_p = -27   # mV
     thetaT2_p = -102  # mV
     sigmaT1_p = -10   # mV
     sigmaT2_p = 15    # mV
 
     # T-type Ca2+ current q-gate
     thetax_q = -85   # mV
     kx_q = 5.8       # mV
     tau0_q = 0e-3    # s
     tau1_q = 400e-3  # s
     thetaT1_q = -50  # mV
     thetaT2_q = -50  # mV
     sigmaT1_q = -15  # mV
     sigmaT2_q = 16   # mV
 
     # L-type Ca2+ current c-gate
     thetax_c = -30.6  # mV
     kx_c = -5         # mV
     tau0_c = 45e-3    # s
     tau1_c = 10e-3    # s
     thetaT1_c = -27   # mV
     thetaT2_c = -50   # mV
     sigmaT1_c = -20   # mV
     sigmaT2_c = 15    # mV
 
     # L-type Ca2+ current d1-gate
     thetax_d1 = -60   # mV
     kx_d1 = 7.5       # mV
     tau0_d1 = 400e-3  # s
     tau1_d1 = 500e-3  # s
     thetaT1_d1 = -40  # mV
     thetaT2_d1 = -20  # mV
     sigmaT1_d1 = -15  # mV
     sigmaT2_d1 = 20   # mV
 
     # L-type Ca2+ current d2-gate
     thetax_d2 = 0.1e-6  # M
     kx_d2 = 0.02e-6     # M
     tau_d2 = 130e-3     # s
 
     # A-type K+ current a-gate
     thetax_a = -45   # mV
     kx_a = -14.7     # mV
     tau0_a = 1e-3    # s
     tau1_a = 1e-3    # s
     thetaT_a = -40   # mV
     sigmaT_a = -0.5  # mV
 
     # A-type K+ current b-gate
     thetax_b = -90   # mV
     kx_b = 7.5       # mV
     tau0_b = 0e-3    # s
     tau1_b = 200e-3  # s
     thetaT1_b = -60  # mV
     thetaT2_b = -40  # mV
     sigmaT1_b = -30  # mV
     sigmaT2_b = 10   # mV
 
     # Ca2+-activated K+ current r-gate
     thetax_r = 0.17e-6  # M
     kx_r = -0.08e-6     # M
     tau_r = 2e-3        # s
 
+    area = 2.86e-9  # Cell membrane area (m2)
+
     # ------------------------------ States names & descriptions ------------------------------
     states = {
         'm': 'iNa activation gate',
         'h': 'iNa inactivation gate',
         'n': 'iKd gate',
         'a': 'iA activation gate',
         'b': 'iA inactivation gate',
         'p': 'iCaT activation gate',
         'q': 'iCaT inactivation gate',
         'c': 'iCaL activation gate',
         'd1': 'iCaL inactivation gate 1',
         'd2': 'iCaL inactivation gate 2',
         'r': 'iCaK gate',
         'Cai': 'submembrane Calcium concentration (M)'
     }
 
     def __new__(cls):
         cls.deff = cls.getEffectiveDepth(cls.Cai0, cls.Vm0)  # m
         cls.current_to_molar_rate_Ca = cls.currentToConcentrationRate(Z_Ca, cls.deff)
         return super(OtsukaSTN, cls).__new__(cls)
 
     @classmethod
     def getPltScheme(cls):
         pltscheme = super().getPltScheme()
         pltscheme['[Ca^{2+}]_i'] = ['Cai']
         return pltscheme
 
     @classmethod
     def getPltVars(cls, wrapleft='df["', wrapright='"]'):
         pltvars = super().getPltVars(wrapleft, wrapright)
         pltvars['Cai'] = {
             'desc': 'submembrane Ca2+ concentration',
             'label': '[Ca^{2+}]_i',
             'unit': 'uM',
             'factor': 1e6
         }
         return pltvars
 
     @classmethod
     def titrationFunc(cls, *args, **kwargs):
         return cls.isSilenced(*args, **kwargs)
 
     @classmethod
     def getEffectiveDepth(cls, Cai, Vm):
         ''' Compute effective depth that matches a given membrane potential
             and intracellular Calcium concentration.
 
             :return: effective depth (m)
         '''
         iCaT = cls.iCaT(cls.pinf(Vm), cls.qinf(Vm), Vm, Cai)  # mA/m2
         iCaL = cls.iCaL(cls.cinf(Vm), cls.d1inf(Vm), cls.d2inf(Cai), Vm, Cai)  # mA/m2
         return -(iCaT + iCaL) / (Z_Ca * FARADAY * Cai / cls.taur_Cai) * 1e-6  # m
 
     # ------------------------------ Gating states kinetics ------------------------------
 
     @staticmethod
     def _xinf(var, theta, k):
         ''' Generic function computing the steady-state opening of a
             particular channel gate at a given voltage or ion concentration.
 
             :param var: membrane potential (mV) or ion concentration (mM)
             :param theta: half-(in)activation voltage or concentration (mV or mM)
             :param k: slope parameter of (in)activation function (mV or mM)
             :return: steady-state opening (-)
         '''
         return 1 / (1 + np.exp((var - theta) / k))
 
     @classmethod
     def ainf(cls, Vm):
         return cls._xinf(Vm, cls.thetax_a, cls.kx_a)
 
     @classmethod
     def binf(cls, Vm):
         return cls._xinf(Vm, cls.thetax_b, cls.kx_b)
 
     @classmethod
     def cinf(cls, Vm):
         return cls._xinf(Vm, cls.thetax_c, cls.kx_c)
 
     @classmethod
     def d1inf(cls, Vm):
         return cls._xinf(Vm, cls.thetax_d1, cls.kx_d1)
 
     @classmethod
     def d2inf(cls, Cai):
         return cls._xinf(Cai, cls.thetax_d2, cls.kx_d2)
 
     @classmethod
     def minf(cls, Vm):
         return cls._xinf(Vm, cls.thetax_m, cls.kx_m)
 
     @classmethod
     def hinf(cls, Vm):
         return cls._xinf(Vm, cls.thetax_h, cls.kx_h)
 
     @classmethod
     def ninf(cls, Vm):
         return cls._xinf(Vm, cls.thetax_n, cls.kx_n)
 
     @classmethod
     def pinf(cls, Vm):
         return cls._xinf(Vm, cls.thetax_p, cls.kx_p)
 
     @classmethod
     def qinf(cls, Vm):
         return cls._xinf(Vm, cls.thetax_q, cls.kx_q)
 
     @classmethod
     def rinf(cls, Cai):
         return cls._xinf(Cai, cls.thetax_r, cls.kx_r)
 
     @staticmethod
     def _taux1(Vm, theta, sigma, tau0, tau1):
         ''' Generic function computing the voltage-dependent, activation/inactivation time constant
             of a particular ion channel at a given voltage (first variant).
 
             :param Vm: membrane potential (mV)
             :param theta: voltage at which (in)activation time constant is half-maximal (mV)
             :param sigma: slope parameter of (in)activation time constant function (mV)
             :param tau0: minimal time constant (s)
             :param tau1: modulated time constant (s)
             :return: (in)activation time constant (s)
         '''
         return tau0 + tau1 / (1 + np.exp(-(Vm - theta) / sigma))
 
     @classmethod
     def taua(cls, Vm):
         return cls._taux1(Vm, cls.thetaT_a, cls.sigmaT_a, cls.tau0_a, cls.tau1_a)
 
     @classmethod
     def taum(cls, Vm):
         return cls._taux1(Vm, cls.thetaT_m, cls.sigmaT_m, cls.tau0_m, cls.tau1_m)
 
     @staticmethod
     def _taux2(Vm, theta1, theta2, sigma1, sigma2, tau0, tau1):
         ''' Generic function computing the voltage-dependent, activation/inactivation time constant
             of a particular ion channel at a given voltage (second variant).
 
             :param Vm: membrane potential (mV)
             :param theta: voltage at which (in)activation time constant is half-maximal (mV)
             :param sigma: slope parameter of (in)activation time constant function (mV)
             :param tau0: minimal time constant (s)
             :param tau1: modulated time constant (s)
             :return: (in)activation time constant (s)
         '''
         return tau0 + tau1 / (np.exp(-(Vm - theta1) / sigma1) + np.exp(-(Vm - theta2) / sigma2))
 
     @classmethod
     def taub(cls, Vm):
         return cls._taux2(Vm, cls.thetaT1_b, cls.thetaT2_b, cls.sigmaT1_b, cls.sigmaT2_b,
                           cls.tau0_b, cls.tau1_b)
 
     @classmethod
     def tauc(cls, Vm):
         return cls._taux2(Vm, cls.thetaT1_c, cls.thetaT2_c, cls.sigmaT1_c, cls.sigmaT2_c,
                           cls.tau0_c, cls.tau1_c)
 
     @classmethod
     def taud1(cls, Vm):
         return cls._taux2(Vm, cls.thetaT1_d1, cls.thetaT2_d1, cls.sigmaT1_d1, cls.sigmaT2_d1,
                           cls.tau0_d1, cls.tau1_d1)
 
     @classmethod
     def tauh(cls, Vm):
         return cls._taux2(Vm, cls.thetaT1_h, cls.thetaT2_h, cls.sigmaT1_h, cls.sigmaT2_h,
                           cls.tau0_h, cls.tau1_h)
 
     @classmethod
     def taun(cls, Vm):
         return cls._taux2(Vm, cls.thetaT1_n, cls.thetaT2_n, cls.sigmaT1_n, cls.sigmaT2_n,
                           cls.tau0_n, cls.tau1_n)
 
     @classmethod
     def taup(cls, Vm):
         return cls._taux2(Vm, cls.thetaT1_p, cls.thetaT2_p, cls.sigmaT1_p, cls.sigmaT2_p,
                           cls.tau0_p, cls.tau1_p)
 
     @classmethod
     def tauq(cls, Vm):
         return cls._taux2(Vm, cls.thetaT1_q, cls.thetaT2_q, cls.sigmaT1_q, cls.sigmaT2_q,
                           cls.tau0_q, cls.tau1_q)
 
     # ------------------------------ States derivatives ------------------------------
 
     @classmethod
     def derCai(cls, p, q, c, d1, d2, Cai, Vm):
         iCa_tot = cls.iCaT(p, q, Vm, Cai) + cls.iCaL(c, d1, d2, Vm, Cai)
         return - cls.current_to_molar_rate_Ca * iCa_tot - Cai / cls.taur_Cai  # M/s
 
     @classmethod
     def derStates(cls):
         return {
             'a': lambda Vm, x: (cls.ainf(Vm) - x['a']) / cls.taua(Vm),
             'b': lambda Vm, x: (cls.binf(Vm) - x['b']) / cls.taub(Vm),
             'c': lambda Vm, x: (cls.cinf(Vm) - x['c']) / cls.tauc(Vm),
             'd1': lambda Vm, x: (cls.d1inf(Vm) - x['d1']) / cls.taud1(Vm),
             'd2': lambda Vm, x: (cls.d2inf(x['Cai']) - x['d2']) / cls.tau_d2,
             'm': lambda Vm, x: (cls.minf(Vm) - x['m']) / cls.taum(Vm),
             'h': lambda Vm, x: (cls.hinf(Vm) - x['h']) / cls.tauh(Vm),
             'n': lambda Vm, x: (cls.ninf(Vm) - x['n']) / cls.taun(Vm),
             'p': lambda Vm, x: (cls.pinf(Vm) - x['p']) / cls.taup(Vm),
             'q': lambda Vm, x: (cls.qinf(Vm) - x['q']) / cls.tauq(Vm),
             'r': lambda Vm, x: (cls.rinf(x['Cai']) - x['r']) / cls.tau_r,
             'Cai': lambda Vm, x: cls.derCai(x['p'], x['q'], x['c'], x['d1'], x['d2'], x['Cai'], Vm)
         }
 
     # ------------------------------ Steady states ------------------------------
 
     @classmethod
     def Caiinf(cls, p, q, c, d1, Vm):
         ''' Steady-state intracellular Calcium concentration '''
         return findModifiedEq(
             cls.Cai0,
             lambda Cai, p, q, c, d1, Vm: cls.derCai(p, q, c, d1, cls.d2inf(Cai), Cai, Vm),
             p, q, c, d1, Vm)
 
     @classmethod
     def steadyStates(cls):
         lambda_dict = {
             'a': lambda Vm: cls.ainf(Vm),
             'b': lambda Vm: cls.binf(Vm),
             'c': lambda Vm: cls.cinf(Vm),
             'd1': lambda Vm: cls.d1inf(Vm),
             'm': lambda Vm: cls.minf(Vm),
             'h': lambda Vm: cls.hinf(Vm),
             'n': lambda Vm: cls.ninf(Vm),
             'p': lambda Vm: cls.pinf(Vm),
             'q': lambda Vm: cls.qinf(Vm),
         }
         lambda_dict['Cai'] = lambda Vm: cls.Caiinf(
             lambda_dict['p'](Vm),
             lambda_dict['q'](Vm),
             lambda_dict['c'](Vm),
             lambda_dict['d1'](Vm),
             Vm)
         lambda_dict['d2'] = lambda Vm: cls.d2inf(lambda_dict['Cai'](Vm))
         lambda_dict['r'] = lambda Vm: cls.rinf(lambda_dict['Cai'](Vm))
         return lambda_dict
 
     # ------------------------------ Membrane currents ------------------------------
 
     @classmethod
     def iNa(cls, m, h, Vm):
         ''' Sodium current '''
         return cls.gNabar * m**3 * h * (Vm - cls.ENa)  # mA/m2
 
     @classmethod
     def iKd(cls, n, Vm):
         ''' delayed-rectifier Potassium current '''
         return cls.gKdbar * n**4 * (Vm - cls.EK)  # mA/m2
 
     @classmethod
     def iA(cls, a, b, Vm):
         ''' A-type Potassium current '''
         return cls.gAbar * a**2 * b * (Vm - cls.EK)  # mA/m2
 
     @classmethod
     def iCaT(cls, p, q, Vm, Cai):
         ''' low-threshold (T-type) Calcium current '''
         return cls.gCaTbar * p**2 * q * (Vm - cls.nernst(Z_Ca, Cai, cls.Cao, cls.T))  # mA/m2
 
     @classmethod
     def iCaL(cls, c, d1, d2, Vm, Cai):
         ''' high-threshold (L-type) Calcium current '''
         return cls.gCaLbar * c**2 * d1 * d2 * (Vm - cls.nernst(Z_Ca, Cai, cls.Cao, cls.T))  # mA/m2
 
     @classmethod
     def iKCa(cls, r, Vm):
         ''' Calcium-activated Potassium current '''
         return cls.gKCabar * r**2 * (Vm - cls.EK)  # mA/m2
 
     @classmethod
     def iLeak(cls, Vm):
         ''' non-specific leakage current '''
         return cls.gLeak * (Vm - cls.ELeak)  # mA/m2
 
     @classmethod
     def currents(cls):
         return {
             'iNa': lambda Vm, x: cls.iNa(x['m'], x['h'], Vm),
             'iKd': lambda Vm, x: cls.iKd(x['n'], Vm),
             'iA': lambda Vm, x: cls.iA(x['a'], x['b'], Vm),
             'iCaT': lambda Vm, x: cls.iCaT(x['p'], x['q'], Vm, x['Cai']),
             'iCaL': lambda Vm, x: cls.iCaL(x['c'], x['d1'], x['d2'], Vm, x['Cai']),
             'iKCa': lambda Vm, x: cls.iKCa(x['r'], Vm),
             'iLeak': lambda Vm, _: cls.iLeak(Vm)
         }
 
     # ------------------------------ Other methods ------------------------------
 
     @staticmethod
     def getLowIntensities():
         ''' Return an array of acoustic intensities (W/m2) used to study the STN neuron in
             Tarnaud, T., Joseph, W., Martens, L., and Tanghe, E. (2018). Computational Modeling
             of Ultrasonic Subthalamic Nucleus Stimulation. IEEE Trans Biomed Eng.
         '''
         return np.hstack((
             np.arange(10, 101, 10),
             np.arange(101, 131, 1),
             np.array([140])
         ))  # W/m2
diff --git a/PySONIC/neurons/thalamic.py b/PySONIC/neurons/thalamic.py
index 632d145..db3345e 100644
--- a/PySONIC/neurons/thalamic.py
+++ b/PySONIC/neurons/thalamic.py
@@ -1,361 +1,361 @@
 # -*- coding: utf-8 -*-
 # @Author: Theo Lemaire
 # @Email: theo.lemaire@epfl.ch
 # @Date:   2017-07-31 15:20:54
 # @Last Modified by:   Theo Lemaire
-# @Last Modified time: 2020-01-14 15:25:51
+# @Last Modified time: 2020-01-20 15:20:03
 
 import numpy as np
 from ..core import PointNeuron
 from ..constants import Z_Ca
 
 
 class Thalamic(PointNeuron):
     ''' Generic thalamic neuron
 
         Reference:
         *Plaksin, M., Kimmel, E., and Shoham, S. (2016). Cell-Type-Selective Effects of
         Intramembrane Cavitation as a Unifying Theoretical Framework for Ultrasonic
         Neuromodulation. eNeuro 3.*
     '''
 
     # ------------------------------ Biophysical parameters ------------------------------
 
     # Resting parameters
     Cm0 = 1e-2  # Membrane capacitance (F/m2)
 
     # Reversal potentials (mV)
     ENa = 50.0   # Sodium
     EK = -90.0   # Potassium
     ECa = 120.0  # Calcium
 
     # ------------------------------ Gating states kinetics ------------------------------
 
     @classmethod
     def alpham(cls, Vm):
         return 0.32 * cls.vtrap(13 - (Vm - cls.VT), 4) * 1e3  # s-1
 
     @classmethod
     def betam(cls, Vm):
         return 0.28 * cls.vtrap((Vm - cls.VT) - 40, 5) * 1e3  # s-1
 
     @classmethod
     def alphah(cls, Vm):
         return 0.128 * np.exp(-((Vm - cls.VT) - 17) / 18) * 1e3  # s-1
 
     @classmethod
     def betah(cls, Vm):
         return 4 / (1 + np.exp(-((Vm - cls.VT) - 40) / 5)) * 1e3  # s-1
 
     @classmethod
     def alphan(cls, Vm):
         return 0.032 * cls.vtrap(15 - (Vm - cls.VT), 5) * 1e3  # s-1
 
     @classmethod
     def betan(cls, Vm):
         return 0.5 * np.exp(-((Vm - cls.VT) - 10) / 40) * 1e3  # s-1
 
     # ------------------------------ States derivatives ------------------------------
 
     @classmethod
     def derStates(cls):
         return {
             'm': lambda Vm, x: cls.alpham(Vm) * (1 - x['m']) - cls.betam(Vm) * x['m'],
             'h': lambda Vm, x: cls.alphah(Vm) * (1 - x['h']) - cls.betah(Vm) * x['h'],
             'n': lambda Vm, x: cls.alphan(Vm) * (1 - x['n']) - cls.betan(Vm) * x['n'],
             's': lambda Vm, x: (cls.sinf(Vm) - x['s']) / cls.taus(Vm),
             'u': lambda Vm, x: (cls.uinf(Vm) - x['u']) / cls.tauu(Vm)
         }
 
     # ------------------------------ Steady states ------------------------------
 
     @classmethod
     def steadyStates(cls):
         return {
             'm': lambda Vm: cls.alpham(Vm) / (cls.alpham(Vm) + cls.betam(Vm)),
             'h': lambda Vm: cls.alphah(Vm) / (cls.alphah(Vm) + cls.betah(Vm)),
             'n': lambda Vm: cls.alphan(Vm) / (cls.alphan(Vm) + cls.betan(Vm)),
             's': lambda Vm: cls.sinf(Vm),
             'u': lambda Vm: cls.uinf(Vm)
         }
 
     # ------------------------------ Membrane currents ------------------------------
 
     @classmethod
     def iNa(cls, m, h, Vm):
         ''' Sodium current '''
         return cls.gNabar * m**3 * h * (Vm - cls.ENa)  # mA/m2
 
     @classmethod
     def iKd(cls, n, Vm):
         ''' delayed-rectifier Potassium current '''
         return cls.gKdbar * n**4 * (Vm - cls.EK)
 
     @classmethod
     def iCaT(cls, s, u, Vm):
         ''' low-threshold (Ts-type) Calcium current '''
         return cls.gCaTbar * s**2 * u * (Vm - cls.ECa)  # mA/m2
 
     @classmethod
     def iLeak(cls, Vm):
         ''' non-specific leakage current '''
         return cls.gLeak * (Vm - cls.ELeak)  # mA/m2
 
     @classmethod
     def currents(cls):
         return {
             'iNa': lambda Vm, states: cls.iNa(states['m'], states['h'], Vm),
             'iKd': lambda Vm, states: cls.iKd(states['n'], Vm),
             'iCaT': lambda Vm, states: cls.iCaT(states['s'], states['u'], Vm),
             'iLeak': lambda Vm, _: cls.iLeak(Vm)
         }
 
 
 class ThalamicRE(Thalamic):
     ''' Thalamic reticular neuron
 
         References:
         *Destexhe, A., Contreras, D., Steriade, M., Sejnowski, T.J., and Huguenard, J.R. (1996).
         In vivo, in vitro, and computational analysis of dendritic calcium currents in thalamic
         reticular neurons. J. Neurosci. 16, 169–185.*
 
         *Huguenard, J.R., and Prince, D.A. (1992). A novel T-type current underlies prolonged
         Ca(2+)-dependent burst firing in GABAergic neurons of rat thalamic reticular nucleus.
         J. Neurosci. 12, 3804–3817.*
 
     '''
 
     # Neuron name
     name = 'RE'
 
     # ------------------------------ Biophysical parameters ------------------------------
 
     # Resting parameters
     Vm0 = -89.5  # Membrane potential (mV)
 
     # Reversal potentials (mV)
     ELeak = -90.0  # Non-specific leakage
 
     # Maximal channel conductances (S/m2)
     gNabar = 2000.0  # Sodium
     gKdbar = 200.0   # Delayed-rectifier Potassium
     gCaTbar = 30.0   # Low-threshold Calcium
     gLeak = 0.5      # Non-specific leakage
 
     # Additional parameters
-    VT = -67.0    # Spike threshold adjustment parameter (mV)
-    area = 14e-9  # Cell membrane area (m2)
+    VT = -67.0       # Spike threshold adjustment parameter (mV)
+    area = 14.00e-9  # Cell membrane area (m2)
 
     # ------------------------------ States names & descriptions ------------------------------
     states = {
         'm': 'iNa activation gate',
         'h': 'iNa inactivation gate',
         'n': 'iKd gate',
         's': 'iCaT activation gate',
         'u': 'iCaT inactivation gate'
     }
 
     # ------------------------------ Gating states kinetics ------------------------------
 
     @staticmethod
     def sinf(Vm):
         return 1.0 / (1.0 + np.exp(-(Vm + 52.0) / 7.4))
 
     @staticmethod
     def taus(Vm):
         return (1 + 0.33 / (np.exp((Vm + 27.0) / 10.0) + np.exp(-(Vm + 102.0) / 15.0))) * 1e-3  # s
 
     @staticmethod
     def uinf(Vm):
         return 1.0 / (1.0 + np.exp((Vm + 80.0) / 5.0))
 
     @staticmethod
     def tauu(Vm):
         return (28.3 + 0.33 / (
             np.exp((Vm + 48.0) / 4.0) + np.exp(-(Vm + 407.0) / 50.0))) * 1e-3  # s
 
 
 class ThalamoCortical(Thalamic):
     ''' Thalamo-cortical neuron
 
         References:
         *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac, Y.,
         Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for different
         classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.*
 
         *Destexhe, A., Bal, T., McCormick, D.A., and Sejnowski, T.J. (1996). Ionic mechanisms
         underlying synchronized oscillations and propagating waves in a model of ferret
         thalamic slices. J. Neurophysiol. 76, 2049–2070.*
 
         Model of Ca2+ buffering and contribution from iCaT derived from:
         *McCormick, D.A., and Huguenard, J.R. (1992). A model of the electrophysiological
         properties of thalamocortical relay neurons. J. Neurophysiol. 68, 1384–1400.*
     '''
 
     # Neuron name
     name = 'TC'
 
     # ------------------------------ Biophysical parameters ------------------------------
 
     # Resting parameters
     # Vm0 = -63.4  # Membrane potential (mV)
     Vm0 = -61.93  # Membrane potential (mV)
 
     # Reversal potentials (mV)
     EH = -40.0     # Mixed cationic current
     ELeak = -70.0  # Non-specific leakage
 
     # Maximal channel conductances (S/m2)
     gNabar = 900.0  # Sodium
     gKdbar = 100.0  # Delayed-rectifier Potassium
     gCaTbar = 20.0  # Low-threshold Calcium
     gKLeak = 0.138  # Leakage Potassium
     gHbar = 0.175   # Mixed cationic current
     gLeak = 0.1     # Non-specific leakage
 
     # Additional parameters
     VT = -52.0       # Spike threshold adjustment parameter (mV)
     Vx = 0.0         # Voltage-dependence uniform shift factor at 36°C (mV)
     taur_Cai = 5e-3  # decay time constant for intracellular Ca2+ dissolution (s)
     Cai_min = 50e-9  # minimal intracellular Calcium concentration (M)
     deff = 100e-9    # effective depth beneath membrane for intracellular [Ca2+] calculation
     nCa = 4          # number of Calcium binding sites on regulating factor
     k1 = 2.5e22      # intracellular Ca2+ regulation factor (M-4 s-1)
     k2 = 0.4         # intracellular Ca2+ regulation factor (s-1)
     k3 = 100.0       # intracellular Ca2+ regulation factor (s-1)
     k4 = 1.0         # intracellular Ca2+ regulation factor (s-1)
-    area = 29e-9     # Cell membrane area (m2)
+    area = 29.00e-9  # Cell membrane area (m2)
 
     # ------------------------------ States names & descriptions ------------------------------
     states = {
         'm': 'iNa activation gate',
         'h': 'iNa inactivation gate',
         'n': 'iKd gate',
         's': 'iCaT activation gate',
         'u': 'iCaT inactivation gate',
         'Cai': 'submembrane Ca2+ concentration (M)',
         'P0': 'proportion of unbound iH regulating factor',
         'O': 'iH gate open state',
         'C': 'iH gate closed state',
     }
 
     def __new__(cls):
         cls.current_to_molar_rate_Ca = cls.currentToConcentrationRate(Z_Ca, cls.deff)
         return super(ThalamoCortical, cls).__new__(cls)
 
     @staticmethod
     def OL(O, C):
         ''' O-gate locked-open probability '''
         return 1 - O - C
 
     @classmethod
     def getPltScheme(cls):
         pltscheme = super().getPltScheme()
         pltscheme['i_{H}\\ kin.'] = ['O', 'OL', 'P0']
         pltscheme['[Ca^{2+}]_i'] = ['Cai']
         return pltscheme
 
     @classmethod
     def getPltVars(cls, wrapleft='df["', wrapright='"]'):
         return {**super().getPltVars(wrapleft, wrapright), **{
             'Cai': {
                 'desc': 'sumbmembrane Ca2+ concentration',
                 'label': '[Ca^{2+}]_i',
                 'unit': 'uM',
                 'factor': 1e6
             },
             'OL': {
                 'desc': 'iH O-gate locked-opening',
                 'label': 'O_L',
                 'bounds': (-0.1, 1.1),
                 'func': 'OL({0}O{1}, {0}C{1})'.format(wrapleft, wrapright)
             },
             'P0': {
                 'desc': 'iH regulating factor activation',
                 'label': 'P_0',
                 'bounds': (-0.1, 1.1)
             }
         }}
 
     # ------------------------------ Gating states kinetics ------------------------------
 
     @classmethod
     def sinf(cls, Vm):
         return 1.0 / (1.0 + np.exp(-(Vm + cls.Vx + 57.0) / 6.2))
 
     @classmethod
     def taus(cls, Vm):
         x = np.exp(-(Vm + cls.Vx + 132.0) / 16.7) + np.exp((Vm + cls.Vx + 16.8) / 18.2)
         return 1.0 / 3.7 * (0.612 + 1.0 / x) * 1e-3  # s
 
     @classmethod
     def uinf(cls, Vm):
         return 1.0 / (1.0 + np.exp((Vm + cls.Vx + 81.0) / 4.0))
 
     @classmethod
     def tauu(cls, Vm):
         if Vm + cls.Vx < -80.0:
             return 1.0 / 3.7 * np.exp((Vm + cls.Vx + 467.0) / 66.6) * 1e-3  # s
         else:
             return 1 / 3.7 * (np.exp(-(Vm + cls.Vx + 22) / 10.5) + 28.0) * 1e-3  # s
 
     @staticmethod
     def oinf(Vm):
         return 1.0 / (1.0 + np.exp((Vm + 75.0) / 5.5))
 
     @staticmethod
     def tauo(Vm):
         return 1 / (np.exp(-14.59 - 0.086 * Vm) + np.exp(-1.87 + 0.0701 * Vm)) * 1e-3  # s
 
     @classmethod
     def alphao(cls, Vm):
         return cls.oinf(Vm) / cls.tauo(Vm)  # s-1
 
     @classmethod
     def betao(cls, Vm):
         return (1 - cls.oinf(Vm)) / cls.tauo(Vm)  # s-1
 
     # ------------------------------ States derivatives ------------------------------
 
     @classmethod
     def derStates(cls):
         return {**super().derStates(), **{
             'Cai': lambda Vm, x: ((cls.Cai_min - x['Cai']) / cls.taur_Cai -
                                   cls.current_to_molar_rate_Ca * cls.iCaT(x['s'], x['u'], Vm)),  # M/s
             'P0': lambda _, x: cls.k2 * (1 - x['P0']) - cls.k1 * x['P0'] * x['Cai']**cls.nCa,
             'O': lambda Vm, x: (cls.alphao(Vm) * x['C'] - cls.betao(Vm) * x['O'] -
                                 cls.k3 * x['O'] * (1 - x['P0']) + cls.k4 * (1 - x['O'] - x['C'])),
             'C': lambda Vm, x: cls.betao(Vm) * x['O'] - cls.alphao(Vm) * x['C'],
         }}
 
     # ------------------------------ Steady states ------------------------------
 
     @classmethod
     def steadyStates(cls):
         lambda_dict = super().steadyStates()
         lambda_dict['Cai'] = lambda Vm: (cls.Cai_min - cls.taur_Cai * cls.current_to_molar_rate_Ca *
                                          cls.iCaT(cls.sinf(Vm), cls.uinf(Vm), Vm))  # M
         lambda_dict['P0'] = lambda Vm: cls.k2 / (cls.k2 + cls.k1 * lambda_dict['Cai'](Vm)**cls.nCa)
         lambda_dict['O'] = lambda Vm: (cls.k4 / (cls.k3 * (1 - lambda_dict['P0'](Vm)) +
                                        cls.k4 * (1 + cls.betao(Vm) / cls.alphao(Vm))))
         lambda_dict['C'] = lambda Vm: cls.betao(Vm) / cls.alphao(Vm) * lambda_dict['O'](Vm)
         return lambda_dict
 
     # ------------------------------ Membrane currents ------------------------------
 
     @classmethod
     def iKLeak(cls, Vm):
         ''' Potassium leakage current '''
         return cls.gKLeak * (Vm - cls.EK)  # mA/m2
 
     @classmethod
     def iH(cls, O, C, Vm):
         ''' outward mixed cationic current '''
         return cls.gHbar * (O + 2 * cls.OL(O, C)) * (Vm - cls.EH)  # mA/m2
 
     @classmethod
     def currents(cls):
         return {**super().currents(), **{
             'iKLeak': lambda Vm, x: cls.iKLeak(Vm),
             'iH': lambda Vm, x: cls.iH(x['O'], x['C'], Vm)
         }}