642-Embedding Flows in a Competition on a Manifold.sagews

Embedding Flows in a Competition on a Manifold

Introduction

Different dynamics are produced, seemingly by the same underlying system, in different conditions (Huys, others), motivating the question, here in neural systems, of how a homogeneous substrate, such as a sheet of firing rate neurons, can be switched to different dynamics.

Syntax Error: %hide

Example of Phase Flow Competition

For a first concrete example, we will take the first four equations 3.5 (and the competition component from 3.4) from Pillai 2008 corresponding to a double Fitzhugh-Nagumo excitator system in which the phase flow composed of $(\xi_1,\xi_2)$ competes with that of $(\xi_3,\xi_4)$:

$\dot{ \xi_1 } = ( 1-\sum_{i=1}^4 \xi_i^2) \xi_1 + \mu (\xi_2+\xi_1-\xi_1^3/3) \tau - c_1 \xi_1 \sum_i^4 \xi_i^2$

$\dot{ \xi_2 } = ( 1-\sum_{i=1}^4 \xi_i^2) \xi_2 -  \mu (\xi_1 - a_1 - I_1 - b_1 \xi_2 ) / \tau - c_1  \xi_2 \sum_i^4 \xi_i^2$

$\dot{ \xi_3 } = ( 1-\sum_{i=1}^4 \xi_i^2) \xi_3 + \mu (\xi_4+\xi_3-\xi_3^3/3) \tau - c_2 \xi_3 \sum_i^4 \xi_i^2$

$\dot{ \xi_4 } = ( 1-\sum_{i=1}^4 \xi_i^2) \xi_4 -  \mu (\xi_3 - a_2 - I_2 - b_2 \xi_4 ) / \tau - c_2 \xi_4 \sum_i^4 \xi_i^2$

in which $a_i$, $b_i$, and $I_i$ produce different phase flows, i.e. a combination of monostable, bistable or limit cycle.

Fig 1: Controlling Flows

The parameters $c_1$ and $c_2$ control the competition between the $(\xi_1,\xi_2)$ and $(\xi_3,\xi_4)$ dynamics. As $c_i$ is raised (lowered) the $i^{th}$ functional subnetwork's activity is inhibited (enhanced).

 

Syntax Error: %hide

Formalities of Functional and Neural Systems

In terms of order parameters, $\xi_i$ we have

$$\xi_i = f(\vec{\xi}) \xi_i + \mu g(\vec{\xi})$$

and in terms firing rate nodes we have

$$q_i = - q_i + \sum_j [ z_{ij} S(q_j) ] q_i$$

or

$$\vec{ q } = - \vec{ q } + [\textbf{ Z } \vec{ S } ( \vec{ q } )  ] \vec{ q }$$

where $\textbf{Z}$ is a connectivity matrix. Based on the coordinate transform $\textbf{w}$, from which we have $\vec{ q } = \textbf{ w } \vec{ \xi }$, $\vec{ \xi } = \textbf{ w }^\dagger \vec{ q }$ and an adjoint system $\textbf{w}^\dagger \textbf{w} = \textbf{I}$, we write

$$\textbf{ w } \vec{ \xi } = [ -1 + \textbf{ Z } \vec{ S } ( \textbf{ w } \vec{ \xi } )  ] \textbf{ w } \vec{ \xi }$$

$$\textbf{w}^\dagger \textbf{ w } \vec{ \xi } = [ - \textbf{w}^\dagger  + \textbf{w}^\dagger \textbf{ Z } \vec{ S } ( \textbf{ w } \vec{ \xi } )  ] \textbf{ w } \vec{ \xi }$$

and expand $\textbf{ Z }$ into $z_0 + \mu \textbf{ Z }$: 

$$\vec{ \xi } = [ - \textbf{w}^\dagger  + \textbf{w}^\dagger ( z_0 + \mu \textbf{ Z } ) \vec{ S } ( \textbf{ w } \vec{ \xi } ) ] \textbf{ w } \vec{ \xi }$$

$$\vec{ \xi } = [ - \textbf{w}^\dagger  + z_0 \textbf{w}^\dagger \vec{ S } ( \textbf{ w } \vec{ \xi } ) + \mu \textbf{w}^\dagger \textbf{ Z } \vec{ S } ( \textbf{ w } \vec{ \xi } ) ] \textbf{ w } \vec{ \xi }$$.

At this point, let the neural interaction function (taken to be a sigmoid) be approximated as $S(q) = a_1 + a_2 q + a_3 q^2$ (vector multiplication is element-wise unless otherwise denoted) and expand the form: 

$$\vec{ \xi } = [ - \textbf{w}^\dagger  + z_0 \textbf{w}^\dagger ( \vec{ a_1 } + \vec{ a_2 } \textbf{ w } \vec{ \xi } + \vec{ a_3 } ( \textbf{ w } \vec{ \xi } )^2 ) ] \textbf{ w } \vec{ \xi } + \mu \textbf{w}^\dagger \textbf{ Z } [ \vec{ a_1 } + \vec{ a_2 } \textbf{ w } \vec{ \xi } + \vec{ a_3 } ( \textbf{ w } \vec{ \xi } )^2 ] \textbf{ w } \vec{ \xi }$$

In the formulation of functional system dynamics above, flows are a static property of the vector field and each flow has a corresponding coefficient that determines its role in competition with other flows. In the neural network formulation, asymmetries in the connectivity matrix $\textbf{Z}$ produce flows. In order to recruit multiple flows in a neural network, we first assume that, as in the functional form, the capability to produce a particular flow or set of flows is a static property of the vector field, i.e. the connectivity matrix does not change. The neural network analog of the biased competition present in the equations used next is a linear combination of baseline excitability, $I_{j,i}$ in the $i^{th}$ neuron,  corresponding to the $j^{th}$ phase flow:

$$\dot{ q }_i = - ( c_1 I_{1,i} + c_2 I_{2,i} ) - q_i + q_i \sum_j z_{ij} S( q_j )$$

which after the derivation above takes the form

$$\dot{\vec{\xi}} = f(\vec{\xi}) \vec{\xi} + \mu g(\vec{\xi}) - \textbf{w}^\dagger ( c_1 \vec{I_1} + c_2 \vec{I_2} )$$

