Polynomializing+Quadratizing Hodgkin Huxley Equations produces Runtime Warning

Hi,

I am trying to convert the Hodgkin Huxley equations to quadratic-bilinear format. Upon calling polynomialize_and_quadratize, I get the following warning message at the start of evaluation:

RuntimeWarning: Substituting new variables failed to produce the expected calculations, so the calculation was done in an alternative mode. If you see this message, please let us know in the Issues: https://github.com/AndreyBychkov/QBee/issues warnings.warn("Substituting new variables failed to produce the expected calculations,"

The equations are defined by the following function:

def HH_qbee():

C_m, g_Na, g_K, g_L, E_Na, E_K, E_L = parameters('C_m, g_Na, g_K, g_L, E_Na, E_K, E_L')
V, m, h, n, I = functions('V, m, h, n, I')

dVdt = (g_Na * m**3 * h * (V - E_Na) - g_K * n ** 4 * (V - E_K) - g_L * (V - E_L) + I) / C_m
# Calculate the gating variable derivatives
dmdt = 0.1 * (V + 40) / (1 - sympy.exp(-0.1 * (V + 40))) * (1 - m) - 4 * sympy.exp(-0.0556 * (V + 65)) * m
dhdt = 0.07 * sympy.exp(-0.05 * (V + 65)) * (1 - h) - 1 / (1 + sympy.exp(-0.1 * (V + 35))) * h
dndt = 0.01 * (V + 55) / (1 - sympy.exp(-0.1 * (V + 55))) * (1 - n) - 0.125 * sympy.exp(-0.0125 * (V + 65)) * n

return [(V, dVdt), (m, dmdt), (h, dhdt), (n, dndt)]

And the qbee function is called as follows:

polynomialize_and_quadratize(HH_qbee(), input_free=True)

I have not gotten the function to successfully terminate (I have let it run for > 1 hour). I have tried passing calc_upper_bound to True as well as providing pruning functions (by_nodes_processed and by_elapsed_time), modfying the call to the qbee function as follows (e.g. for elapsed time):

from functools import partial pruning = partial(pruning_by_elapsed_time, start_t = time(), max_t=60) polynomialize_and_quadratize(HH_qbee(), input_free=True, pruning_functions=pruning)

I do get a message "current best order is 7" very early on in evaluation. This system size would be more than sufficient for my purposes. Is there a way to force early termination?

Thanks for your help!

Hi,

You can try two things: 1) Ctrl+C when "current best order is 7" is shown. It should stop searching and should return a result quickly; 2) Add to pruning functions like this:

upper_bound = partial(pruning_by_vars_number, nvars=7)
timeout = partial(pruning_by_elapsed_time, start_t=time(), max_t=60)  # 60 seconds
res = polynomialize_and_quadratize(HH_qbee(), input_free=True, pruning_functions=[upper_bound, timeout, *default_pruning_rules])

The warning is mostly for us to find edgy cases during polynomialization, so thank you for submitting it :)

I will also update the answer tomorrow!

