<jupyter_start><jupyter_text>Up states: Stability of SynapticScaling learning rule Soldado-Magraner, Laje & Buonomano (2021)  Author: <NAME>  Neural subsystem Load script:  (RUN FIRST `upstates-Neural subsystem stability.ipynb`)<jupyter_code>%%capture
load('upstates-Neural subsystem stability.sage.py')<jupyter_output><empty_output><jupyter_text>Synaptic subsystem: SynapticScaling  SynapticScaling equations<jupyter_code>var('dWEEdt,dWEIdt,dWIEdt,dWIIdt')
var('alpha_EE,alpha_EI,alpha_IE,alpha_II')
f_WEE = dWEEdt == alpha_EE*W_EE*(E_set-E)
f_WEI = dWEIdt == -alpha_EI*W_EI*(E_set-E)
f_WIE = dWIEdt == alpha_IE*W_IE*(I_set-I)
f_WII = dWIIdt == -alpha_II*W_II*(I_set-I)
show(f_WEE)
show(f_WEI)
show(f_WIE)
show(f_WII)
synapticFixedPoint = solve([f_WEE.subs(dWEEdt==0),f_WEI.subs(dWEIdt==0),f_WIE.subs(dWIEdt==0),f_WII.subs(dWIIdt==0)],E,I)
synapticFixedPoint_up = synapticFixedPoint[0]<jupyter_output><empty_output><jupyter_text>Only one solution compatible with suprathreshold activity:  Fixed point (Up state):  $E=E_{set}$, $I=I_{set}$ (the weights can take any value subject to $E_{ss}=E_{set}$ and $I_{ss}=I_{set}$)<jupyter_code>show(synapticFixedPoint_up)<jupyter_output><empty_output><jupyter_text>Quasi-steady-state approximation1. Full system neural+synaptic2. Nondimensionalization3. Normalization4. Approximation (fast neural subsystem)  Full system<jupyter_code>show(f_E)
show(f_I)
show(f_WEE)
show(f_WEI)
show(f_WIE)
show(f_WII)<jupyter_output><empty_output><jupyter_text>Nondimensionalization  Definition of nondimensional variables  Define new (nondimensional) variables $e$, $i$, etc. and scaling parameters $E_0$, $I_0$, etc:  $e = E/E_0$  $i = I/I_0$  $\tau = t/\tau_0$  $w_{EE} = W_{EE}/W_{EE0}$  ... etc<jupyter_code>var('e,i')
var('E_0,I_0,tau_0')
var('w_EE,w_EI,w_IE,w_II')
var('W_EE0,W_EI0,W_IE0,W_II0')
var('dedtau,didtau')
var('dwEEdtau,dwEIdtau,dwIEdtau,dwIIdtau');<jupyter_output><empty_output><jupyter_text>Replacement rules:<jupyter_code>nondim = [E==E_0*e,dEdt==(E_0/tau_0)*dedtau,I==I_0*i,dIdt==(I_0/tau_0)*didtau,W_EE==W_EE0*w_EE,W_EI==W_EI0*w_EI,W_IE==W_IE0*w_IE,W_II==W_II0*w_II,dWEEdt==(W_EE0/tau_0)*dwEEdtau,dWEIdt==(W_EI0/tau_0)*dwEIdtau,dWIEdt==(W_IE0/tau_0)*dwIEdtau,dWIIdt==(W_II0/tau_0)*dwIIdtau]
show(nondim)<jupyter_output><empty_output><jupyter_text>Full system after replacement:<jupyter_code>f_e = f_E.subs(nondim)*tau_E/E_0
f_i = f_I.subs(nondim)*tau_I/I_0
f_wEE = f_WEE.subs(nondim)*tau_0/W_EE0
f_wEI = f_WEI.subs(nondim)*tau_0/W_EI0
f_wIE = f_WIE.subs(nondim)*tau_0/W_IE0
f_wII = f_WII.subs(nondim)*tau_0/W_II0
show(f_e.expand())
show(f_i.expand())
show(f_wEE.expand())
show(f_wEI.expand())
show(f_wIE.expand())
show(f_wII.expand())<jupyter_output><empty_output><jupyter_text>Normalization of equations  Choose scaling parameters such that nondimensional variables are of order 1:  Choose scaling parameters such that nondimensional variables are of order 1:  $\tau_0 = 1/(\alpha_{EE} E_{set})$  $E_0 = E_{set}$  $I_0 = I_{set}$  $W_{EE0} = 1/g_E$  $W_{II0} = R/g_E$  $W_{EI0} = 1/(R g_I)$  $W_{IE0} = 1/g_I$  Define new parameters:  $\epsilon_E = \tau_E/\tau_0$  $\epsilon_I = \tau_I/\tau_0$  $\alpha_2 = \alpha_{EI}/\alpha_{EE}$  $\alpha_3 = \alpha_{IE}/\alpha_{EE}$  $\alpha_4 = \alpha_{II}/\alpha_{EE}$  $R = E_{set}/I_{set}$  $g = g_I/g_E$  $\theta_E = (g_E/E_{set}) Thr_E$  $\theta_I = (g_I/I_{set}) Thr_I$<jupyter_code>var('epsilon_E,epsilon_I')
var('theta_E,theta_I,R,g,alpha_2,alpha_3,alpha_4')
normaliz = [tau_0==1/(alpha_EE*E_set),tau_E==epsilon_E*tau_0,tau_I==epsilon_I*tau_0,E_0==E_set,I_0==I_set,W_EE0==1/g_E,W_EI0==R/g_E,W_IE0==1/(R*g_I),W_II0==1/g_I,I_set==E_set/R,Theta_E==theta_E*E_set/g_E,Theta_I==theta_I*I_set/g_I,g_E==g*g_I,alpha_EI==alpha_2*alpha_EE,alpha_IE==alpha_3*alpha_EE,alpha_II==alpha_4*alpha_EE]
show(normaliz)<jupyter_output><empty_output><jupyter_text>Replacement rules: numerical values, parameter definitions, nondimensionalized weights<jupyter_code>parameters = [epsilon_E==tau_E/tau_0,epsilon_I==tau_I/tau_0,g==g_E/g_I,R==E_set/I_set,theta_E==(g_E/E_set)*Theta_E,theta_I==(g_I/I_set)*Theta_I,W_EE0==1/g_E,W_EI0==R/g_E,W_IE0==1/(R*g_I),W_II0==1/g_I,alpha_2==alpha_EI/alpha_EE,alpha_3==alpha_IE/alpha_EE,alpha_4==alpha_II/alpha_EE]
weights = [w_EE==W_EE/W_EE0,w_IE==W_IE/W_IE0,w_EI==W_EI/W_EI0,w_II==W_II/W_II0]<jupyter_output><empty_output><jupyter_text>Nondimensionalized, normalized equations<jupyter_code>f_e_norm = f_e.subs(normaliz).subs(normaliz).subs(normaliz).expand()
f_i_norm = f_i.subs(normaliz).subs(normaliz).subs(normaliz).expand()
f_wEE_norm = f_wEE.subs(normaliz).subs(normaliz).subs(normaliz).expand().factor()
f_wEI_norm = f_wEI.subs(normaliz).subs(normaliz).expand().factor()
f_wIE_norm = f_wIE.subs(normaliz).subs(normaliz).expand().factor()
f_wII_norm = f_wII.subs(normaliz).subs(normaliz).subs(normaliz).expand().factor()
show(f_e_norm)
show(f_i_norm)
show(f_wEE_norm)
show(f_wEI_norm)
show(f_wIE_norm)
show(f_wII_norm)<jupyter_output><empty_output><jupyter_text>Approximation: Quasi-steady-state solutions for the neural subsystem Neural dynamics is much faster than synaptic dynamics, thus we assume E and I are in "equilibrium" (~instantaneous convergence to quasi-steady-state values, as if the weights were fixed) while the synaptic weights evolve according to their slow dynamics.  That is, we consider  $\tau_E \ll \tau_0$ and $\tau_I \ll \tau_0$  thus  $\epsilon_E \ll 1$ and $\epsilon_E \ll 1$  thus  $\epsilon_E\frac{de}{d\tau} \sim 0$ and $\epsilon_I\frac{di}{d\tau} \sim 0$:<jupyter_code>neuralFixedPoint_nondim = solve([f_e_norm.subs(dedtau==0),f_i_norm.subs(didtau==0)],e,i)
e_ss = neuralFixedPoint_nondim[0][0]
i_ss = neuralFixedPoint_nondim[0][1]
show(e_ss)
show(i_ss)<jupyter_output><empty_output><jupyter_text>Stability of the Up state  Jacobian of the nondimensionalized synaptic dynamics in the QSS approximation If $\frac{dw_{EE}}{d\tau} = f_1(w_{EE},w_{EI}...)$, ...  then the Jacobian matrix is$J_{qss} = \begin{pmatrix}\frac{\partial f_1}{\partial w_{EE}} & \frac{\partial f_1}{\partial w_{EI}} & ... \\... & ... & ...\end{pmatrix}$  Stability analysis: we should evaluate $J_{qss}$ at the QSS fixed point and compute its eigenvalues. If at least one eigenvalue has positive real part, then the fixed point is unstable.  Jacobian matrix The Jacobian matrix $J$ is the matrix of all the first-order partial derivatives of the learning rule with respect to the weights. In terms of the QSS synaptic system:  $\displaystyle \frac{dw_{EE}}{d \tau} = f_1(w_{EE},w_{EI},\ldots)$  $\displaystyle \frac{dw_{EI}}{d \tau} = f_2(w_{EE},w_{EI},\ldots)$  ... etc  then  \begin{array}{lll}J_{11} = \displaystyle \frac{\partial f_1}{\partial w_{EE}} & J_{12} = \displaystyle \frac{\partial f_1}{\partial w_{EI}} & \ldots \\J_{21} = \displaystyle \frac{\partial f_2}{\partial w_{EE}} & J_{22} = \displaystyle \frac{\partial f_2}{\partial w_{EI}} & \ldots \\\vdots & \vdots &\end{array}  Instead of substituting $e = e_{ss}(w_{EE},w_{EI},w_{IE},w_{II})$ and $i = i_{ss}(w_{EE},w_{EI},w_{IE},w_{II})$ directly into $f_1(e,i), \ldots$ and then differentiating, we will use the chain rule:  Actual computation using the chain rule The learning rule is expressed as a composition of functions $e$ and $i$ of the weights:  $\displaystyle \frac{dw_{EE}}{d \tau} = f_1(e,i)$  $\displaystyle \frac{dw_{EI}}{d \tau} = f_2(e,i)$ ... etc  Using the chain rule, the Jacobian matrix elements can be computed as  $J_{11} = \displaystyle \frac{\partial f_1}{\partial w_{EE}} = \frac{\partial f_1}{\partial e} \frac{\partial e}{\partial w_{EE}} + \frac{\partial f_1}{\partial i} \frac{\partial i}{\partial w_{EE}}$  $J_{12} = \displaystyle \frac{\partial f_1}{\partial w_{EI}} = \frac{\partial f_1}{\partial e} \frac{\partial e}{\partial w_{EI}} + \frac{\partial f_1}{\partial i} \frac{\partial i}{\partial w_{EI}}$  $J_{21} = \displaystyle \frac{\partial f_2}{\partial w_{EE}} = \frac{\partial f_2}{\partial e} \frac{\partial e}{\partial w_{EE}} + \frac{\partial f_2}{\partial i} \frac{\partial i}{\partial w_{EE}}$  $J_{22} = \displaystyle \frac{\partial f_2}{\partial w_{EI}} = \frac{\partial f_2}{\partial e} \frac{\partial e}{\partial w_{EI}} + \frac{\partial f_2}{\partial i} \frac{\partial i}{\partial w_{EI}}$ If the learning rule includes explicit dependence on the weights (as in this case) we need a more general expression:  $\frac{dw_{EE}}{d \tau} = f_1(e,i,w_{EE},w_{EI},\ldots)$  $\frac{dw_{EI}}{d \tau} = f_2(e,i,w_{EE},w_{EI},\ldots)$ ... etc  first substitute  $w_{EE} \rightarrow p_{EE}$  $w_{EI} \rightarrow p_{EI}, \ldots$  to get  $\displaystyle \frac{dw_{EE}}{d \tau} = f_1(e,i,p_{EE},p_{EI},\ldots)$  $\displaystyle \frac{dw_{EI}}{d \tau} = f_2(e,i,p_{EE},p_{EI},\ldots)$  and then take the partial derivatives:  $J_{11} = \displaystyle \frac{\partial f_1}{\partial w_{EE}} = \frac{\partial f_1}{\partial e} \frac{\partial e}{\partial w_{EE}} + \frac{\partial f_1}{\partial i} \frac{\partial i}{\partial w_{EE}} + \frac{\partial f_1}{\partial p_{EE}}$ (because $\frac{\partial p_{EE}}{\partial w_{EE}}=1$ and $\frac{\partial p_{EI}}{\partial w_{EE}}=0 \ldots$)  $J_{12} = \displaystyle \frac{\partial f_1}{\partial w_{EI}} = \frac{\partial f_1}{\partial e} \frac{\partial e}{\partial w_{EI}} + \frac{\partial f_1}{\partial i} \frac{\partial i}{\partial w_{EI}} + \frac{\partial f_1}{\partial p_{EI}}$ (because $\frac{\partial p_{EE}}{\partial w_{EI}}=0$ and $\frac{\partial p_{EI}}{\partial w_{EI}}=1 \ldots$)  ...etc<jupyter_code>var('dedwEE_ph,didwEE_ph') #placeholders
var('dedwEI_ph,didwEI_ph')
var('dedwIE_ph,didwIE_ph')
var('dedwII_ph,didwII_ph')
J_11 = diff(dwEEdtau.subs(f_wEE_norm),e)*dedwEE_ph + diff(dwEEdtau.subs(f_wEE_norm),i)*didwEE_ph + diff(dwEEdtau.subs(f_wEE_norm),w_EE)
J_12 = diff(dwEEdtau.subs(f_wEE_norm),e)*dedwEI_ph + diff(dwEEdtau.subs(f_wEE_norm),i)*didwEI_ph + diff(dwEEdtau.subs(f_wEE_norm),w_EI)
J_13 = diff(dwEEdtau.subs(f_wEE_norm),e)*dedwIE_ph + diff(dwEEdtau.subs(f_wEE_norm),i)*didwIE_ph + diff(dwEEdtau.subs(f_wEE_norm),w_IE)
J_14 = diff(dwEEdtau.subs(f_wEE_norm),e)*dedwII_ph + diff(dwEEdtau.subs(f_wEE_norm),i)*didwII_ph + diff(dwEEdtau.subs(f_wEE_norm),w_II)
J_21 = diff(dwEIdtau.subs(f_wEI_norm),e)*dedwEE_ph + diff(dwEIdtau.subs(f_wEI_norm),i)*didwEE_ph + diff(dwEIdtau.subs(f_wEI_norm),w_EE)
J_22 = diff(dwEIdtau.subs(f_wEI_norm),e)*dedwEI_ph + diff(dwEIdtau.subs(f_wEI_norm),i)*didwEI_ph + diff(dwEIdtau.subs(f_wEI_norm),w_EI)
J_23 = diff(dwEIdtau.subs(f_wEI_norm),e)*dedwIE_ph + diff(dwEIdtau.subs(f_wEI_norm),i)*didwIE_ph + diff(dwEIdtau.subs(f_wEI_norm),w_IE)
J_24 = diff(dwEIdtau.subs(f_wEI_norm),e)*dedwII_ph + diff(dwEIdtau.subs(f_wEI_norm),i)*didwII_ph + diff(dwEIdtau.subs(f_wEI_norm),w_II)
J_31 = diff(dwIEdtau.subs(f_wIE_norm),e)*dedwEE_ph + diff(dwIEdtau.subs(f_wIE_norm),i)*didwEE_ph + diff(dwIEdtau.subs(f_wIE_norm),w_EE)
J_32 = diff(dwIEdtau.subs(f_wIE_norm),e)*dedwEI_ph + diff(dwIEdtau.subs(f_wIE_norm),i)*didwEI_ph + diff(dwIEdtau.subs(f_wIE_norm),w_EI)
J_33 = diff(dwIEdtau.subs(f_wIE_norm),e)*dedwIE_ph + diff(dwIEdtau.subs(f_wIE_norm),i)*didwIE_ph + diff(dwIEdtau.subs(f_wIE_norm),w_IE)
J_34 = diff(dwIEdtau.subs(f_wIE_norm),e)*dedwII_ph + diff(dwIEdtau.subs(f_wIE_norm),i)*didwII_ph + diff(dwIEdtau.subs(f_wIE_norm),w_II)
J_41 = diff(dwIIdtau.subs(f_wII_norm),e)*dedwEE_ph + diff(dwIIdtau.subs(f_wII_norm),i)*didwEE_ph + diff(dwIIdtau.subs(f_wII_norm),w_EE)
J_42 = diff(dwIIdtau.subs(f_wII_norm),e)*dedwEI_ph + diff(dwIIdtau.subs(f_wII_norm),i)*didwEI_ph + diff(dwIIdtau.subs(f_wII_norm),w_EI)
J_43 = diff(dwIIdtau.subs(f_wII_norm),e)*dedwIE_ph + diff(dwIIdtau.subs(f_wII_norm),i)*didwIE_ph + diff(dwIIdtau.subs(f_wII_norm),w_IE)
J_44 = diff(dwIIdtau.subs(f_wII_norm),e)*dedwII_ph + diff(dwIIdtau.subs(f_wII_norm),i)*didwII_ph + diff(dwIIdtau.subs(f_wII_norm),w_II)
J = matrix([[J_11,J_12,J_13,J_14],[J_21,J_22,J_23,J_24],[J_31,J_32,J_33,J_34],[J_41,J_42,J_43,J_44]])
show(J)<jupyter_output><empty_output><jupyter_text>In order to compute the eigenvalues of the Jacobian matrix we need:  - the (quasi)steady-state values $e=e_{ss}(w_{EE},w_{EI},\ldots)$ and $i=i_{ss}(w_{EE},w_{EI},\ldots)$  - the partial derivatives $\displaystyle \frac{de}{dw_{EE}} \ldots$  - the fixed point values $w_{EEup}, w_{EIup},\ldots$  Quasi-steady-state values of $e$ and $i$<jupyter_code>show(e_ss)
show(i_ss)<jupyter_output><empty_output><jupyter_text>Partial derivatives<jupyter_code>dedwEE = diff(e.subs(e_ss),w_EE).factor()
dedwEI = diff(e.subs(e_ss),w_EI).factor()
dedwIE = diff(e.subs(e_ss),w_IE).factor()
dedwII = diff(e.subs(e_ss),w_II).factor()
didwEE = diff(i.subs(i_ss),w_EE).factor()
didwEI = diff(i.subs(i_ss),w_EI).factor()
didwIE = diff(i.subs(i_ss),w_IE).factor()
didwII = diff(i.subs(i_ss),w_II).factor()
dxdw = [dedwEE_ph==dedwEE,dedwEI_ph==dedwEI,dedwIE_ph==dedwIE,dedwII_ph==dedwII,didwEE_ph==didwEE,didwEI_ph==didwEI,didwIE_ph==didwIE,didwII_ph==didwII]<jupyter_output><empty_output><jupyter_text>Fixed points Fixed points are defined by $\frac{d w_{EE}}{d\tau}=0$, etc.  1. Fixed points without substitution: (i.e. in terms of $e$ and $i$)<jupyter_code>synapticFixedPoint_qss_nosubs = solve([f_wEE_norm.subs(dwEEdtau==0),f_wEI_norm.subs(dwEIdtau==0),f_wIE_norm.subs(dwIEdtau==0),f_wII_norm.subs(dwIIdtau==0)],e,i)
synapticFixedPoint_qss_nosubs_up = synapticFixedPoint_qss_nosubs[0]
show(synapticFixedPoint_qss_nosubs_up)<jupyter_output><empty_output><jupyter_text>2. Fixed points with substitution (i.e. explicit values of the weights).  Solve for the weight values at the Up state:<jupyter_code>synapticFixedPoint_qss_subs = solve([e_ss.subs(e==1),i_ss.subs(i==1)],w_EI,w_II)
synapticFixedPoint_qss_subs_up = synapticFixedPoint_qss_subs[0]
show(synapticFixedPoint_qss_subs_up)<jupyter_output><empty_output><jupyter_text>Interpretation of the above expressions: the weights converge to a plane attractor (a 2D plane in 4D phase space space---two weights are free and the other two are linear functions of the first two), so potentially many different combinations of weight values are a solution.  Confirm obtained solution:<jupyter_code>show(e_ss.subs(synapticFixedPoint_qss_subs_up).full_simplify())
show(i_ss.subs(synapticFixedPoint_qss_subs_up).full_simplify())<jupyter_output><empty_output><jupyter_text>Fixed-point values of the original weights at the Up state:<jupyter_code>synapticFixedPoint_qss_subs_up_orig = [W_EIup,W_IIup]
show(synapticFixedPoint_qss_subs_up_orig)<jupyter_output><empty_output><jupyter_text>Evaluation of the Jacobian matrix at the Up state fixed point The Jacobian matrix evaluated at the Up state:<jupyter_code>J_up = J.subs(synapticFixedPoint_qss_nosubs_up).subs(dxdw).subs(synapticFixedPoint_qss_subs_up).simplify()
show(J_up)<jupyter_output><empty_output><jupyter_text>The Jacobian matrix at the first solution has two dependent rows, then two of the eigenvalues will be zero:<jupyter_code>J_up.eigenvalues()<jupyter_output><empty_output><jupyter_text>Recall that the steady-state solutions of the learning rule form a 2D plane in 4D phase space, and thus two out of the four eigenvalues are expected to be zero (meaning there is no dynamics along the plane).  If both nonzero eigenvalues have negative real part then the learning rule is stable; if any of them has positive real part, then the rule is unstable.  Linear stability of the directions with nonzero eigenvalues Take the nonzero eigenvalues and check whether any of them has positive real part:<jupyter_code>J_up_eigvalues = J_up.eigenvalues()
J_up_eigvalues1_aux = J_up_eigvalues[0]
J_up_eigvalues2_aux = J_up_eigvalues[1]
show(J_up_eigvalues1_aux)
show(J_up_eigvalues2_aux)<jupyter_output><empty_output><jupyter_text>Rewrite eigenvalues Define $A,B,C$ and write the eigenvalues in terms of them:  (this section might depend on the particular implementation of Sage, as specific terms and factors are extracted from the eigenvalue expression)<jupyter_code>A = sum([J_up_eigvalues1_aux.numerator().operands()[k] for k in [0..15]]).expand().factor() # all terms in the numerator but the square root
B = J_up_eigvalues1_aux.numerator().operands()[16].operands()[0] # argument of the square root
C = J_up_eigvalues1_aux.denominator() # denominator
show(A)
show(B)
show(C)
BA2diff = (B - A^2).expand().factor()
show(BA2diff)<jupyter_output><empty_output><jupyter_text>BA2diff is equal to C up to a negative factor (see below). Define D as the opposite:<jupyter_code>D = -(BA2diff/C).expand().factor()
show(D)<jupyter_output><empty_output><jupyter_text>Show that D is positive definite by substituting $w_{EE}$ and $w_{IE}$:<jupyter_code>aux = solve([e_ss.subs(e==1),i_ss.subs(i==1)],w_EE,w_IE)
synapticFixedPoint_qss_subs_up_v2 = aux[0]
show(synapticFixedPoint_qss_subs_up_v2)
D_v2 = D.subs(synapticFixedPoint_qss_subs_up_v2).expand().factor()
show(D_v2)<jupyter_output><empty_output><jupyter_text>which is the sum of positive definite terms.  Then the eigenvalues now can be expressed as:  $$\lambda_{\pm} = \frac{\left( A \pm \sqrt{A^2 - DC} \right)}{C} \hspace{1cm} \mbox{(J_up_eigvalue_build)}$$  Confirm that the difference between J_up_eigvalue and J_up_eigvalue_build is zero:<jupyter_code>J_up_eigvalues1_build = (A + sqrt(A^2 - D*C))/C
show((J_up_eigvalues1_build - J_up_eigvalues[0]).expand())<jupyter_output><empty_output><jupyter_text>Note also that $C$ must be positive because it is related to one of the stability conditions for the neural system.  First express $C$ in terms of the original parameters:<jupyter_code>C_orig = C.subs(weights).subs(parameters).subs(parameters).expand()
show(C_orig)<jupyter_output><empty_output><jupyter_text>Now evaluate the stability condition for the neural system at the Up state:<jupyter_code>show(neural_stable_detcond_v2)
neural_stable_detcond_v3 = (neural_stable_detcond_v2*2*E_set*g_I/I_set^2).expand() # multiply by positive factors only to preserve the inequality
show(neural_stable_detcond_v3)<jupyter_output><empty_output><jupyter_text>Note that the left-hand side is equal to $C$:<jupyter_code>neural_stable_detcond_v3.lhs() - C_orig<jupyter_output><empty_output><jupyter_text>Then we know the denominator is positive because the neural system must be stable ($C>0$).  Analytic expression for the instability condition of the rule Given that $C$ must be positive and $D$ is positive, the stability of the learning rule now depends on the sign of $A$: if $A<0$ then the rule is stable because $\sqrt{A^2 - DC} < |A|$ and then both eigenvalues are negative: $A - \sqrt{A^2 - DC} < A + \sqrt{A^2 - DC} < 0$.Express $A$ in terms of $W_{EE}$ and $W_{IE}$ in order to compare it with the neural subsystem conditions:<jupyter_code>aux0 = solve(synapticFixedPoint_qss_subs_up_orig[1],W_IE)[0]
A_orig = A.subs(weights).subs(parameters).subs(parameters).expand()
A_orig_v2 = A_orig.subs(aux0).expand()
synaptic_stable_cond = (A_orig*I_set^2*alpha_EE).expand() < 0
show(synaptic_stable_cond)<jupyter_output><empty_output><jupyter_text>Analysis of the stability condition Express the synaptic stability condition in terms of $W_{EE}$ and $W_{II}$:<jupyter_code>synaptic_stable_cond_v2 = (A_orig_v2*I_set*alpha_EE).expand() < 0
synaptic_stable_cond_v3 = (I_set*alpha_II*W_II*g_I + Theta_I*g_I*alpha_IE)*(W_EE*g_E - 1) < (E_set*alpha_EE*W_EE*g_E - Theta_E*g_E*alpha_EI)*(W_II*g_I + 1) + (E_set*alpha_EI - I_set*alpha_IE)*(W_EE*g_E - 1)*(W_II*g_I + 1)
show(synaptic_stable_cond_v2)
show(synaptic_stable_cond_v3)
show(((synaptic_stable_cond_v3.lhs() - synaptic_stable_cond_v3.rhs()) - synaptic_stable_cond_v2.lhs()).expand())
show(synaptic_stable_cond_v3)<jupyter_output><empty_output><jupyter_text>The above expression can be written as:  $(W_{EE}g_E-1) a < (W_{II}g_I+1) (a^\prime + (W_{EE}g_E-1)b)$  where  $a = I_{set}W_{II} g_I \alpha_{II} + \Theta_I g_I \alpha_{IE}$  $a^\prime = E_{set}W_{EE}g_E \alpha_{EE} - \Theta_E g_E \alpha_{EI}$  $b = E_{set}\alpha_{EI} - I_{set}\alpha_{IE}$  Taking into account that  $I_{set} > E_{set}$  $g_I > g_E$  $\Theta_I > \Theta_E$  and that $a$ is the sum of positive terms while $a^\prime$ is the difference of positive terms, then for $\alpha$s of the same order we can say that  $a > a^\prime$  (despite $W_{II} \leq W_{EE}$)  and  $(W_{II}g_I+1)\tau_E > (W_{EE}g_E - 1)\tau_I$)  both helping the condition to hold but in addition it is likely that $b<0$ which has the opposite effect.  Thus it is hard to say whether the condition is satisfied or not. Case: all learning rates are equal<jupyter_code>var('alpha')
synaptic_stable_cond_v4 = (synaptic_stable_cond_v3.subs(alpha_EE==alpha,alpha_EI==alpha,alpha_IE==alpha,alpha_II==alpha)/alpha).expand()
synaptic_stable_cond_v4b = (I_set - E_set)*(W_II*g_I + 1)*(W_EE*g_E - 1) + (I_set*W_II + Theta_I)*g_I*(W_EE*g_E - 1) < (E_set*W_EE - Theta_E)*g_E*(W_II*g_I + 1)
show(synaptic_stable_cond_v4)
show(synaptic_stable_cond_v4b)
# confirm
show(((synaptic_stable_cond_v4.lhs()-synaptic_stable_cond_v4.rhs())-(synaptic_stable_cond_v4b.lhs()-synaptic_stable_cond_v4b.rhs())).expand())<jupyter_output><empty_output><jupyter_text>Semi-nondimensional version<jupyter_code>synaptic_stable_cond_v5 = (I_set*(alpha_4)*W_II*g_I + Theta_I*g_I*(alpha_3))*(W_EE*g_E - 1) < (E_set*W_EE*g_E + ((E_set*(W_EE*g_E - 1)) - Theta_E*g_E)*alpha_2 - I_set*(alpha_3)*(W_EE*g_E - 1))*(W_II*g_I + 1)
show(synaptic_stable_cond_v5)
# Confirm
show((synaptic_stable_cond_v5.subs(parameters)*alpha_EE - synaptic_stable_cond_v3).expand())<jupyter_output><empty_output><jupyter_text>If $(W_{EE}g_E - 1) > 0$, then decreasing $\alpha_4$ and/or $\alpha_3$ help satisfying the condition and thus making the rule stable.  On the other hand, the role of $\alpha_2$ is less clear.  Attempt to maximize the overlap between the stability region of the learning rule and the stability region of the neural subsystem 1. Compute the horizontal asymptote of the stability condition of the learning rule:<jupyter_code>aux = synaptic_stable_cond.lhs() - synaptic_stable_cond.rhs()
synaptic_stable_cond_border = solve(aux,W_IE)[0].factor()
synaptic_stable_cond_horiz_asymptote = synaptic_stable_cond_border.limit(W_EE=oo).subs(normaliz).factor()
show(synaptic_stable_cond_border)
show(synaptic_stable_cond_horiz_asymptote)<jupyter_output><empty_output><jupyter_text>A higher horizontal asymptote would increase the overlap, which can be achieved by increasing $\alpha_2$ and/or decreasing $\alpha_{3,4}$ 2. Additionally, we can compute the coefficient of $\alpha_3$ in the horizontal asymptote. A positive coefficient would make the asymptote to move upwards when increasing $\alpha_3$<jupyter_code>alpha_3_coeff = diff(synaptic_stable_cond_horiz_asymptote.rhs(),alpha_3)
show(alpha_3_coeff)<jupyter_output><empty_output><jupyter_text>The denominator is positive definite. The coefficient is then positive if $\alpha_2$ is large enough and/or $\alpha_4$ is small enough. 3. Simplified analysis: reduce the number of free learning rates (excitatory and inhibitory learning rates only)  $\alpha_2 \equiv \alpha_E$  $\alpha_3 = \alpha_4 \equiv \alpha_I$  The horizontal asymptote in this case is:<jupyter_code>var('alpha_E,alpha_I')
synaptic_stable_cond_horiz_asymptote_EI = synaptic_stable_cond_horiz_asymptote.subs(alpha_2==alpha_E,alpha_3==alpha_I,alpha_4==alpha_I).factor()
show(synaptic_stable_cond_horiz_asymptote_EI)<jupyter_output><empty_output><jupyter_text>which is higher if $\alpha_I$ is small enough (slow evolution of the inhibitory synapses).  Relationship between the stability condition and the paradoxical condition The stability condition for this learning rule is a homographic function in the $(W_{EE},W_{IE})$ space with horizontal and vertical asymptotes. For high $W_{IE}$ values, the region of stability is to the left of the vertical asymptote. Here we show that, when the excitatory learning rates $\alpha_{EE}$ and $\alpha_{EI}$ are zero, the vertical asymptote coincides exactly with the border of the paradoxical condition ($W_{EE}>1/g_E$), which can be construed as an inconsistency/contradiction between the stability of the rule and the existence of the paradoxical effect.<jupyter_code>aux = solve(synaptic_stable_cond_border,W_EE)[0]
synaptic_stable_cond_vert_asymptote = aux.limit(W_IE=oo).subs(normaliz).subs(parameters).subs(alpha_EE==0,alpha_EI==0).factor()
show(aux)
show(synaptic_stable_cond_vert_asymptote)<jupyter_output><empty_output><jupyter_text>Numerical analysis Stability and instability conditions in terms of the free weights $W_{EE}$ and $W_{IE}$<jupyter_code>show(neural_stable_detcond_v2)
show(neural_stable_trcond_v2)
show(paradox_cond)
show(up_exist_cond_2_v2)
show(positive_WEI_cond)
show(positive_WII_cond)
show(synaptic_stable_cond)<jupyter_output><empty_output><jupyter_text>Paradoxical conditions<jupyter_code>W_EEup = W_EE==5
W_IEup = W_IE==10
alphas = [alpha_EE==0.02,alpha_EI==0.02,alpha_IE==0.002,alpha_II==0.002] # label: alphaE=alphaI
#alphas = [alpha_EE==0.02,alpha_EI==0.02,alpha_IE==0.002,alpha_II==0.002] # label: alphaE>alphaI
#alphas = [alpha_EE==0.02,alpha_EI==0.02,alpha_IE==0.0002,alpha_II==0.0002] # label: alphaE>>alphaI
#alphas = [alpha_EE==0.002,alpha_EI==0.002,alpha_IE==0.02,alpha_II==0.02] # label: alphaE<alphaI
#alphas = [alpha_EE==0.0002,alpha_EI==0.0002,alpha_IE==0.02,alpha_II==0.02] # label: alphaE<<alphaI
#alphas = [alpha_EE==0.02,alpha_EI==0.01,alpha_IE==0.002,alpha_II==0.002] # label: alphaEE=0.02-alphaEI=0.01-alphaI=0.002
#alphas = [alpha_EE==0.02,alpha_EI==0.01,alpha_IE==0.0002,alpha_II==0.005] # label: alphaEE=0.02-alphaEI=0.01-alphaIE=0.0002-alphaII=0.005
upstate_orig = [W_EEup,
                 W_EI==W_EI.subs(synapticFixedPoint_qss_subs_up_orig).subs(values_paradoxical).subs(W_EEup),
                 W_IEup,
                 W_II==W_II.subs(synapticFixedPoint_qss_subs_up_orig).subs(values_paradoxical).subs(W_IEup)]