where $f(\vec{\xi})=\textbf{w}^\dagger [ - 1  + z_0 ( a_1 + a_2 \textbf{ w } \vec{ \xi } + a_3 ( \textbf{ w } \vec{ \xi } )^2 ) ] \textbf{ w }$ and $g(\vec{\xi})= \textbf{w}^\dagger \textbf{ Z } [ a_1 +  a_2 \textbf{ w } \vec{ \xi } + a_3 (\textbf{ w } \vec{ \xi } )^2 ] \textbf{ w } \vec{ \xi }$. Comparing this form to that of the functional dynamics

$$\dot{\vec{\xi}} = f(\vec{\xi}) \vec{\xi} + \mu g(\vec{\xi}) - \vec{c} \vec{\xi} \sum_i \xi_i^2$$

we should be able to fit $f(.)$ from the network formulation to that of the functional form in order to create the desired manifold and similarly $g(.)+ \dots$ to produce the flows and flow competition. In the procedure to fit the network formulation, we want to minimize the angle between the functional vector field, $\vec{V}_f (\vec{\xi},\vec{c})$, and the network vector field, $\vec{V}_n (\vec{\xi},\vec{c})$,

$$\theta = \| \arccos \frac{ \langle \vec{V}_f, \vec{V}_n \rangle } {\| \vec{V}_f \| \| \vec{V}_n \| } \|$$

($\arccos$ is evaluated element-wise) evaluated over some region of the space spanned by $(\vec{\xi}, \vec{c})$ by manipulating parameters of the network vector field $\textbf{Z}$, $z_0$, $a_i$, and $I_{j,i}$. In particular, $z_0$ and $a_i$ will be set when fitting the manifold, while $\textbf{Z}$ and $I_{j,i}$ will be set in order to produce the flows and their competition.

Syntax Error: %hide

Embedding a Single Flow

First, we'll attempt to embed a simpler excitator limit-cycle flow into a network of five firing rate neurons. From the above, we have

$$\dot{ \xi_1 } = ( 1- \xi_1^2 - \xi_2^2) \xi_1 + \mu (\xi_2+\xi_1-\xi_1^3/3) \tau$$

$$\dot{ \xi_2 } = ( 1-  \xi_1^2 - \xi_2^2 ) \xi_2 -  \mu (\xi_1 - a_1 - I_1 - b_1 \xi_2 ) / \tau$$

whose behavior we want to map into

$$\dot{ q_i } = I_i - q_i + q_i \sum_{j=1}^{5} z_{ij} S( q_j )$$

by fitting the form $\dot{\vec{\xi}} = f(\vec{\xi}) \vec{\xi} + \mu g(\vec{\xi})$ (where $f(\vec{\xi})=\textbf{w}^\dagger [ - 1  + z_0 ( a_1 + a_2 \textbf{ w } \vec{ \xi } + a_3 ( \textbf{ w } \vec{ \xi } )^2 ) ] \textbf{ w }$ and $g(\vec{\xi})= \textbf{w}^\dagger \textbf{ Z } [ a_1 +  a_2 \textbf{ w } \vec{ \xi } + a_3 (\textbf{ w } \vec{ \xi } )^2 ] \textbf{ w } \vec{ \xi }$) as shown above. In the following we put subscripts $\xi$ or $q$ to denote the form from which the expression is taken, i.e. $f_\xi (.)$ vs $f_q(.)$.

First, we select a coordinate transform $\textbf{w}$ and one of its possible adjoints:

# matrix with dimensions 5x2
w = matrix( [[1,0], [0,1], [0,0], [0,0], [0,0]] )

# and an adjoint
w_adj = matrix( [[1,0,1,1,0], [0,1,1,1,1]] )
w_adj * w