UPD: It works slowly because of the polynomialization part (I haven't optimized it enough yet), so item "2" is unlikely to work. "1" should work, although discarding other branches seems to take a long time.

There is a bug there (#21), I will fix it soon

Hi,

I have fixed the issue with polynomialization (qbee 0.8.2), so you can get a polynomialized system early. Here is one them:

\begin{array}{ll}
\text{Introduced variables:}     \\                                                                                     
w_0 = 1/(1 - 0.0183156388887342 \exp(-0.1 V)) \\
w_1 = 1/(1 + 0.0301973834223185 \exp(-0.1 V)) \\
w_2 = \exp(-0.05 V) \\
w_3 = \exp(-0.0556 V) \\
w_4 = \exp(-0.0125 V) \\
w_5 = 1/(1 - 0.00408677143846407 \exp(-0.1 V)) \\
w_6 = \exp(-0.1 V) \\
\\
V' = E_K g_K n^4/C_m + E_L g_L/C_m - e_{Na} g_{Na} h m^3/C_m + I/C_m - V g_K n^4/C_m - V g_L/C_m + V g_{Na} h m^3/C_m \\
m' = -0.1 V m w_0 + 0.1 V w_0 - 4.0 m w_0 - 0.107775422421581 m w_3 + 4.0 w_0 \\
h' = -h w_1 - 0.00271419454822054 h w_2 + 0.00271419454822054 w_2 \\
n' = -0.01 V n w_5 + 0.01 V w_5 - 0.055468413760135 n w_4 - 0.55 n w_5 + 0.55 w_5 \\
w_0' = -0.00183156388887342 w_0^2 w_6 (E_K g_K n^4/C_m + E_L g_L/C_m - e_{Na} g_{Na} h m^3/C_m + I/C_m - V g_K n^4/C_m - V g_L/C_m + V g_{Na} h m^3/C_m) \\
w_1' = 0.00301973834223185 w_1^2 w_6 (E_K g_K n^4/C_m + E_L g_L/C_m - e_{Na} g_{Na} h m^3/C_m + I/C_m - V g_K n^4/C_m - V g_L/C_m + V g_{Na} h m^3/C_m) \\
w_2' = -0.05 w_2 (E_K g_K n^4/C_m + E_L g_L/C_m - e_{Na} g_{Na} h m^3/C_m + I/C_m - V g_K n^4/C_m - V g_L/C_m + V g_{Na} h m^3/C_m) \\
w_3' = -0.0556 w_3 (E_K g_K n^4/C_m + E_L g_L/C_m - e_{Na} g_{Na} h m^3/C_m + I/C_m - V g_K n^4/C_m - V g_L/C_m + V g_{Na} h m^3/C_m) \\
w_4' = -0.0125 w_4 (E_K g_K n^4/C_m + E_L g_L/C_m - e_{Na} g_{Na} h m^3/C_m + I/C_m - V g_K n^4/C_m - V g_L/C_m + V g_{Na} h m^3/C_m) \\
w_5' = -0.000408677143846407 w_5^2 w_6 (E_K g_K n^4/C_m + E_L g_L/C_m - e_{Na} g_{Na} h m^3/C_m + I/C_m - V g_K n^4/C_m - V g_L/C_m + V g_{Na} h m^3/C_m) \\
w_6' = -0.1 w_6 (E_K g_K n^4/C_m + E_L g_L/C_m - e_{Na} g_{Na} h m^3/C_m + I/C_m - V g_K n^4/C_m - V g_L/C_m + V g_{Na} h m^3/C_m)
\end{array}

I also need to warn you: there is no input-free quadratization since we have expressions like w_5' = w_5^2·w_6·I + ... which makes it impossible. If it's possible, try to introduce derivatives of I similar to this:

C_m, g_Na, g_K, g_L, E_Na, E_K, E_L = parameters('C_m, g_Na, g_K, g_L, E_Na, E_K, E_L')
V, m, h, n, I = functions('V, m, h, n, I')

dVdt = (g_Na * m ** 3 * h * (V - E_Na) - g_K * n ** 4 * (V - E_K) - g_L * (V - E_L) + I) / C_m
# Calculate the gating variable derivatives
dmdt = (V + 40) / (1 - sympy.exp(-sp.Rational(1, 10) * (V + 40))) * (1 - m) - 4 * sympy.exp(-0.0556 * (V + 65)) * m
dhdt = 0.07 * sympy.exp(-0.05 * (V + 65)) * (1 - h) - 1 / (1 + sympy.exp(-0.1 * (V + 35))) * h
dndt = 0.01 * (V + 55) / (1 - sympy.exp(-0.1 * (V + 55))) * (1 - n) - 0.125 * sympy.exp(-0.0125 * (V + 65)) * n

system = [(V, dVdt), (m, dmdt), (h, dhdt), (n, dndt)]

quad_upper_bound = partial(pruning_by_vars_number, nvars=30)  # We need some upper bound
res = polynomialize_and_quadratize(system, polynomialization_upper_bound=7,
                                       pruning_functions=[*default_pruning_rules, quad_upper_bound],
                                       calc_upper_bound=False,
                                       input_free=False, input_der_orders={I: 2})  # Introduced I' and I''
res.print()
print('-'*50)
res.print(sympy.latex)

It will take a while to compute since such systems with inputs are hard. I would be grateful if you could let us know if you were able to find a quadratization and how long it took.

Hi Andrey,

Thanks so much for the pointers. I have successfully quadratized the system with 32 additional variables. This solution is found almost immediately, and subsequent search does not seem to improve on this. I have one final question - the solution is found upon setting input_free=False. The resulting equations include three references to the input I - \ddot{I} = 0, and then two equations that are listed as (True, True) which presumably refers to I, \dot{I}. As none of the quadratized equations contain any derivatives of I, so I'm a little confused what the augmentation of the state space is accomplishing in this case. For reference, the final quadratized system is given below:

[Eq(y4', w{14}w{4}/8000 - w{14}y4/8000 + 11w{23}/1600 + 11839643642w{4}y4/17075871241891 - 11y4y5/1600), Eq(y5', 10958101147w{13}/4875198079284 + 660041270w{14}w{1}/16150677373043 - 660041270w{14}w{7}/16150677373043 - 3954442325w{1}w{4}/17444525053772 - 10958101147w{7}y5/4875198079284), Eq(y6', -1440481505984w{0}w{10}/19661906345921 + 80422150287w{0}w{15}/43909006273577 - 67774191020w{0}w{6}/34333852436323 - 80422150287w{10}w{20}/43909006273577 + 1440481505984w{12}/19661906345921), True, Eq(I'', 0), Eq(w{0}', 2880963011968w{0}2/19661906345921 - 2880963011968w{10}w{12}/19661906345921 + 160844300574w{12}w{15}/43909006273577 - 160844300574w{12}w{26}/43909006273577 - 135548382040w{12}w{6}/34333852436323), Eq(w{1}', 1320082540w{13}w{14}/16150677373043 - 3954442325w{13}w{4}/8722262526886 + 10958101147*w{1}2/2437599039642 - 10958101147w{1}w{7}/2437599039642 - 1320082540w{21}w{7}/16150677373043), Eq(w{2}', nw{14}/50 + 11ny5/10 - w{14}w{2}/50 - 1313409468016w{18}/11839255704117 - 11w{2}y5/10), Eq(w{3}', 712051904w{19}w{3}/8687630971209 - 712051904w{3}2/8687630971209 + 2998559000w{3}w{5}/22095387391931 + 857976458667w{3}w{8}/28412278198676 - 2998559000w{3}y3/22095387391931 + w{8}y3/20), Eq(w{4}', w{14}w{18}/8000 - 81w{14}w{4}/8000 + w{14}y4/100 + 11w{18}y5/1600 - 891w{23}/1600 + 11839643642*w{4}2/17075871241891 - 656704734008w{4}y4/11839255704117 + 11y4y5/20), Eq(w{5}', hw{19}/20 - w{19} + 2998559000w{5}2/22095387391931 - 2046295261173901452120001w{5}y3/718022453831513814780171480 + 88201785971*y3*2/32496486307080), Eq(w{6}', 139w{10}w{6}/625 - 2639w{15}w{6}/25000 + w{15}y0/10 + 139w{26}w{6}/25000 - 2639w{27}/625 + 167344234781w{6}2/27926482009678 - 2639069740595w{6}y0/24486749217014 + 4y0y6), Eq(w{7}', -3954442325w{13}w{18}/8722262526886 + 11w{13}/20 + w{14}w{1}/100 - w{14}w{7}/100 - 656704734008w{1}w{4}/11839255704117 + 10958101147w{1}w{7}/2437599039642 + 1320082540w{21}w{7}/16150677373043 - 1320082540w{28}w{7}/16150677373043 - 10958101147w{7}2/2437599039642 - 11w{7}y5/20), Eq(w{8}', 1424103808w{19}w{8}/8687630971209 - 1424103808w{3}w{8}/8687630971209 + 88201785971w{3}y2/32496486307080 + 857976458667*w{8}*2/14206139099338 - w{8}y2 - 88201785971w{8}y3/32496486307080), Eq(w{9}', -4w{0}w{10} + w{0}w{15}/10 - 2639069740595w{0}w{6}/24486749217014 + 2880963011968w{10}w{12}/19661906345921 - w{10}w{20}/10 + 160844300574w{12}w{26}/43909006273577 + 4w{12} - 135548382040w{27}w{9}/34333852436323 - 160844300574w{29}w{9}/43909006273577 - 2880963011968w{9}2/19661906345921), Eq(w{10}', 1440481505984w{0}w{10}/19661906345921 + 4w{0} + 80422150287w{10}w{20}/43909006273577 - 67774191020w{10}w{27}/34333852436323 - 80422150287w{10}w{29}/43909006273577 - 1440481505984w{12}w{22}/19661906345921 + w{20}/10 - 2639069740595w{27}/24486749217014 - w{29}/10 - 4w{9}), Eq(w{11}', -12w{10}w{22} + 3w{15}w{22}/10 - 3w{22}w{26}/10 - 7917209221785w{22}w{6}/24486749217014 + 12w{22}y6), Eq(w{12}', 4321444517952w{0}w{12}/19661906345921 + 241266450861w{12}w{20}/43909006273577 - 203322573060w{12}w{27}/34333852436323 - 241266450861w{12}w{29}/43909006273577 - 4321444517952w{12}w{9}/19661906345921), Eq(w{13}', 10958101147w{13}w{1}/1625066026428 + 1980123810w{13}w{21}/16150677373043 - 11863326975w{13}w{23}/17444525053772 - 1980123810w{13}w{28}/16150677373043 - 10958101147w{13}w{7}/1625066026428), Eq(w{14}', -36w{14}w{16} + 10958101147w{14}w{1}/4875198079284 + 660041270w{14}w{21}/16150677373043 - 3954442325w{14}w{23}/17444525053772 + 120w{14}w{24} - 660041270w{14}w{28}/16150677373043 - 10958101147w{14}w{7}/4875198079284 - 3w{14}/10 - 2772w{16}y5 - 6000w{24}y5 - 679487560596473y5/41645219175935), Eq(w{15}', 1440481505984w{0}w{15}/19661906345921 - 1440481505984w{10}w{20}/19661906345921 - 6000w{10}w{25} - 36w{15}w{16} + 80422150287w{15}w{20}/43909006273577 - 67774191020w{15}w{27}/34333852436323 - 3w{15}/10 - 2772w{16}y6 - 80422150287w{20}w{26}/43909006273577 + 120w{25}w{26} - 679487560596473y6/41645219175935), Eq(w{16}', -w{14}w{16}/25 + w{14}w{17}/25 - 11w{16}y5/5 - 2626818936032w{17}w{4}/11839255704117 + 11w{17}y5/5), Eq(w{17}', -3w{14}w{17}/100 + 3w{14}w{2}/100 - 33w{17}y5/20 - 656704734008w{2}w{4}/3946418568039 + 33w{2}y5/20), Eq(w{18}', -161w{14}w{18}/8000 + w{14}w{4}/50 + 11839643642w{18}w{4}/17075871241891 - 1771w{18}y5/1600 + 11w{23}w{2}/1600 + w{23}w{30}/8000 + 11w{23}/10 - 1313409468016*w{4}2/11839255704117), Eq(w{19}', 712051904*w{19}2/8687630971209 - 712051904w{19}w{3}/8687630971209 + 2998559000w{19}w{5}/22095387391931 + 857976458667w{19}w{8}/28412278198676 - 2046295261173901452120001w{3}w{5}/718022453831513814780171480 + 88201785971w{3}y3/32496486307080 + w{5}w{8}/20 - w{8}y3), Eq(w{20}', -2772w{0}w{16} - 679487560596473w{0}/41645219175935 + 2880963011968w{12}w{15}/19661906345921 - 2880963011968w{12}w{26}/19661906345921 - 36w{16}w{20} + 160844300574w{20}2/43909006273577 - 135548382040w{20}w{27}/34333852436323 - 160844300574w{20}w{29}/43909006273577 - 3w{20}/10 + 120w{25}w{29} - 6000w{25}w{9}), Eq(w{21}', 10958101147w{13}w{14}/2437599039642 - 36w{17}w{28} - 2772w{17}w{7} - 6000w{1}w{24} - 679487560596473w{1}/41645219175935 + 1320082540*w{21}2/16150677373043 - 3954442325w{21}w{23}/8722262526886 + 120w{21}w{24} - 1320082540w{21}w{28}/16150677373043 - 10958101147w{21}w{7}/2437599039642 - 3w{21}/10), Eq(w{22}', 8w{10} - w{15}w{22}/5 - 2639069740595w{22}y0/12243374608507 - 8w{22}y6 + w{26}/5), Eq(w{23}', -10958101147w{13}w{18}/4875198079284 + 10958101147w{13}w{4}/4875198079284 - 81w{14}w{23}/8000 + w{18}w{21}/8000 - 891w{1}w{4}/1600 + 11w{1}y4/20 + 660041270w{21}w{23}/16150677373043 + w{21}y4/100 - 3954442325w{23}2/17444525053772 - 660041270w{23}w{28}/16150677373043 + 11839643642w{23}w{4}/17075871241891 - 656704734008w{23}y4/11839255704117 + 11w{4}w{7}/1600), Eq(w{24}', -12w{10}w{25} - 88201785971w{11}w{5}/32496486307080 + 88201785971w{11}y3/32496486307080 + 3w{15}w{25}/10 - w{24}y2 - 3w{25}w{26}/10 - 7917209221785w{25}w{6}/24486749217014 + 12w{25}y6), Eq(w{25}', 8hw{10} + hw{26}/5 - w{15}w{25}/5 - 88201785971w{22}w{5}/32496486307080 + 88201785971w{22}y3/32496486307080 - 2639069740595w{25}y0/12243374608507 - w{25}y2 - 8w{25}y6), Eq(w{26}', -2772w{10}w{16} + 1440481505984w{10}w{20}/19661906345921 - 6000w{10}w{24} - 1440481505984w{10}w{29}/19661906345921 - 679487560596473*w{10}/41645219175935 + w{15}2/10 - w{15}w{26}/10 - 2639069740595w{15}w{6}/24486749217014 - 36w{16}w{26} + 80422150287w{20}w{26}/43909006273577 + 4w{20} + 120w{24}w{26} - 67774191020w{26}w{27}/34333852436323 - 80422150287w{26}w{29}/43909006273577 - 3w{26}/10 - 4w{29}), Eq(w{27}', -2639w{0}w{6}/625 + 4w{0}y0 + 139w{10}w{27}/625 + 1440481505984w{12}w{6}/19661906345921 - 2639w{15}w{27}/25000 + 80422150287w{20}w{27}/43909006273577 + w{20}y0/10 + 139w{26}w{27}/25000 - 67774191020w{27}2/34333852436323 - 80422150287w{27}w{29}/43909006273577 + 167344234781w{27}w{6}/27926482009678 - 1440481505984w{27}w{9}/19661906345921 - 2639069740595w{27}y0/24486749217014), Eq(w{28}', 11w{14}w{1}/20 + w{14}w{21}/100 - 656704734008w{14}w{23}/11839255704117 - w{14}w{28}/100 - 11w{14}w{7}/20 - 36w{16}w{28} - 2772w{16}w{7} + 1320082540w{21}w{28}/16150677373043 + 10958101147w{21}w{7}/2437599039642 - 3954442325w{23}w{28}/8722262526886 + 120w{24}w{28} - 6000w{24}w{7} - 1320082540*w{28}2/16150677373043 - 10958101147w{28}w{7}/2437599039642 - 3w{28}/10 - 679487560596473w{7}/41645219175935), Eq(w{29}', 4w{0}w{15} - 4w{10}w{20} + 2880963011968w{12}w{26}/19661906345921 + w{15}w{20}/10 - 2639069740595w{15}w{27}/24486749217014 - 36w{16}w{29} - 2772w{16}w{9} - w{20}w{26}/10 + 160844300574w{20}w{29}/43909006273577 + 120w{24}w{29} - 6000w{24}w{9} - 135548382040w{27}w{29}/34333852436323 - 160844300574*w{29}*2/43909006273577 - 2880963011968w{29}w{9}/19661906345921 - 3w{29}/10 - 679487560596473w{9}/41645219175935), Eq(w{30}', 11nw{14}/10 + nw{31}/50 - 11w{14}w{2}/10 - w{14}w{30}/50 - 36w{16}w{30} - 2772w{17}*2 - 6000w{24}w{2} + 120w{24}w{30} - 679487560596473w{2}/41645219175935 - 1313409468016w{30}y4/11839255704117 - 3w{30}/10), Eq(w{31}', -5544w{14}w{16} + 10958101147w{14}w{21}/4875198079284 - 12000w{14}w{24} - 10958101147w{14}w{28}/4875198079284 - 1358975121192946w{14}/41645219175935 - 72w{16}w{31} + 660041270w{21}w{31}/16150677373043 - 3954442325w{23}w{31}/17444525053772 + 240w{24}w{31} - 660041270w{28}w{31}/16150677373043 - 3*w{31}/5)]

Hi,

Thanks again for the feedback! Yes, you're right about (True, True) - I fixed that, and I'm also surprised at the absence of I in the quadratization (hope it's just a representation error). I'll investigate it and hope to resolve it this week.

Speaking of which, can I ask you to send me the code so I can reproduce everything accurately?

Hi Andrey,

Code that can be used to reproduce the results is given here: https://github.com/akumar01/QBee/blob/master/hh_share.py. I actually could not get both polynomialization and quadratization to work with qbee. I had been developing my own polynomialization code (that is not particularly sophisticated, so no attempt at optimality) and used that as an input to quadratize.

@akumar01 Sorry for the silence and thank you again for reporting all the issues! We are working on fixing our polynomialization for your model. Meanwhile, I have tried your code, and printed out qbeqs, here is what I got:

(V(_t), 120.0*(V(_t) - 50.0)*h(_t)*m(_t)**3 - 36.0*(V(_t) + 77.0)*n(_t)**4 - 0.3*V(_t) - 16.3161)
(m(_t), (1 - m(_t))*(0.1*V(_t) + 4.0)*y6(_t) - 0.107775422421581*m(_t)*y0(_t))
(h(_t), 0.00271419454822054*(1 - h(_t))*y3(_t) - h(_t)*y2(_t))
(n(_t), (1 - n(_t))*(0.01*V(_t) + 0.55)*y5(_t) - 0.055468413760135*n(_t)*y4(_t))
(y0(_t), -0.0556*((1 - m(_t))*(0.1*V(_t) + 4.0)*y6(_t) - 0.107775422421581*m(_t)*y0(_t))*y0(_t))
(y1(_t), -0.1*((1 - m(_t))*(0.1*V(_t) + 4.0)*y6(_t) - 0.107775422421581*m(_t)*y0(_t))*y1(_t))
(y2(_t), -0.0301973834223185*(0.00271419454822054*(1 - h(_t))*y3(_t) - h(_t)*y2(_t))*y2(_t)**2)
(y3(_t), -0.05*(0.00271419454822054*(1 - h(_t))*y3(_t) - h(_t)*y2(_t))*y3(_t))
(y4(_t), -0.0125*((1 - n(_t))*(0.01*V(_t) + 0.55)*y5(_t) - 0.055468413760135*n(_t)*y4(_t))*y4(_t))
(y5(_t), 0.00408677143846407*((1 - n(_t))*(0.01*V(_t) + 0.55)*y5(_t) - 0.055468413760135*n(_t)*y4(_t))*y5(_t)**2)
(y6(_t), 0.0183156388887342*((1 - m(_t))*(0.1*V(_t) + 4.0)*y6(_t) - 0.107775422421581*m(_t)*y0(_t))*y6(_t)**2)

As far as I can tell, these equations do not contain I at all, so when you further run QBee, no derivatives of I emerge in the RHD because they were not where. Could it happen that I was eliminated by your code?

@akumar01 I have looked carefully at your system, after by-hand polynomialization I get a system with really a lot of nonlinearities. Qbee has been running on it for a couple of hours so far without finding anything. Using some other tricks I think I can produce a quadratization of overall dimension 49 (maybe a bit less). Would this be within the range you think you would find useful? Note that the new variables will involve division, I am not sure how critical is this in your setup.

AndreyBychkov / QBee

Polynomializing+Quadratizing Hodgkin Huxley Equations produces Runtime Warning #20