.. currentmodule:: brian2

.. Morris_Lecar_1981:

Example: Morris_Lecar_1981
==========================


        .. only:: html

            .. |launchbinder| image:: file:///usr/share/doc/python-brian-doc/docs/badge.svg
            .. _launchbinder: https://mybinder.org/v2/gh/brian-team/brian2-binder/master?filepath=examples/frompapers/Morris_Lecar_1981.ipynb

            .. note::
               You can launch an interactive, editable version of this
               example without installing any local files
               using the Binder service (although note that at some times this
               may be slow or fail to open): |launchbinder|_

        

Morris-Lecar model

Reproduces Fig. 9 of:

Catherine Morris and Harold Lecar.
“Voltage Oscillations in the Barnacle Giant Muscle Fiber.”
Biophysical Journal 35, no. 1 (1981): 193–213.

::

    from brian2 import *
    set_device('cpp_standalone')
    defaultclock.dt = 0.01*ms
    
    g_L = 2*mS
    g_Ca = 4*mS
    g_K = 8*mS
    V_L = -50*mV
    V_Ca = 100*mV
    V_K = -70*mV
    lambda_n__max = 1.0/(15*ms)
    V_1 = 10*mV
    V_2 = 15*mV  # Note that Figure caption says -15 which seems to be a typo
    V_3 = -1*mV
    V_4 = 14.5*mV
    C = 20*uF
    
    
    # V,N-reduced system (Eq. 9 in article), note that the variables M and N (and lambda_N, etc.)
    # have been renamed to m and n to better match the Hodgkin-Huxley convention, and because N has
    # a reserved meaning in Brian (number of neurons)
    eqs = '''
    dV/dt = (-g_L*(V - V_L) - g_Ca*m_inf*(V - V_Ca) - g_K*n*(V - V_K) + I)/C : volt
    dn/dt = lambda_n*(n_inf - n) : 1
    m_inf = 0.5*(1 + tanh((V - V_1)/V_2)) : 1
    n_inf = 0.5*(1 + tanh((V - V_3)/V_4)) : 1
    lambda_n = lambda_n__max*cosh((V - V_3)/(2*V_4)) : Hz
    I : amp
    '''
    
    neuron = NeuronGroup(17, eqs, method='exponential_euler')
    neuron.I = (np.arange(17)*25+100)*uA
    neuron.V = V_L
    neuron.n = 'n_inf'
    mon = StateMonitor(neuron, ['V', 'n'], record=True)
    
    run_time = 220*ms
    run(run_time)
    fig, (ax1, ax2) = plt.subplots(1, 2, gridspec_kw={'right': 0.95, 'bottom': 0.15},
                                   figsize=(6.4, 3.2))
    fig.subplots_adjust(wspace=0.4)
    for line_no, idx in enumerate([0, 4, 12, 15]):
        color = 'C%d' % line_no
        ax1.plot(mon.t/ms, mon.V[idx]/mV, color=color)
        ax1.text(225, mon.V[idx][-1]/mV, '%.0f' % (neuron.I[idx]/uA), color=color)
    ax1.set(xlim=(0, 220), ylim=(-50, 50), xlabel='time (ms)')
    ax1.set_ylabel('V (mV)', rotation=0)
    ax1.spines['right'].set_visible(False)
    ax1.spines['top'].set_visible(False)
    
    # dV/dt nullclines
    V = linspace(-50, 50, 100)*mV
    for line_no, (idx, color) in enumerate([(0, 'C0'), (4, 'C1'), (8, 'C4'), (12, 'C2'), (16, 'C5')]):
        n_null = (g_L*(V - V_L) + g_Ca*0.5*(1 + tanh((V - V_1)/V_2))*(V - V_Ca) - neuron.I[idx])/(-g_K*(V - V_K))
        ax2.plot(V/mV, n_null, color=color)
        ax2.text(V[20+5*line_no]/mV, n_null[20+5*line_no]+0.01, '%.0f' % (neuron.I[idx]/uA), color=color)
    # dn/dt nullcline
    n_null = 0.5*(1 + tanh((V - V_3)/V_4))
    ax2.plot(V/mV, n_null, color='k')
    ax2.set(xlim=(-50, 50), ylim=(0, 1), xlabel='V (mV)')
    ax2.set_ylabel('n', rotation=0)
    ax2.spines['right'].set_visible(False)
    ax2.spines['top'].set_visible(False)
    plt.show()
    