\newcommand{\Bold}[1]{\mathbf{#1}}\left($$\begin{array}{rr} 1 & 0 \\ 0 & 1 \end{array}$$\right)

Creating the Manifold

We describe the functional manifold as

$$f_\xi (\vec{\xi}) = 1-\xi_1^2-\xi_2^2$$

# create symbolic vectors for q and xi space
xi = mkSymVec('xi',2)
q = mkSymVec('q',5)

# parameters
var('z0 a1')
#a1 = mkSymVec('a1',5)
a2 = mkSymVec('a2',5)
a3 = mkSymVec('a3',5)

# xi manifold
f_xi = vector([(1-xi[0]**2-xi[1]**2)*xi[0],xi[1]*(1-xi[0]**2-xi[1]**2)])
f_xi

\newcommand{\Bold}[1]{\mathbf{#1}}\left(-{(\xi_{1}^{2} + \xi_{2}^{2} - 1)} \xi_{1},-{(\xi_{1}^{2} + \xi_{2}^{2} - 1)} \xi_{2}\right)

 

and the neural network manifold

$$f_q (\vec{\xi}) = \textbf{w}^\dagger [ - 1  + z_0 ( a_1 + \vec{a}_2 \textbf{ w } \vec{ \xi } + \vec{a}_3 ( \textbf{ w } \vec{ \xi } )^2 ) ]$$

# q manifold
f_q = []
xi_sq = (w*xi).pairwise_product(w*xi)
for i in range(0,5):
    f_q.append( (-1 + z0*( a1 + a2.dot_product(w*xi) + a3.dot_product(xi_sq)))*((w*xi)[i]))
f_q = w_adj*vector(f_q)
f_q

\newcommand{\Bold}[1]{\mathbf{#1}}\left({({(a_{31} \xi_{1}^{2} + a_{32} \xi_{2}^{2} + a_{21} \xi_{1} + a_{22} \xi_{2} + a_{1})} z_{0} - 1)} \xi_{1},{({(a_{31} \xi_{1}^{2} + a_{32} \xi_{2}^{2} + a_{21} \xi_{1} + a_{22} \xi_{2} + a_{1})} z_{0} - 1)} \xi_{2}\right)

Next, we take the angles between the corresponding components of the two vector fields, and sum the square of their errors,

$$\Theta = \sum_i^2 \theta_i^2 = \sum_i^2 [ \frac{ \int_\xi f_{q,i} f_{\xi,i} d\xi }{ \sqrt{\int_\xi f_{q,i}^2 d\xi} \sqrt{ \int_\xi f_{\xi,i}^2 d\xi} } ]^2$$

where $d\xi$ is a relevant region of the phase space. Elsewhere (Pillai 2008, Huys), fitting has done by numerically evaluating the vector field and finding a numerical angle between high-dimensional vectors. However, because the expressions for the vector fields used here are known, we can integrate to obtain an expression of the sum squared angle in terms of the fitting parameters:  (click %hide to see code)

<html><span class="math">\newcommand{\Bold}[1]{\mathbf{#1}}\arccos\left(-\frac{2}{3147} \, \frac{{({(574 \, a_{1} - 945 \, a_{22} + 1064 \, a_{31} + 1608 \, a_{32})} z_{0} + {(574 \, a_{1} + 945 \, a_{22} + 1064 \, a_{31} + 1608 \, a_{32})} z_{0} - 1148)} \sqrt{\frac{1049}{35

\sqrt{\frac{32}{315} \, {(140 \, {(2 \, a_{1} + 3 \, a_{22})} a_{31} + 168 \, {(3 \, a_{1} + 5 \, a_{22} + 4 \, a_{31})} a_{32} + 105 \, a_{1}^{2} + 315 \, a_{1} a_{22} + 140 \, a_{21}^{2} + 252 \, a_{22}^{2} + 336 \, a_{31}^{2} + 720 \, a_{32}^{2})} z_{0}^{2} + \frac{32}{315} \, {(140 \, {(2 \, a_{1} - 3 \, a_{22})} a_{31} + 168 \, {(3 \, a_{1} - 5 \, a_{22} + 4 \, a_{31})} a_{32} + 105 \, a_{1}^{2} - 315 \, a_{1} a_{22} + 140 \, a_{21}^{2} + 252 \, a_{22}^{2} + 336 \, a_{31}^{2} + 720 \, a_{32}^{2})} z_{0}^{2} - \frac{32}{45} \, {(30 \, a_{1} - 45 \, a_{22} + 40 \, a_{31} + 72 \, a_{32})} z_{0} - \frac{32}{45} \, {(30 \, a_{1} + 45 \, a_{22} + 40 \, a_{31} + 72 \, a_{32})} z_{0} + \frac{64}{3}}}\right)^{2} + \arccos\left(-\frac{4}{3147} \, \frac{{({(287 \, a_{1} - 357 \, a_{22} + 804 \, a_{31} + 532 \, a_{32})} z_{0} + {(287 \, a_{1} + 357 \, a_{22} + 804 \, a_{31} + 532 \, a_{32})} z_{0} - 574)} \sqrt{\frac{1049}{35}}}{\sqrt{\frac{32}{315} \, {(504 \, {(a_{1} + a_{22})} a_{31} + 28 \, {(10 \, a_{1} + 15 \, a_{22} + 24 \, a_{31})} a_{32} + 105 \, a_{1}^{2} + 210 \, a_{1} a_{22} + 252 \, a_{21}^{2} + 140 \, a_{22}^{2} + 720 \, a_{31}^{2} + 336 \, a_{32}^{2})} z_{0}^{2} + \frac{32}{315} \, {(504 \, {(a_{1} - a_{22})} a_{31} + 28 \, {(10 \, a_{1} - 15 \, a_{22} + 24 \, a_{31})} a_{32} + 105 \, a_{1}^{2} - 210 \, a_{1} a_{22} + 252 \, a_{21}^{2} + 140 \, a_{22}^{2} + 720 \, a_{31}^{2} + 336 \, a_{32}^{2})} z_{0}^{2} - \frac{64}{45} \, {(15 \, a_{1} - 15 \, a_{22} + 36 \, a_{31} + 20 \, a_{32})} z_{0} - \frac{64}{45} \, {(15 \, a_{1} + 15 \, a_{22} + 36 \, a_{31} + 20 \, a_{32})} z_{0} + \frac{64}{3}}}\right)^{2} }}}

Next, we can minimize the expression for the sum squared angle to a value near zero (using a guess for starting):

tmin = list(minimize(f_Theta,[0,0,0,0,1,1]))

Optimization terminated successfully. Current function value: 0.000000 Iterations: 19 Function evaluations: 22 Gradient evaluations: 22

Just to check, we evalute the expression at the parameters found:

print 'f_Theta = %g' % f_Theta.subs(a1=tmin[0],a21=tmin[1],a22=tmin[2],a31=tmin[3],a32=tmin[4],z0=tmin[5]).numerical_approx()

f_Theta = 7.80465e-12

Now, we verify the foregoing method works by sampling the two vector fields over a region of space, and compute the angle of the high-dimensional vector (click %hide to see code).

Theta = 2.5922e-06

Close enough.

Creating the Flow

The next step is to fit

$$g_q(\vec{\xi}) = \textbf{w}^\dagger \textbf{ Z } [ a_1 +  \vec{a}_2 \textbf{ w } \vec{ \xi } + \vec{a}_3 (\textbf{ w } \vec{ \xi } )^2 ] \textbf{ w } \vec{ \xi }$$

to

$$g_\xi(\vec{\xi})=\begin{pmatrix} (\xi_2+\xi_1-\xi_1^3/3) \tau  \\ (\xi_1 - \alpha - \beta \xi_2 ) / \tau  \end{pmatrix}$$

As before, we create the functional flow

var('alpha beta tau')
g_xi = vector( [ (xi2+xi1-(xi1**3)/3)/tau, -(xi1-alpha-beta*xi2)/tau ])
g_xi

\newcommand{\Bold}[1]{\mathbf{#1}}\left(-\frac{1}{3} \, \frac{{(\xi_{1}^{3} - 3 \, \xi_{1} - 3 \, \xi_{2})}}{\tau},\frac{{(\beta \xi_{2} + \alpha - \xi_{1})}}{\tau}\right)

and the neural flow

Z = mkSymMat('z',5,5)
B = mkSymVec('b',5)
g_q = []
for i in range(0,5):
    g_q.append( a1 + a2.dot_product(w*xi) + a3.dot_product( (w*xi).pairwise_product(w*xi) ) )
g_q = w_adj * Z * vector(g_q).pairwise_product(w*xi) + w_adj*B
g_q

\newcommand{\Bold}[1]{\mathbf{#1}}\left({(z_{12} + z_{32} + z_{42})} {(a_{31} \xi_{1}^{2} + a_{32} \xi_{2}^{2} + a_{21} \xi_{1} + a_{22} \xi_{2} + a_{1})} \xi_{2} + {(z_{11} + z_{31} + z_{41})} {(a_{31} \xi_{1}^{2} + a_{32} \xi_{2}^{2} + a_{21} \xi_{1} + a_{22} \xi_{2} + a_{1})} \xi_{1} + b_{1} + b_{3} + b_{4},{(z_{22} + z_{32} + z_{42} + z_{52})} {(a_{31} \xi_{1}^{2} + a_{32} \xi_{2}^{2} + a_{21} \xi_{1} + a_{22} \xi_{2} + a_{1})} \xi_{2} + {(z_{21} + z_{31} + z_{41} + z_{51})} {(a_{31} \xi_{1}^{2} + a_{32} \xi_{2}^{2} + a_{21} \xi_{1} + a_{22} \xi_{2} + a_{1})} \xi_{1} + b_{2} + b_{3} + b_{4} + b_{5}\right)

substitute our previous values for $a_{ij}$ and set $\alpha$, $\beta$ and $\tau$ in the functional flow,

g_q_ = map( lambda f: f.subs(a1=tmin[0],a21=tmin[1],a22=tmin[2],a31=tmin[3],a32=tmin[4]), g_q )
g_xi_ = map( lambda f: f.subs(alpha=.2,beta=2,tau=3), g_xi)

compute $\Theta$,

g_theta = map( lambda f,g: ftheta(f,g,xi,a,b), g_q_, g_xi_ )
g_Theta = sum( map( lambda f: f**2, g_theta ) )

minimize using free variables

g_varlist = g_Theta.variables()
g_varlist, len(g_varlist)

\newcommand{\Bold}[1]{\mathbf{#1}}((b1, b2, b3, b4, b5, z11, z12, z21, z22, z31, z32, z41, z42, z51, z52), 15)

g_ic = []
for i in range(0,len(g_varlist)):
    g_ic.append(-.2)
minimize(g_Theta,g_ic)

Optimization terminated successfully. Current function value: 0.713578 Iterations: 44 Function evaluations: 45 Gradient evaluations: 45

\newcommand{\Bold}[1]{\mathbf{#1}}\left(-2.58441011722,3.67728797475,1.29287785753,1.29287785753,3.67728797475,-1.66308902101,-3.62578243786,1.93102203328,-0.686072587586,0.467933012272,-4.11185502545,0.467933012272,-4.11185502545,1.93102203328,-0.686072587586\right)

var('t1 t2 t3')
eq1 = g_q_[0].expand() == g_xi_[0]
eq1.subs(z12=0,z32=0,z31=0,z41=0)

\newcommand{\Bold}[1]{\mathbf{#1}}-2.88384694956 \, \xi_{1}^{3} z_{11} - 2.88384694956 \, \xi_{1}^{2} \xi_{2} z_{42} - 2.883849661 \, \xi_{1} \xi_{2}^{2} z_{11} - 2.883849661 \, \xi_{2}^{3} z_{42} + 3.16144543614 \, \xi_{1} z_{11} + 3.16144543614 \, \xi_{2} z_{42} = -\frac{1}{9} \, \xi_{1}^{3} + \frac{1}{3} \, \xi_{1} + \frac{1}{3} \, \xi_{2}

g_q_[0].variables()

\newcommand{\Bold}[1]{\mathbf{#1}}\left(\xi_{1}, \xi_{2}, z_{11}, z_{12}, z_{31}, z_{32}, z_{41}, z_{42}\right)

Most of that was wrong, so Now, I'm gonna try again to translate the neural form to the functional form

N = 4
M = 2

q = mkSymVec('q',N)
xi = mkSymVec('xi',M)

w = mkSymMat('w',N,M)
wi = mkSymMat('wi',M,N) # inverse
Z = mkSymMat('z',N,N)
a = mkSymVec('a',N)
b = mkSymVec('b',N)
c = mkSymVec('c',N)
I = mkSymVec('I',N)

# ddt q_i = q_i ( -1 + sum_{j=1}^N ( z_{ij} S( q_j ) ) )
input_vec = Z*(a+b.pairwise_product(q)+c.pairwise_product(q.pairwise_product(q))).transpose()
vq = []
for i in range(0,N):
    vq.append( I[i] + q[i]*(-1 + input_vec[i][0]) )
vector(vq).transpose()

[((c4*q4^2 + b4*q4 + a4)*z14 + (c3*q3^2 + b3*q3 + a3)*z13 + (c2*q2^2 + b2*q2 + a2)*z12 + (c1*q1^2 + b1*q1 + a1)*z11 - 1)*q1 + I1] [((c4*q4^2 + b4*q4 + a4)*z24 + (c3*q3^2 + b3*q3 + a3)*z23 + (c2*q2^2 + b2*q2 + a2)*z22 + (c1*q1^2 + b1*q1 + a1)*z21 - 1)*q2 + I2] [((c4*q4^2 + b4*q4 + a4)*z34 + (c3*q3^2 + b3*q3 + a3)*z33 + (c2*q2^2 + b2*q2 + a2)*z32 + (c1*q1^2 + b1*q1 + a1)*z31 - 1)*q3 + I3] [((c4*q4^2 + b4*q4 + a4)*z44 + (c3*q3^2 + b3*q3 + a3)*z43 + (c2*q2^2 + b2*q2 + a2)*z42 + (c1*q1^2 + b1*q1 + a1)*z41 - 1)*q4 + I4]

qsub = w*xi.transpose(); qsub

\newcommand{\Bold}[1]{\mathbf{#1}}\left($$\begin{array}{r} w_{11} \xi_{1} + w_{12} \xi_{2} \\ w_{21} \xi_{1} + w_{22} \xi_{2} \\ w_{31} \xi_{1} + w_{32} \xi_{2} \\ w_{41} \xi_{1} + w_{42} \xi_{2} \end{array}$$\right)

vxiq = []
for i in range(0,N):
    vxiq.append( vq[i].subs( q1=qsub[0][0], q2=qsub[1][0], q3=qsub[2][0], q4=qsub[3][0] ) )
vector(vxiq).transpose()

\newcommand{\Bold}[1]{\mathbf{#1}}\left($$\begin{array}{r} {(w_{11} \xi_{1} + w_{12} \xi_{2})} {({({(w_{11} \xi_{1} + w_{12} \xi_{2})}^{2} c_{1} + {(w_{11} \xi_{1} + w_{12} \xi_{2})} b_{1} + a_{1})} z_{11} + {({(w_{21} \xi_{1} + w_{22} \xi_{2})}^{2} c_{2} + {(w_{21} \xi_{1} + w_{22} \xi_{2})} b_{2} + a_{2})} z_{12} + {({(w_{31} \xi_{1} + w_{32} \xi_{2})}^{2} c_{3} + {(w_{31} \xi_{1} + w_{32} \xi_{2})} b_{3} + a_{3})} z_{13} + {({(w_{41} \xi_{1} + w_{42} \xi_{2})}^{2} c_{4} + {(w_{41} \xi_{1} + w_{42} \xi_{2})} b_{4} + a_{4})} z_{14} - 1)} + I_{1} \\ {(w_{21} \xi_{1} + w_{22} \xi_{2})} {({({(w_{11} \xi_{1} + w_{12} \xi_{2})}^{2} c_{1} + {(w_{11} \xi_{1} + w_{12} \xi_{2})} b_{1} + a_{1})} z_{21} + {({(w_{21} \xi_{1} + w_{22} \xi_{2})}^{2} c_{2} + {(w_{21} \xi_{1} + w_{22} \xi_{2})} b_{2} + a_{2})} z_{22} + {({(w_{31} \xi_{1} + w_{32} \xi_{2})}^{2} c_{3} + {(w_{31} \xi_{1} + w_{32} \xi_{2})} b_{3} + a_{3})} z_{23} + {({(w_{41} \xi_{1} + w_{42} \xi_{2})}^{2} c_{4} + {(w_{41} \xi_{1} + w_{42} \xi_{2})} b_{4} + a_{4})} z_{24} - 1)} + I_{2} \\ {(w_{31} \xi_{1} + w_{32} \xi_{2})} {({({(w_{11} \xi_{1} + w_{12} \xi_{2})}^{2} c_{1} + {(w_{11} \xi_{1} + w_{12} \xi_{2})} b_{1} + a_{1})} z_{31} + {({(w_{21} \xi_{1} + w_{22} \xi_{2})}^{2} c_{2} + {(w_{21} \xi_{1} + w_{22} \xi_{2})} b_{2} + a_{2})} z_{32} + {({(w_{31} \xi_{1} + w_{32} \xi_{2})}^{2} c_{3} + {(w_{31} \xi_{1} + w_{32} \xi_{2})} b_{3} + a_{3})} z_{33} + {({(w_{41} \xi_{1} + w_{42} \xi_{2})}^{2} c_{4} + {(w_{41} \xi_{1} + w_{42} \xi_{2})} b_{4} + a_{4})} z_{34} - 1)} + I_{3} \\ {(w_{41} \xi_{1} + w_{42} \xi_{2})} {({({(w_{11} \xi_{1} + w_{12} \xi_{2})}^{2} c_{1} + {(w_{11} \xi_{1} + w_{12} \xi_{2})} b_{1} + a_{1})} z_{41} + {({(w_{21} \xi_{1} + w_{22} \xi_{2})}^{2} c_{2} + {(w_{21} \xi_{1} + w_{22} \xi_{2})} b_{2} + a_{2})} z_{42} + {({(w_{31} \xi_{1} + w_{32} \xi_{2})}^{2} c_{3} + {(w_{31} \xi_{1} + w_{32} \xi_{2})} b_{3} + a_{3})} z_{43} + {({(w_{41} \xi_{1} + w_{42} \xi_{2})}^{2} c_{4} + {(w_{41} \xi_{1} + w_{42} \xi_{2})} b_{4} + a_{4})} z_{44} - 1)} + I_{4} \end{array}$$\right)

vxi = wi*(vector(vxiq)).transpose(); vxi

\newcommand{\Bold}[1]{\mathbf{#1}}\left($$\begin{array}{r} {({(w_{11} \xi_{1} + w_{12} \xi_{2})} {({({(w_{11} \xi_{1} + w_{12} \xi_{2})}^{2} c_{1} + {(w_{11} \xi_{1} + w_{12} \xi_{2})} b_{1} + a_{1})} z_{11} + {({(w_{21} \xi_{1} + w_{22} \xi_{2})}^{2} c_{2} + {(w_{21} \xi_{1} + w_{22} \xi_{2})} b_{2} + a_{2})} z_{12} + {({(w_{31} \xi_{1} + w_{32} \xi_{2})}^{2} c_{3} + {(w_{31} \xi_{1} + w_{32} \xi_{2})} b_{3} + a_{3})} z_{13} + {({(w_{41} \xi_{1} + w_{42} \xi_{2})}^{2} c_{4} + {(w_{41} \xi_{1} + w_{42} \xi_{2})} b_{4} + a_{4})} z_{14} - 1)} + I_{1})} \mbox{wi}_{11} + {({(w_{21} \xi_{1} + w_{22} \xi_{2})} {({({(w_{11} \xi_{1} + w_{12} \xi_{2})}^{2} c_{1} + {(w_{11} \xi_{1} + w_{12} \xi_{2})} b_{1} + a_{1})} z_{21} + {({(w_{21} \xi_{1} + w_{22} \xi_{2})}^{2} c_{2} + {(w_{21} \xi_{1} + w_{22} \xi_{2})} b_{2} + a_{2})} z_{22} + {({(w_{31} \xi_{1} + w_{32} \xi_{2})}^{2} c_{3} + {(w_{31} \xi_{1} + w_{32} \xi_{2})} b_{3} + a_{3})} z_{23} + {({(w_{41} \xi_{1} + w_{42} \xi_{2})}^{2} c_{4} + {(w_{41} \xi_{1} + w_{42} \xi_{2})} b_{4} + a_{4})} z_{24} - 1)} + I_{2})} \mbox{wi}_{12} + {({(w_{31} \xi_{1} + w_{32} \xi_{2})} {({({(w_{11} \xi_{1} + w_{12} \xi_{2})}^{2} c_{1} + {(w_{11} \xi_{1} + w_{12} \xi_{2})} b_{1} + a_{1})} z_{31} + {({(w_{21} \xi_{1} + w_{22} \xi_{2})}^{2} c_{2} + {(w_{21} \xi_{1} + w_{22} \xi_{2})} b_{2} + a_{2})} z_{32} + {({(w_{31} \xi_{1} + w_{32} \xi_{2})}^{2} c_{3} + {(w_{31} \xi_{1} + w_{32} \xi_{2})} b_{3} + a_{3})} z_{33} + {({(w_{41} \xi_{1} + w_{42} \xi_{2})}^{2} c_{4} + {(w_{41} \xi_{1} + w_{42} \xi_{2})} b_{4} + a_{4})} z_{34} - 1)} + I_{3})} \mbox{wi}_{13} + {({(w_{41} \xi_{1} + w_{42} \xi_{2})} {({({(w_{11} \xi_{1} + w_{12} \xi_{2})}^{2} c_{1} + {(w_{11} \xi_{1} + w_{12} \xi_{2})} b_{1} + a_{1})} z_{41} + {({(w_{21} \xi_{1} + w_{22} \xi_{2})}^{2} c_{2} + {(w_{21} \xi_{1} + w_{22} \xi_{2})} b_{2} + a_{2})} z_{42} + {({(w_{31} \xi_{1} + w_{32} \xi_{2})}^{2} c_{3} + {(w_{31} \xi_{1} + w_{32} \xi_{2})} b_{3} + a_{3})} z_{43} + {({(w_{41} \xi_{1} + w_{42} \xi_{2})}^{2} c_{4} + {(w_{41} \xi_{1} + w_{42} \xi_{2})} b_{4} + a_{4})} z_{44} - 1)} + I_{4})} \mbox{wi}_{14} \\ {({(w_{11} \xi_{1} + w_{12} \xi_{2})} {({({(w_{11} \xi_{1} + w_{12} \xi_{2})}^{2} c_{1} + {(w_{11} \xi_{1} + w_{12} \xi_{2})} b_{1} + a_{1})} z_{11} + {({(w_{21} \xi_{1} + w_{22} \xi_{2})}^{2} c_{2} + {(w_{21} \xi_{1} + w_{22} \xi_{2})} b_{2} + a_{2})} z_{12} + {({(w_{31} \xi_{1} + w_{32} \xi_{2})}^{2} c_{3} + {(w_{31} \xi_{1} + w_{32} \xi_{2})} b_{3} + a_{3})} z_{13} + {({(w_{41} \xi_{1} + w_{42} \xi_{2})}^{2} c_{4} + {(w_{41} \xi_{1} + w_{42} \xi_{2})} b_{4} + a_{4})} z_{14} - 1)} + I_{1})} \mbox{wi}_{21} + {({(w_{21} \xi_{1} + w_{22} \xi_{2})} {({({(w_{11} \xi_{1} + w_{12} \xi_{2})}^{2} c_{1} + {(w_{11} \xi_{1} + w_{12} \xi_{2})} b_{1} + a_{1})} z_{21} + {({(w_{21} \xi_{1} + w_{22} \xi_{2})}^{2} c_{2} + {(w_{21} \xi_{1} + w_{22} \xi_{2})} b_{2} + a_{2})} z_{22} + {({(w_{31} \xi_{1} + w_{32} \xi_{2})}^{2} c_{3} + {(w_{31} \xi_{1} + w_{32} \xi_{2})} b_{3} + a_{3})} z_{23} + {({(w_{41} \xi_{1} + w_{42} \xi_{2})}^{2} c_{4} + {(w_{41} \xi_{1} + w_{42} \xi_{2})} b_{4} + a_{4})} z_{24} - 1)} + I_{2})} \mbox{wi}_{22} + {({(w_{31} \xi_{1} + w_{32} \xi_{2})} {({({(w_{11} \xi_{1} + w_{12} \xi_{2})}^{2} c_{1} + {(w_{11} \xi_{1} + w_{12} \xi_{2})} b_{1} + a_{1})} z_{31} + {({(w_{21} \xi_{1} + w_{22} \xi_{2})}^{2} c_{2} + {(w_{21} \xi_{1} + w_{22} \xi_{2})} b_{2} + a_{2})} z_{32} + {({(w_{31} \xi_{1} + w_{32} \xi_{2})}^{2} c_{3} + {(w_{31} \xi_{1} + w_{32} \xi_{2})} b_{3} + a_{3})} z_{33} + {({(w_{41} \xi_{1} + w_{42} \xi_{2})}^{2} c_{4} + {(w_{41} \xi_{1} + w_{42} \xi_{2})} b_{4} + a_{4})} z_{34} - 1)} + I_{3})} \mbox{wi}_{23} + {({(w_{41} \xi_{1} + w_{42} \xi_{2})} {({({(w_{11} \xi_{1} + w_{12} \xi_{2})}^{2} c_{1} + {(w_{11} \xi_{1} + w_{12} \xi_{2})} b_{1} + a_{1})} z_{41} + {({(w_{21} \xi_{1} + w_{22} \xi_{2})}^{2} c_{2} + {(w_{21} \xi_{1} + w_{22} \xi_{2})} b_{2} + a_{2})} z_{42} + {({(w_{31} \xi_{1} + w_{32} \xi_{2})}^{2} c_{3} + {(w_{31} \xi_{1} + w_{32} \xi_{2})} b_{3} + a_{3})} z_{43} + {({(w_{41} \xi_{1} + w_{42} \xi_{2})}^{2} c_{4} + {(w_{41} \xi_{1} + w_{42} \xi_{2})} b_{4} + a_{4})} z_{44} - 1)} + I_{4})} \mbox{wi}_{24} \end{array}$$\right)

vxis = []
for i in range(0,M):
    vxis.append( vxi[i][0].subs( b1=0, b2=0, b3=0, b4=0) )
vector(vxis).transpose()

\newcommand{\Bold}[1]{\mathbf{#1}}\left($$\begin{array}{r} {({(w_{11} \xi_{1} + w_{12} \xi_{2})} {({({(w_{11} \xi_{1} + w_{12} \xi_{2})}^{2} c_{1} + a_{1})} z_{11} + {({(w_{21} \xi_{1} + w_{22} \xi_{2})}^{2} c_{2} + a_{2})} z_{12} + {({(w_{31} \xi_{1} + w_{32} \xi_{2})}^{2} c_{3} + a_{3})} z_{13} + {({(w_{41} \xi_{1} + w_{42} \xi_{2})}^{2} c_{4} + a_{4})} z_{14} - 1)} + I_{1})} \mbox{wi}_{11} + {({(w_{21} \xi_{1} + w_{22} \xi_{2})} {({({(w_{11} \xi_{1} + w_{12} \xi_{2})}^{2} c_{1} + a_{1})} z_{21} + {({(w_{21} \xi_{1} + w_{22} \xi_{2})}^{2} c_{2} + a_{2})} z_{22} + {({(w_{31} \xi_{1} + w_{32} \xi_{2})}^{2} c_{3} + a_{3})} z_{23} + {({(w_{41} \xi_{1} + w_{42} \xi_{2})}^{2} c_{4} + a_{4})} z_{24} - 1)} + I_{2})} \mbox{wi}_{12} + {({(w_{31} \xi_{1} + w_{32} \xi_{2})} {({({(w_{11} \xi_{1} + w_{12} \xi_{2})}^{2} c_{1} + a_{1})} z_{31} + {({(w_{21} \xi_{1} + w_{22} \xi_{2})}^{2} c_{2} + a_{2})} z_{32} + {({(w_{31} \xi_{1} + w_{32} \xi_{2})}^{2} c_{3} + a_{3})} z_{33} + {({(w_{41} \xi_{1} + w_{42} \xi_{2})}^{2} c_{4} + a_{4})} z_{34} - 1)} + I_{3})} \mbox{wi}_{13} + {({(w_{41} \xi_{1} + w_{42} \xi_{2})} {({({(w_{11} \xi_{1} + w_{12} \xi_{2})}^{2} c_{1} + a_{1})} z_{41} + {({(w_{21} \xi_{1} + w_{22} \xi_{2})}^{2} c_{2} + a_{2})} z_{42} + {({(w_{31} \xi_{1} + w_{32} \xi_{2})}^{2} c_{3} + a_{3})} z_{43} + {({(w_{41} \xi_{1} + w_{42} \xi_{2})}^{2} c_{4} + a_{4})} z_{44} - 1)} + I_{4})} \mbox{wi}_{14} \\ {({(w_{11} \xi_{1} + w_{12} \xi_{2})} {({({(w_{11} \xi_{1} + w_{12} \xi_{2})}^{2} c_{1} + a_{1})} z_{11} + {({(w_{21} \xi_{1} + w_{22} \xi_{2})}^{2} c_{2} + a_{2})} z_{12} + {({(w_{31} \xi_{1} + w_{32} \xi_{2})}^{2} c_{3} + a_{3})} z_{13} + {({(w_{41} \xi_{1} + w_{42} \xi_{2})}^{2} c_{4} + a_{4})} z_{14} - 1)} + I_{1})} \mbox{wi}_{21} + {({(w_{21} \xi_{1} + w_{22} \xi_{2})} {({({(w_{11} \xi_{1} + w_{12} \xi_{2})}^{2} c_{1} + a_{1})} z_{21} + {({(w_{21} \xi_{1} + w_{22} \xi_{2})}^{2} c_{2} + a_{2})} z_{22} + {({(w_{31} \xi_{1} + w_{32} \xi_{2})}^{2} c_{3} + a_{3})} z_{23} + {({(w_{41} \xi_{1} + w_{42} \xi_{2})}^{2} c_{4} + a_{4})} z_{24} - 1)} + I_{2})} \mbox{wi}_{22} + {({(w_{31} \xi_{1} + w_{32} \xi_{2})} {({({(w_{11} \xi_{1} + w_{12} \xi_{2})}^{2} c_{1} + a_{1})} z_{31} + {({(w_{21} \xi_{1} + w_{22} \xi_{2})}^{2} c_{2} + a_{2})} z_{32} + {({(w_{31} \xi_{1} + w_{32} \xi_{2})}^{2} c_{3} + a_{3})} z_{33} + {({(w_{41} \xi_{1} + w_{42} \xi_{2})}^{2} c_{4} + a_{4})} z_{34} - 1)} + I_{3})} \mbox{wi}_{23} + {({(w_{41} \xi_{1} + w_{42} \xi_{2})} {({({(w_{11} \xi_{1} + w_{12} \xi_{2})}^{2} c_{1} + a_{1})} z_{41} + {({(w_{21} \xi_{1} + w_{22} \xi_{2})}^{2} c_{2} + a_{2})} z_{42} + {({(w_{31} \xi_{1} + w_{32} \xi_{2})}^{2} c_{3} + a_{3})} z_{43} + {({(w_{41} \xi_{1} + w_{42} \xi_{2})}^{2} c_{4} + a_{4})} z_{44} - 1)} + I_{4})} \mbox{wi}_{24} \end{array}$$\right)

W = matrix( [[1,0],[0,1],[1,0],[0,1]] )/sqrt(N/2)
Wi = matrix( [[1,0,1,0],[0,1,0,1]] )/sqrt(N/2)
Wi*W

\newcommand{\Bold}[1]{\mathbf{#1}}\left($$\begin{array}{rr} 1 & 0 \\ 0 & 1 \end{array}$$\right)

vxi2 = []
for i in range(0,M):
    temp = vxis[i].collect( [xi1,xi2][i] ).subs( c1=-4/N^2, c2=-4/N^2, c3=-4/N^2, c4=-4/N^2, a1=2/N, a2=2/N, a3=2/N, a4=2/N, z14=0, z32=0, z34=0, z13=0, z31=0, z33=0, z23=0, z41=0, z43=0, z24=0, z42=0, z44=0).expand()
    for j1 in range(0,N):
        for j2 in range(0,M):
            temp = eval('temp.subs( w' + str(j1+1) + str(j2+1) + ' = W[' + str(j1) + ',' + str(j2) + '], wi' + str(j2+1) + str(j1+1) + ' = Wi[' + str(j2) + ','+ str(j1)+'] )')
    vxi2.append(temp)
vector(vxi2).transpose()

\newcommand{\Bold}[1]{\mathbf{#1}}\left($$\begin{array}{r} -\frac{1}{16} \, \xi_{1}^{3} z_{11} - \frac{1}{16} \, \xi_{1} \xi_{2}^{2} z_{12} + \frac{1}{4} \, \xi_{1} z_{11} + \frac{1}{4} \, \xi_{1} z_{12} + \frac{1}{2} \, \sqrt{2} I_{1} + \frac{1}{2} \, \sqrt{2} I_{3} - \xi_{1} \\ -\frac{1}{16} \, \xi_{1}^{2} \xi_{2} z_{21} - \frac{1}{16} \, \xi_{2}^{3} z_{22} + \frac{1}{4} \, \xi_{2} z_{21} + \frac{1}{4} \, \xi_{2} z_{22} + \frac{1}{2} \, \sqrt{2} I_{2} + \frac{1}{2} \, \sqrt{2} I_{4} - \xi_{2} \end{array}$$\right)

$\dot{ \xi_1 } = ( 1-\sum_{i=1}^4 \xi_i^2) \xi_1 + \mu (\xi_2+\xi_1-\xi_1^3/3) \tau - c_1 \xi_1 \sum_i^4 \xi_i^2$

$\dot{ \xi_2 } = ( 1-\sum_{i=1}^4 \xi_i^2) \xi_2 + \mu ( - \xi_1 + a_1 + I_1 + b_1 \xi_2 ) / \tau - c_1  \xi_2 \sum_i^4 \xi_i^2$

var('alpha tau')
vxi3 = vector([vxi2[0].subs(I1=0,I3=0,z12=1,z11=8), vxi2[1].subs(I2=0,I4=2*alpha/sqrt(2), z22=0, z21=4)]).transpose(); vxi3

\newcommand{\Bold}[1]{\mathbf{#1}}\left($$\begin{array}{r} -\frac{1}{2} \, \xi_{1}^{3} - \frac{1}{16} \, \xi_{1} \xi_{2}^{2} + \frac{5}{4} \, \xi_{1} \\ -\frac{1}{4} \, \xi_{1}^{2} \xi_{2} + \alpha \end{array}$$\right)

solve( vxi3[0][0]==0, xi2 ), solve( vxi3[1][0]==0, xi2 )

\newcommand{\Bold}[1]{\mathbf{#1}}\left(\text{[ xi2 == -2*sqrt(-2*xi1^2 + 5), xi2 == 2*sqrt(-2*xi1^2 + 5) ]}, \text{[ xi2 == 4*alpha/xi1^2 ]}\right)

vxi3[0]

\newcommand{\Bold}[1]{\mathbf{#1}}\left(-\frac{1}{2} \, \xi_{1}^{3} - \frac{1}{16} \, \xi_{1} \xi_{2}^{2} + \frac{5}{4} \, \xi_{1}\right)

import numpy as np

vf1 = ((1-xi1^2-xi2^2)*xi1 + vxi3[0][0])._fast_float_(xi1,xi2)
vf2 = ((1-xi1^2-xi2^2)*xi2 + vxi3[1][0])._fast_float_(xi1,xi2,alpha)

test2 = ode_solver( function= lambda t,y: [vf1(y[0],y[1]),vf2(y[0],y[1],.3)], y_0=[.3,-.3] )
test2.ode_solve( t_span=(0,200), num_points=200 )

sol = []
t,y = test2.solution[0]
sol.append([])
for i in y:
    sol.append([])
for t,y in test2.solution:
    sol[0].append(t)
    for j in range(0,len(y)):
        sol[j+1].append(y[j])

import pylab as p
p.close()
p.plot(sol[0],sol[1],sol[0],sol[2])
p.savefig('model.png')

alpha = .5
sp = plot(-2*sqrt(-2*x^2 + 5),(-2,2)) + plot(2*sqrt(-2*x^2 + 5),(-2,2)) + plot(4*alpha/x^2,(-2,2))
sp.show(ymax=5,xmin=-1.6,xmax=1.6)
2 != 3

verbose 0 (3052: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 44 points. verbose 0 (3052: plot.py, generate_plot_points) Last error message: '' verbose 0 (3052: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 42 points. verbose 0 (3052: plot.py, generate_plot_points) Last error message: ''

\newcommand{\Bold}[1]{\mathbf{#1}}{\rm True}

Z = matrix([[8,0],[0,0]])
I0 = [0,0,0,0]

def testsys(t,y):

    if t>55 and t<65:
        boxcar = 1
    else:
        boxcar = 0
        
    dy1 = (1-y[0]**2-y[1]**2)*y[0] -1/16*( (y[1]**2-4)*Z[0,1] + (y[0]**2-4)*Z[0,0] + 16)*y[0]/2 + sqrt(2)*I0[0]/2 + sqrt(2)*I0[2]/2 - 4*boxcar
    dy2 = (1-y[0]**2-y[1]**2)*y[1] -1/16*( (y[1]**2-4)*Z[1,1] + (y[0]**2-4)*Z[1,0] + 16)*y[1]/9 + sqrt(2)*I0[1]/2 + sqrt(2)*I0[3]/2 + 2*boxcar
    
    return [dy1,dy2]

import pylab as p
test = ode_solver( function=testsys, y_0=[.3,.3] )
test.ode_solve( t_span=(0,100), num_points=2000 )

sol = []
t,y = test.solution[0]
sol.append([])
for i in y:
    sol.append([])
for t,y in test.solution:
    sol[0].append(t)
    for j in range(0,len(y)):
        sol[j+1].append(y[j])

import numpy as n
X,Y = n.meshgrid( n.linspace(-1.5,1.5,20), n.linspace(-1.5,1.5,20) )
U,V = n.zeros([20,20]), n.zeros([20,20])
for i in range(0,20):
    for j in range(0,20):
        vij = testsys(0,[X[i,j],Y[i,j]])
        U[i,j] = vij[0]
        V[i,j] = vij[1]
        
p.close()
p.quiver(X,Y,U,V)
p.plot(sol[1],sol[2])
p.savefig("pplane.png")

ode_solver?