upstate_norm = [w_EE==w_EE.subs(weights).subs(parameters).subs(parameters).subs(values_paradoxical).subs(upstate_orig),
                 w_EI==w_EI.subs(weights).subs(parameters).subs(parameters).subs(values_paradoxical).subs(upstate_orig),
                 w_IE==w_IE.subs(weights).subs(parameters).subs(parameters).subs(values_paradoxical).subs(upstate_orig),
                 w_II==w_II.subs(weights).subs(parameters).subs(parameters).subs(values_paradoxical).subs(upstate_orig)]
show(upstate_orig)
show(upstate_norm)
probe = [W_EE==5,W_IE==10]
positive_WEI_cond_border_pdx = solve(positive_WEI_cond.lhs()==positive_WEI_cond.rhs(),W_EE)[0].subs(values_paradoxical)
positive_WII_cond_border_pdx = solve(positive_WII_cond.lhs()==positive_WII_cond.rhs(),W_IE)[0].subs(values_paradoxical)
paradox_cond_border_pdx = solve(paradox_cond.left()==paradox_cond.right(),W_EE,W_IE)[0][0].subs(values_paradoxical)
neural_stable_detcond_v2_border_pdx = solve(neural_stable_detcond_v2.lhs()==neural_stable_detcond_v2.rhs(),W_IE)[0].subs(values_paradoxical)
neural_stable_trcond_v2_border_pdx = solve(neural_stable_trcond_v2.lhs()==neural_stable_trcond_v2.rhs(),W_IE)[0].subs(values_paradoxical)
synaptic_stable_cond_border_pdx = solve(synaptic_stable_cond.lhs()==synaptic_stable_cond.rhs(),W_IE)[0].subs(values_paradoxical).subs(alphas)
up_exist_cond_2_v2_border_pdx = solve(up_exist_cond_2_v2.lhs()==up_exist_cond_2_v2.rhs(),W_IE)[0].subs(values_paradoxical)
print("PARADOXICAL CONDITIONS")
print("positive_WEI_cond:")
print("     ",positive_WEI_cond)
print("      border: ",positive_WEI_cond_border_pdx)
print("      probe: ",bool(positive_WEI_cond.subs(values_paradoxical).subs(probe)))
print("positive_WII_cond:")
print("     ",positive_WII_cond)
print("      border: ",positive_WII_cond_border_pdx)
print("      probe: ",bool(positive_WII_cond.subs(values_paradoxical).subs(probe)))
print("paradox_cond:")
print("     ",paradox_cond)
print("      border: ",paradox_cond_border_pdx)
print("      probe: ",bool(paradox_cond.subs(values_paradoxical).subs(probe)))
print("neural_stable_detcond_v2:")
print("     ",neural_stable_detcond_v2)
print("      border: ",neural_stable_detcond_v2_border_pdx)
print("      probe: ",bool(neural_stable_detcond_v2.subs(values_paradoxical).subs(probe)))
print("neural_stable_trcond_v2:")
print("     ",neural_stable_trcond_v2)
print("      border: ",neural_stable_trcond_v2_border_pdx)
print("      probe: ",bool(neural_stable_trcond_v2.subs(values_paradoxical).subs(probe)))
print("synaptic_stable_cond:")
print("     ",synaptic_stable_cond)
print("      border: ",synaptic_stable_cond_border_pdx)
print("      probe: ",bool(synaptic_stable_cond.subs(values_paradoxical).subs(alphas).subs(probe)))
print("up_exist_cond_2_v2:")
print("     ",up_exist_cond_2_v2)
print("      border: ",up_exist_cond_2_v2_border_pdx)
print("      probe: ",bool(up_exist_cond_2_v2.subs(values_paradoxical).subs(probe)))
W_EE_max = 10
fig1 = line([[W_EE.subs(positive_WEI_cond_border_pdx),0],[W_EE.subs(positive_WEI_cond_border_pdx),50]],color='blue',linestyle='-',legend_label='positive WEI')
fig2 = line([[0,W_IE.subs(positive_WII_cond_border_pdx)],[W_EE_max,W_IE.subs(positive_WII_cond_border_pdx)]],color='blue',linestyle='--',legend_label='positive WEI')
fig3 = line([[W_EE.subs(paradox_cond_border_pdx),0],[W_EE.subs(paradox_cond_border_pdx),50]],color='magenta',legend_label='paradoxical')
fig4 = plot(W_IE.subs(neural_stable_detcond_v2_border_pdx),(W_EE,0,W_EE_max),ymin=0,color='green',linestyle='-',legend_label='neural detcond')
fig5 = plot(W_IE.subs(neural_stable_trcond_v2_border_pdx),(W_EE,0,W_EE_max),ymin=0,color='green',linestyle='--',legend_label='neural trcond')
fig7 = plot(W_IE.subs(synaptic_stable_cond_border_pdx),(W_EE,0,W_EE_max),detect_poles='no',ymin=0,color='black',linestyle='-',legend_label='synaptic stable')
tt1 = text('(positive\nWEI)', (0.1+W_EE.subs(positive_WEI_cond_border_pdx),18),color='blue',horizontal_alignment='left')
tt2 = text('(positive WII)', (2.5,0.5+W_IE.subs(positive_WII_cond_border_pdx)),color='blue',horizontal_alignment='left')
tt3 = text('(paradoxical)', (0.1+W_EE.subs(paradox_cond_border_pdx),20),color='magenta',horizontal_alignment='left')
tt4 = text('(detcond\nstable)', (4,14),color='green',horizontal_alignment='left')
tt5 = text('(trcond\nstable)', (7,7.5),color='green',horizontal_alignment='left')
#tt7a = text('(synaptic stable)', (7,2),color='black',horizontal_alignment='left')
tt7b = text('(synaptic stable)', (0.6,14),color='black',horizontal_alignment='left',rotation=90)
tt7a = text('(synaptic stable)', (7,10),color='black',horizontal_alignment='left')
#tt7b = text('(synaptic stable)', (1.2,14),color='black',horizontal_alignment='left',rotation=90)
fig = fig1 + fig2 + fig3 + fig4 + fig5 + fig7 + tt1 + tt2 + tt3 + tt4 + tt5 + tt7a + tt7b
fig.xmin(0)
fig.xmax(W_EE_max)
fig.ymin(0)
fig.ymax(20)
fig.axes_labels(['$W_{EE}$', '$W_{IE}$'])
fig.set_legend_options(loc='upper right')
fig.legend(False)
#fig.save('SS_stability_alphaE=alphaI.pdf')
#fig.save('SS_stability_alphaE>alphaI.pdf')
#fig.save('SS_stability_alphaE>>alphaI.pdf')
#fig.save('SS_stability_alphaE<alphaI.pdf')
#fig.save('SS_stability_alphaE<<alphaI.pdf')
#fig.save('SS_stability_alphaEE=0.02-alphaEI=0.01-alphaI=0.002.pdf')
#fig.save('SS_stability_alphaEE=0.02-alphaEI=0.01-alphaIE=0.0002-alphaII=0.005.pdf')
show(fig)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15394266962648392, 'fraction_numerical': 0.021523019013608726, 'mean_word_length': 7.883681051729066, 'pattern_counts': {'":': 0, '<': 171, '<?xml version=': 0, '>': 168, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '38', 'dir': 'jupyter-structured-clean-dedup', 'id': '218076', 'license': 'MIT', 'n_tokens_mistral': 13205, 'n_tokens_neox': 12719, 'n_words': 2031, 'path': 'SageMath code/upstates-SynapticScaling stability.ipynb', 'repo_name': 'saraysoldado/UpDev2021-1', 'size': '31089'}