Beginner questions and bugreport in PyCont for verbosity >= 10

Help
Vincent
2014-05-28
2015-05-11
  • Vincent
    Vincent
    2014-05-28

    Hi,

    I recently tried to use bifurcation analysis methods for my work (systems biology), unfortunately with no success for now, using several softwares/libraries. I gave a try to PyCont in the last days (using last version from github), without much success either.

    I'm trying to reproduced published results (from http://dx.doi.org/10.1016/j.bpj.2013.02.012), where a model exhibit bistable behaviour. When I try to compute Equilibrium-Point Curve, I get this error :

    Computing curve...
    Traceback (most recent call last):
      File "PyCont_Gerard.py", line 40, in <module>             
        PyCont['EQ1'].forward()
      File "/home/labestiol/pydstool/PyDSTool/PyCont/Continuation.py", line 1243, in forward
        self._compute(v0=self.initdirec)
      File "/home/labestiol/pydstool/PyDSTool/PyCont/Continuation.py", line 1047, in _compute
        raise PyDSTool_ExistError("Problem in _compute: Failed to compute tangent vector.")
    PyDSTool.errors.PyDSTool_ExistError: 'Problem in _compute: Failed to compute tangent vector.'
    

    Trying to debug this, I increased the verbosity, and found this bug :

    Computing curve...
    Traceback (most recent call last):
      File "PyCont_Gerard.py", line 40, in <module>
        PyCont['EQ1'].forward()
      File "/home/labestiol/pydstool/PyDSTool/PyCont/Continuation.py", line 1243, in forward
        self._compute(v0=self.initdirec)
      File "/home/labestiol/pydstool/PyDSTool/PyCont/Continuation.py", line 1038, in _compute
        print("%d: %lf" % (d,log10(cond(r_[c_[J_coords, J_params], [perpvec]]))))
    NameError: global name 'cond' is not defined
    

    Hope this will help you guys, this library seems really nice, and I still have hope I'll be able to use it soon :D
    And just in case my case seems like a classic beginner mistake, here is my code :

    from PyDSTool import *
    from matplotlib import pyplot as plt
    
    pars = {'cdk': 0.2, 'K1S': 0.002, 'K2S': 0.04, 'K1M': 0.04, 'K2M': 0.002, 'k1S': 3.6, 'V2S': 1, 'k1M': 1, 'V2M': 0.2}
    icdict = {'Sp': 0.102045, 'S': 0.897955, 'Mp': 0.000394, 'M': 0.999606}
    
    u1str = 'k1S*cdk*S/(K1S*(1+(M/K1M)) + S) - V2S*Sp/(K2S*(1+(Mp/K2M)) + Sp)'
    u2str = 'V2S*Sp/(K2S*(1+(Mp/K2M)) + Sp) - k1S*cdk*S/(K1S*(1+(M/K1M)) + S)'
    u3str = 'k1M*cdk*M/(K1M*(1+(S/K1S)) + M) - V2M*Mp/(K2M*(1+(Sp/K2S)) + Mp)'
    u4str = 'V2M*Mp/(K2M*(1+(Sp/K2S)) + Mp) - k1M*cdk*M/(K1M*(1+(S/K1S)) + M)'
    
    DSargs = args(name='Gerard')
    DSargs.pars = pars
    DSargs.varspecs = {'Sp': u1str, 'S': u2str, 'Mp': u3str, 'M': u4str}
    DSargs.ics = icdict
    testDS = Generator.Vode_ODEsystem(DSargs)
    
    PyCont = ContClass(testDS)
    
    PCargs = args(name='EQ1', type='EP-C')
    PCargs.freepars = ['cdk']
    PCargs.StepSize = 0.01
    PCargs.MaxNumPoints = 500
    PCargs.MaxStepSize = 0.1
    PCargs.MinStepSize = 1e-5
    PCargs.LocBifPoints = 'All'
    PCargs.verbosity = 2
    PCargs.SaveEigen = True
    PyCont.newCurve(PCargs)
    
    print('Computing curve...')
    start = clock()
    PyCont['EQ1'].forward()
    print('done in %.3f seconds!' % (clock()-start))
    

    The choice of the [min|max] stepsize still looks a little like black magic to me, so no clue if this choice is good or not. My initial values are supposed to be a at equilibrium.
    Anyway, I'll keep looking !

    Thanks

     
    Last edit: Vincent 2014-05-28
  • Rob Clewley
    Rob Clewley
    2014-05-28

    Hi Vincent,

    First, there seem to be some missing parameters in your model description. It looks like all of the numerators in the second fraction of each ODE RHS are missing definitions.

    When not used as a temporary float to store the condition number, the symbol 'cond' should be a function from the numpy.linalg module. I suspect that you have found a bug in the code that is only run when verbosity >= 10, and that it should be prefixed as 'linalg.cond' on that line and one a few below. Try making that change and see what information you get back when PyCont fails to converge. Drew can hopefully also confirm that this is the right change to make here.

     
  • Vincent
    Vincent
    2014-05-28

    Hi Rob,

    First of all thanks for the answer.
    I tried the fix you suggested, and it's not working :

    File "/home/labestiol/pydstool/PyDSTool/PyCont/Continuation.py", line 1038, in _compute
        print("%d: %lf" % (d,log10(linalg.cond(r_[c_[J_coords, J_params], [perpvec]]))))
    AttributeError: 'module' object has no attribute 'cond'
    

    I also triple-check, and I don't see any parameter definition missing. I believe this kind of error would trigger something like "ValueError: Undeclared or illegal token K1S in spec string M", which it doesn't.

    My problem (with bifurcation analysis in general, same thing with xpp, auto and others) seems to always come back to a convergence failure. I initially thought the problem was with initial conditions not being at equilibrium, but even with that it doesn't work. Do you have any idea how to tackle this, or do you have some links which explain how to solve this problem ?

     
  • Rob Clewley
    Rob Clewley
    2014-05-28

    Drew will be able to help with improving convergence. I had a copy-paste problem that dropped some '*' symbols, so I corrected that.

    Also, it turns out that scipy.linalg doesn't contain cond, only numpy.linalg (despite online claims to the contrary - maybe there's a version issue here). So, restore 'linalg.cond' to 'cond' and add 'from numpy.linalg import cond' at the end of the import list. That at least lets me run your code, to find:

    Computing curve...
    2: 20.576081
    3: 20.820803
    4: 20.006843
    5: 19.686496
    /sw/lib/python2.7/site-packages/numpy/linalg/linalg.py:1425: RuntimeWarning: divide by zero encountered in double_scalars
      return s[0]/s[-1]
    6: inf
    
     
  • Drew LaMar
    Drew LaMar
    2014-05-28

    Taking a look at your system, you have a lot of symmetry going on. In particular, the right-hand side of the differential equation for S is just the negative of the right-hand side of the differential equation for Sp. This gives you a singular jacobian matrix, and thus continuation cannot occur. I suggest you look at the reduced 2D (S,M) system by substituting Sp = C1-S and Mp = C2-M (from the fact that dSp/dt + dS/dt = 0 and dMp+dM = 0), where the constants C1 and C2 are found from the initial conditions.

    Hope that helps!

    • Drew
     
  • Vincent
    Vincent
    2014-05-28

    Thanks A LOT guys !

    Rob correction did the trick. I now have tons of debugging information :)
    And Drew suggestion did solve the convergence issue. I'm able to reproduce the results :)

    There's no word to describe how much I love you guys.

     
  • Wiber
    Wiber
    2015-05-10

    Hi Vincent,

    would you mind posting the updated and working file?! This would be wonderful and help me a lot - thanks!

     
  • Vincent
    Vincent
    2015-05-10

    Hi Wiber,

    I did as Drew told me, and reduced the model. Here is the reduced definition :

    pars = {'cdk': 0, 'K1S': 0.002, 'K2S': 0.04, 'K1M': 0.04, 'K2M': 0.002, 'k1S': 3.6, 'V2S': 1, 'k1M': 1, 'V2M': 0.2, 'TS': 1, 'TM': 1}
    icdict = {'Sp': 0, 'Mp': 0}
    
    u1str = 'k1S*cdk*(TS-Sp)/(K1S*(1+((TM-Mp)/K1M)) + (TS-Sp)) - V2S*Sp/(K2S*(1+(Mp/K2M)) + Sp)'
    u3str = 'k1M*cdk*(TM-Mp)/(K1M*(1+((TS-Sp)/K1S)) + (TM-Mp)) - V2M*Mp/(K2M*(1+(Sp/K2S)) + Mp)'
    
    DSargs = args(name='Gerard')
    DSargs.pars = pars
    DSargs.varspecs = {'Sp': u1str, 'Mp': u3str}
    DSargs.ics = icdict
    testDS = Generator.Vode_ODEsystem(DSargs)
    

    The rest of the code is the same

     
  • Wiber
    Wiber
    2015-05-11

    Hi Vincent,

    Thanks a lot for your immediate reply and the updated code! That allows to reproduce the figures 1b and 1c easily, great!

     
    Last edit: Wiber 2015-05-11