First attempt at creating a hybrid system, events occur, but does not transition between "sub"-dynamical-systems

[Matthew Egbert]

Hello,

Following the example from the PLoS Computational Biology paper appendix (PyDSTool/tests/IF squarespike/model.py), I have tried to implement a minimal hybrid system. I have created two "sub"-dynamical-systems that are connected via a simple zero-crossing event.

The program is simple. A variable G decreases at a constant rate from an initial value of 100. When G-H == 0 (H=50.0 is a parameter) the system should transition to a second "sub-DS". In the first sub-DS, the variable X changes sinusoidally, and in the second it is unchanging.

The program output suggests that the event (G-H==0) occurs, but the dynamics of the system do not appear to change.

Also, I noticed what I think may be a minor, but possibly misleading error in the example. I could be wrong, but in the example, the 'varspec' for the threshold-variable 'excited' is set to '1' in the spiking sub-system. But shouldn't the varspec for this variable always be '0' as it is not really a variable, but more like a flag that does not change except when the system transitions from one sub-DS to the other? The way it is in the example, this variable increases steadily, which does not cause any error in the output, but might be slightly misleading for those trying to learn from the example. Or am I completely mistaken??

Any help as to what I am doing wrong would be much appreciated.

Many thanks,
Matthew

from pylab import *
import PyDSTool as ds

pars = {
    'H' : 50.0,
}

## this allows for a hybrid model such that two different dynamical
## systems are used depending on whether G is < or > H
event_args = {
    'name' : 'thresh_event',
    'eventtol': 1E-4,
    'eventdelay': 1E-1,
    'starttime' : 0,
    'term ': True}

threshold_event = ds.makeZeroCrossEvent(
    'G-H',0, ## have experimented with [-1,0,1]
    event_args,
    varnames = ['G'],
    parnames = ['H'],
)

## variable and ODE definitions
varspecs_above_thresh = {
    'FLAG' : '0', ## have experimented with [0,1] see comment above
                  ## about possible mistake in example
    'X' : 'sin(t)', # above thresh, X should be sinusoidal, and below it should be unchanging
    'G' : '-1.0', # starts high, steadily decreases until passes threshold of H
}
varspecs_below_thresh = varspecs_above_thresh.copy()
varspecs_below_thresh['FLAG'] = '0'
varspecs_below_thresh['X'] = "0"

above_thresh_args = {
    'name' : 'above_thresh_model',
    'tdomain' : [0,Inf],
    'varspecs' : varspecs_above_thresh,
    'pars' : pars,
    'ics' : { 'FLAG' : 1.0 },
    'xdomain' : {'FLAG' : 1.0 },
    'events' : threshold_event,
}

below_thresh_args = {
    'name' : 'below_thresh_model',
    'tdomain' : [0,Inf],
    'varspecs' : varspecs_below_thresh,
    'pars' : pars,
    'ics' : { 'FLAG' : 0.0 },
    'xdomain' : {'FLAG' : 0.0},
}

above_thresh_model = ds.embed(ds.Generator.Vode_ODEsystem(above_thresh_args),
                              name='above_thresh_model')
below_thresh_model = ds.embed(ds.Generator.Vode_ODEsystem(below_thresh_args),
                              name='below_thresh_model')

# model interface
submodel_names = ['below_thresh_model','above_thresh_model']
above_thresh_MIF = ds.intModelInterface(above_thresh_model)
## changes to below_threshold_model when thresh_event event occurs
above_thresh_info = ds.makeModelInfoEntry(above_thresh_MIF,
                                          submodel_names,
                                          [('thresh_event','below_thresh_model')]) 

modelInfoDict = ds.makeModelInfo([above_thresh_info])
mod_args = {'name': 'hybrid_model',
            'modelInfo' : modelInfoDict}
IFmodel = ds.Model.HybridModel(mod_args)

### calculate a trajectory
traj = IFmodel.compute( trajname= 'test', tdata =[0,100],
                        ics={'G': 100, 'X': 0, 'FLAG' : 1}
)
pts = IFmodel.sample('test',dt=0.05)

if len(pts.labels.getLabels()) > 0 :
    print(pts.labels.by_label['Event:thresh_event'])

### Plot results
figure(figsize=(14,5))
for i,var in enumerate(['X','G','FLAG']) :
    plt.subplot2grid((3,1),(i,0))
    plot(pts['t'],pts[var])
    ylabel(var)
    if var == 'G' :
        axhline(pars['H'],color='r')
plt.show()

[Matthew Egbert]

I'm sorry - I should have put this in the 'help' forum. Can it be moved by an administrator?

M

[Matthew Egbert]

Ah! I have found the primary error -- there is a space in the key of this line

'term ': True}
     ^

Which meant that the event was not terminal, and so the event occured, but did not trigger a transition.

M