Embedding Flows

151 days ago by mmwoodman

%hide #auto # this code should be run when the worksheet is loaded! def mkSymMat(el_name,m,n): A = [] for i in range(0,m): A.append([]) for j in range(0,n): A[i].append(var(el_name+str(i+1)+str(j+1))) return matrix(A) def mkSymVec(el_name,m): V = [] for i in range(0,m): V.append(var(el_name+str(i+1))) return(vector(V)) def rfold(op,ls): l = len(ls) if l>1: return op( rfold(op,ls[:(l-1)]) , ls[l-1] ) else: return ls[0] def lfold(op,ls): l = len(ls) if l>1: return op( ls[0], lfold(op,ls[1:l]) ) else: return ls[l-1] 
       

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.
%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.
%hide a1, b1, I1 = 0.2, 2, 0 a2, b2, I2 = 1,1,0 c1, c2 = 0, 1 tau, mu = 3, .2 def fsys(t,y): # sum of squares of state variables sumsq = 0.0 for i in range(0,4): sumsq += y[i]**2 # vector field components, 1st excitator dx1 = (1-sumsq)*y[0] + mu*( y[1]+y[0]-(y[0]**3)/3 )*tau - (t/1000)*y[0]*sumsq dx2 = (1-sumsq)*y[1] - mu*( y[0] -a1 -I1 -b1*y[1] )/tau - (t/1000)*y[1]*sumsq # 2nd excitator dx3 = (1-sumsq)*y[2] + mu*( y[3]+y[2]-(y[2]**3)/3 )*tau - (1-t/1000)*y[2]*sumsq dx4 = (1-sumsq)*y[3] - mu*( y[2] -a2 -I2 -b2*y[3] )/tau - (1-t/1000)*y[3]*sumsq return [dx1,dx2,dx3,dx4] S = ode_solver( function=fsys, y_0=[1,1,1,1] ) S.ode_solve( t_span=(0,1000), num_points=1000 ) sol = [] t,y = S.solution[0] sol.append([]) for i in y: sol.append([]) for t,y in S.solution: sol[0].append(t) for j in range(0,len(y)): sol[j+1].append(y[j]) t, c1, c2 = [], [],[] for i in range(1,1001): t.append(i*1.0) c1.append( i/1000.0 ) c2.append(1-i/1000.0) 
       
%hide import pylab as p p.close() p.subplot(211) p.plot(t,c1,label='$c_1$') p.plot(t,c2,label='$c_2$') p.ylabel("$c_i$") p.legend(loc=7) p.subplot(212) p.plot(sol[0],sol[1],label=r"$\xi_1$") p.plot(sol[0],sol[3],label=r"$\xi_3$") p.ylabel("activity") p.xlabel("time") p.legend(loc=8) p.savefig("fig1.png") 
       

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).

 
%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.
%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)
%hide # first define functional dot product wrt vec u, with limits a,b that # are scalars or vectors w/ length = length(u) def fdot(f,g,u,a,b): prod = f*g if isinstance(a,list) and isinstance(b,list): for i in range( 0, len(u) ): prod = integrate(prod,u[i],a[i],b[i]) else: for i in range( 0, len(u) ): prod = integrate(prod,u[i],a,b) return(prod) # then a norm def fnorm(f,u,a,b): return sqrt(fdot(f,f,u,a,b)) # now a theta def ftheta(f,g,u,a,b): return acos( fdot(f,g,u,a,b) / (fnorm(f,u,a,b)*fnorm(g,u,a,b)) ) # now, compute Theta over a region of (-10,10)^2 a,b = -2, 2 theta = map( lambda f,g: ftheta(f,g,xi,a,b), f_q, f_xi ) Theta = sum( map( lambda f: f**2, theta ) ) Theta 
        
\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(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 'Theta = %g' % Theta.subs(a1=tmin[0],a21=tmin[1],a22=tmin[2],a31=tmin[3],a32=tmin[4],z0=tmin[5]).numerical_approx() 
        
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).
%hide # create fast versions of the two vector fields fast_f_xi = map(lambda f:fast_float(f,'xi1','xi2'),f_xi) fast_f_q = map(lambda f:fast_float(f,'xi1','xi2'), map( lambda f: f.subs(a1=tmin[0],a21=tmin[1],a22=tmin[2],a31=tmin[3],a32=tmin[4],z0=tmin[5]), f_q ) ) # evaluate vec fields over and npts x npts grid npts = 40 n_Theta_args = [ 0.0, 0.0, 0.0 ] import numpy as np n_xi1, n_xi2 = np.linspace(a,b,npts), np.linspace(a,b,npts) vf_q_u, vf_q_v, vf_xi_u, vf_xi_v = np.zeros([npts,npts]), np.zeros([npts,npts]), np.zeros([npts,npts]), np.zeros([npts,npts]) for i in range(0,npts): for j in range(0,npts): vf_q_u[i,j] = fast_f_q[0]( n_xi1[i], n_xi2[j] ) vf_q_v[i,j] = fast_f_q[1]( n_xi1[i], n_xi2[j] ) vf_xi_u[i,j] = fast_f_xi[0]( n_xi1[i], n_xi2[j] ) vf_xi_v[i,j] = fast_f_xi[1]( n_xi1[i], n_xi2[j] ) for i in range(0,npts): for j in range(0,npts): n_Theta_args[0] += vf_q_u[i,j]*vf_xi_u[i,j] + vf_q_v[i,j]*vf_xi_v[i,j] n_Theta_args[1] += vf_q_u[i,j]**2 + vf_q_v[i,j]**2 n_Theta_args[2] += vf_xi_u[i,j]**2 + vf_xi_v[i,j]**2 n_Theta = acos( n_Theta_args[0] / (sqrt( n_Theta_args[1] )*sqrt( n_Theta_args[2] )) ) print 'Theta = %g' % n_Theta 
        
