HH_model_Cintegrator_testbounds.py

"""A demonstration of a model for a single Hodgkin-Huxley membrane
potential for an oscillatory cortical neuron. Tests bounds events
and dynamic update of bounds between integration runs.

   Robert Clewley, September 2006.
"""
from __future__ import print_function

# textually substitute 'Dopri' for 'Radau' in this file to use Radau

from PyDSTool import *
from time import clock


# ------------------------------------------------------------


def makeHHneuron(name, par_args, ic_args, evs=None, extra_terms='', nobuild=False):
    # extra_terms must not introduce new variables!
    vfn_str = '(I'+extra_terms+'-ionic(v,m,h,n))/C'
    mfn_str = 'ma(v)*(1-m)-mb(v)*m'
    nfn_str = 'na(v)*(1-n)-nb(v)*n'
    hfn_str = 'ha(v)*(1-h)-hb(v)*h'
    aux_str = 'm*m*m*h'

    auxdict = {'ionic': (['vv', 'mm', 'hh', 'nn'],
                              'gna*mm*mm*mm*hh*(vv-vna) + gk*nn*nn*nn*nn*(vv-vk) + gl*(vv-vl)'),
               'ma': (['v'], '0.32*(v+54)/(1-exp(-(v+54)/4))'),
               'mb': (['v'], '0.28*(v+27)/(exp((v+27)/5)-1)'),
               'ha': (['v'], '.128*exp(-(50+v)/18)'),
               'hb': (['v'], '4/(1+exp(-(v+27)/5))'),
               'na': (['v'], '.032*(v+52)/(1-exp(-(v+52)/5))'),
               'nb': (['v'], '.5*exp(-(57+v)/40)')}

    DSargs = {}
    DSargs['varspecs'] = {'v': vfn_str, 'm': mfn_str,
                          'h': hfn_str, 'n': nfn_str,
                          'v_bd0': 'getbound("v",0)',
                          'v_bd1': 'getbound("v",1)',
                          'h_bd0': 'getbound("h",0)'
                             }
    DSargs['pars'] = par_args
    DSargs['auxvars'] = ['v_bd0','v_bd1','h_bd0']
    DSargs['fnspecs'] = auxdict
    DSargs['xdomain'] = {'v': [-100, 50], 'm': [0,1], 'h': [0,1], 'n': [0,1]}
    DSargs['algparams'] = {'init_step': 0.1, 'refine': 0, 'max_step': 0.5}
    DSargs['checklevel'] = 0
    DSargs['ics'] = ic_args
    DSargs['name'] = name
    # nobuild=True so that we can test additional library inclusion
    if nobuild:
        DSargs['nobuild']=True
    if evs is not None:
        DSargs['events'] = evs
    return Generator.Dopri_ODEsystem(DSargs)


# ------------------------------------------------------------


if __name__=='__main__':
    print('-------- Test: Hodgkin-Huxley system')
    par_args = {'gna': 100, 'gk': 80, 'gl': 0.1,
                'vna': 50, 'vk': -100, 'vl': -67,
                'I': 1.75, 'C': 1.0}
    ic_args = {'v':-70.0, 'm': 0, 'h': 1, 'n': 0}

    # test single terminal event first
    print("Testing bounds terminal event and its sampling")

    bd_event = Events.makeZeroCrossEvent('v-getbound("v",1)', 1,
                                         {'name': 'bd_ev',
                                          'eventtol': 1e-4,
                                          'precise': True,
                                          'term': True}, ['v'],
                                         targetlang='c')
    HH = makeHHneuron('HHtest_bdev', par_args, ic_args, [bd_event])
    HH.set(tdata=[0, 25])
    HHtraj = HH.compute('test_term')
    trajdata = HHtraj.sample(dt=0.05, tlo=1.3, thi=20)
    assert HH.getEventTimes()['bd_ev']==[]
    assert trajdata['h_bd0'][0]==0
    assert trajdata['v_bd1'][0]==50
    assert trajdata['v_bd0'][0]==-100

    print('Testing continued integration from t=25, having now set')
    print('voltage domain to be [-100,20]')
    HH.set(xdomain={'v':[-100,20]}, tdata=[0,50])
    HHtraj2 = HH.compute('test_cont', 'c')
    print("Sampled this data up until the event", HH.getEventTimes(), ":")

    plotData = HHtraj2.sample(dt=0.1)
    evt=HH.getEventTimes()['bd_ev']
    yaxislabelstr = 'v'
    plt.ylabel(yaxislabelstr)
    plt.xlabel('t')
    vline=plot(plotData['t'], plotData['v'])
    plot(evt, HHtraj2(evt, 'v'), 'ro')
    show()