KeyError 'EP' occurring when plotting nullclines.

[BrittanyH1997]

Hello. I have a 2D system of ODEs and I am trying to plot the nullclines to then vary my model parameters and plot the effects based on the Zero-crossing events tutorial.

Based on the current parameter set, I have three steady states (nullclines are verified using XPP but I digress). My goal is to plot multiple nullclines as a parameter (or parameters) vary to see what my model does. I have currently gotten through the portion where the nullclines are plotted but based on the settings, my figure only gives the directional field and sort of white space where my nullclines should be and the KeyError 'EP' (meaning end point I would assume based on AUTO literature). I'm not sure if because I am setting [2] at the end of the fp_coord() function is working, and depending on the model parameter(s) that is(are) varied, I may get only one intersection by the nullclines as opposed to three (I want to plot this as well). I also want to extend my nullclines slightly into the negative by changing the domains of A1 and A2, but then I get a math domain error (although in some circumstances, my nullclines do slope into the negatives, so I would also be interested in fixing this.)

Any assistance would be fantastic. Thanks so much. You will find my code attached.

import PyDSTool as dst
from PyDSTool import args
import numpy as np
from matplotlib import pyplot as plt
from PyDSTool.Toolbox import phaseplane as pp

pars = {'iota':3.09523809523809535, 'upsilon':27.5, 'vartheta':3.09523809523809535, 'varphi1':10, 'varphi2':10}
icdict = {'A1':0.1, 'A2':0.1}
A1str = '((1+((10*(A1**2.2))/((A1**2.2)+1)))*iota)-(upsilon*A1)+(varphi1*A2)'
A2str = '((1+((10*(A2**2.2))/((A2**2.2)+1)))*vartheta)-(upsilon*A1)+(varphi2*A1)'

event_A1_A2 = dst.makeZeroCrossEvent('A1-A2', 0,
                                   {'name': 'event_A1_A2',
                                   'eventtol':1e-6,
                                   'term': False,
                                   'active': True},
                                   varnames=['A1','A2'],
                                   targetlang='python')

DSargs = args(name='twoDsignalmodel')
DSargs.events = [event_A1_A2]
DSargs.pars = pars
DSargs.tdata = [0,3]
DSargs.algparams = {'max_pts': 1000, 'init_step': 0.01, 'stiff': True}
DSargs.varspecs = {'A1': A1str, 'A2': A2str}
DSargs.xdomain = {'A1': [0.0,2.5], 'A2': [0.0,2.5]}
DSargs.fnspecs = {'Jacobian': (['t','A1','A2'],
                              """[[(10*2.2*(A1**(2.2-1)))/(((A1**2.2)+1)**2)*iota-upsilon, varphi1],
                                  [varphi2, (10*2.2*(A2**(2.2-1)))/(((A2**2.2)+1)**2)*vartheta-upsilon]]""")}
DSargs.ics = icdict
model = dst.Vode_ODEsystem(DSargs)

traj = model.compute('test_traj')
pts = traj.sample()
evs = traj.getEvents('event_A1_A2')

plt.figure(1)
plt.plot(pts['t'], pts['A1'], 'r', linewidth=2)
plt.plot(pts['t'], pts['A2'], 'b', linewidth=2)
plt.show()

plt.figure(2)
pp.plot_PP_vf(model, 'A1', 'A2', scale_exp=-1)
fp_coord= pp.find_fixedpoints(model, n=10, eps=1e-8)[2]
fp = pp.fixedpoint_2D(model, dst.Point(fp_coord), eps=1e-8)
nulls_A1, nulls_A2 = pp.find_nullclines(model, 'A1', 'A2', n=10, eps=1e-8, max_step=0.01, fps=[fp_coord])

# plots fixed point(s)
pp.plot_PP_fps(fp)

# plots the nullclines
plt.plot(nulls_A1[:,0], nulls_A1[:,1], 'g')
plt.plot(nulls_A2[:,0], nulls_A2[:,1], 'y')

# plots the trajectory
plt.plot(pts['A1'], pts['A2'], 'k-o', linewidth=2)

# plots the event points
plt.plot(evs['A1'], evs['A2'], 'rs')

plt.axis('tight')
plt.title('Phase Plane')
plt.xlabel('A1')
plt.ylabel('A2')
plt.show()