Theta = 1.83617e-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) 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) 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},{(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}\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 ) ) g_Theta 
        
\newcommand{\Bold}[1]{\mathbf{#1}}\arccos\left(-\frac{20}{5455583513727} \, \frac{{(291939899101 \, z_{11} + 4761066364129 \, z_{12} + 291939899101 \, z_{31} + 4761066364129 \, z_{32} + 291939899101 \, z_{41} + 4761066364129 \, z_{42})} \sqrt{\frac{122}{35}}}{\sqrt{\frac{2}{1959696563137922789941005} \, {(\left(3.29385193864 \times 10^{27}\right) \, z_{12} + \left(3.29385193864 \times 10^{27}\right) \, z_{32})} z_{42} + \frac{2}{1959696563137922789941005} \, {(\left(3.29385029392 \times 10^{27}\right) \, z_{11} + \left(3.29385029392 \times 10^{27}\right) \, z_{31})} z_{41} + \frac{1}{1959696563137922789941005} \, {(\left(3.29385029392 \times 10^{27}\right) \, z_{11} - 297055606876.0 \, z_{12})} z_{31} + \frac{1}{1959696563137922789941005} \, {(\left(3.29385029392 \times 10^{27}\right) \, z_{11} + 297055606876.0 \, z_{12})} z_{31} + 1680.79607623 \, z_{11}^{2} + 1680.7969155 \, z_{12}^{2} + 3361.59383101 \, z_{12} z_{32} + 1680.79607623 \, z_{31}^{2} + 1680.7969155 \, z_{32}^{2} + \frac{3293850293924868824397045184}{1959696563137922789941005} \, z_{41}^{2} + 1680.7969155 \, z_{42}^{2}}}\right)^{2} + \arccos\left(-\left(6.57909184807 \times 10^{-13}\right) \, \frac{{(23805322706581 \, z_{21} - 47610663641290 \, z_{22} + 23805322706581 \, z_{31} - 47610663641290 \, z_{32} + 23805322706581 \, z_{41} - 47610663641290 \, z_{42} + 23805322706581 \, z_{51} - 47610663641290 \, z_{52})}}{\sqrt{-{(-1680.7969155 \, z_{22} - 1680.7969155 \, z_{32})} z_{42} + {(1680.7969155 \, z_{22} + 1680.7969155 \, z_{32})} z_{42} - {(-1680.79607623 \, z_{21} - 1680.79607623 \, z_{31})} z_{41} - {(-1680.79607623 \, z_{21} - \left(1.51582450295 \times 10^{-13}\right) \, z_{22})} z_{31} + {(1680.79607623 \, z_{21} + 1680.79607623 \, z_{31})} z_{41} + {(1680.79607623 \, z_{21} - \left(1.51582450295 \times 10^{-13}\right) \, z_{22})} z_{31} - {(-1680.7969155 \, z_{22} - 1680.7969155 \, z_{32} - 1680.7969155 \, z_{42})} z_{52} + {(1680.7969155 \, z_{22} + 1680.7969155 \, z_{32} + 1680.7969155 \, z_{42})} z_{52} - {(-1680.79607623 \, z_{21} - \left(1.51582450295 \times 10^{-13}\right) \, z_{22} - 1680.79607623 \, z_{31} - \left(1.51582450295 \times 10^{-13}\right) \, z_{32} - 1680.79607623 \, z_{41} - \left(1.51582450295 \times 10^{-13}\right) \, z_{42})} z_{51} + {(1680.79607623 \, z_{21} - \left(1.51582450295 \times 10^{-13}\right) \, z_{22} + 1680.79607623 \, z_{31} - \left(1.51582450295 \times 10^{-13}\right) \, z_{32} + 1680.79607623 \, z_{41} - \left(1.51582450295 \times 10^{-13}\right) \, z_{42})} z_{51} + 1680.79607623 \, z_{21}^{2} + 1680.7969155 \, z_{22}^{2} + 3361.59383101 \, z_{22} z_{32} + 1680.79607623 \, z_{31}^{2} + 1680.7969155 \, z_{32}^{2} + 1680.79607623 \, z_{41}^{2} + 1680.7969155 \, z_{42}^{2} + 1680.79607623 \, z_{51}^{2} + 1680.7969155 \, z_{52}^{2}}}\right)^{2}

minimize using free variables
g_varlist = g_Theta.variables() g_varlist, len(g_varlist) 
        
\newcommand{\Bold}[1]{\mathbf{#1}}((z11, z12, z21, z22, z31, z32, z41, z42, z51, z52), 10)
g_ic = [] for i in range(0,len(g_varlist)): g_ic.append(random()) g_min = list( minimize(g_Theta,g_ic) ) 
        
Optimization terminated successfully.
         Current function value: 0.720901
         Iterations: 13
         Function evaluations: 14
         Gradient evaluations: 14

0.72 radians is not too good.