# Species eradications and introductions

## A scaling of the LV food web

The Lotka-Volterra food web has dynamics

\begin{equation}
\frac{dn_{t;i}}{dt} = n_{t;i}
\left(r_i + \sum_j a_{ij} n_{t;j}\right)
\end{equation}

with a steady state

\begin{equation}
\mathbf{n} = A^{-1} (-\mathbf{r})
\end{equation}

Let $N$ be the diagonal matrix formed from $\mathbf{n}$. Then the Jacobian

\begin{equation}
J = \left[ \left. \frac{\partial n_{t;i}'}{\partial
n_{t;j}}\right|_{\mathbf{n}}\right] = NA
\end{equation}

Interested in scaling a system such that $\mathbf{\hat{n}} = 1$, so choose

\begin{equation}
\hat{n}_{t;i} = \frac{n_{t;i}}{n_i}
\end{equation}

so the dynamics

\begin{equation}
\hat{n}_{t;i}' 
= \hat{n}_{t;i} \left(
    r_i + \sum_j a_{ij} \hat{n}_{t;i} n_j
\right) 
\end{equation}

and the Jacobian

\begin{equation}
M = \left[ \left.\frac{\partial \hat{n}_{t;i}'}{\partial
\hat{n}_{t;j}}\right|_{\mathbf{\hat{n}} = [1]} \right] = 
AN = N^{-1}JN
\end{equation}

At the steady state

\begin{equation}
r_{i} = -\sum_{j} a_{ij} n_j = -\sum_j m_{ij}
\end{equation}

The sensitivity matrix of the scaled system is

\begin{equation}
M^{-1} = N^{-1} J^{-1} N 
\end{equation}

## Special properties of the scaled system

1. $M$ is a similarity transform of $J$, therefore the eigenvalues (but not eigenvectors) are equal.

2. $\mathbf{n}$ is strictly positive and $N$ is a diagonal matrix, therefore the signs (but not the magnitudes) of the elements of $J^{-1}$ are preserved.

$M$ can be further scaled so that its elements $m_{ij} = (0,1)$, scaling the eigenvalues and the elements of the sensitivity matrix.

$M$ therefore represents a suite of systems conforming the the transformations above.

### Example of scaling and special properties

Consider an LV food web with the following structure


```
import networkx as nx
import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline

# Define positive and negative interactions
positive_edges_list = [(4,2), (4,3), (2,3), (2,1), (3,1), (3,0)]
negative_edges_list = [(0,0), (1,1), (2,2), (3,3), (4,4),
                       (2,4), (3,4), (3,2), (1,2), (1,3), (0,3)]

# Create digraph
G = nx.DiGraph()
G.add_edges_from(negative_edges_list, color='red', weight=-1)
G.add_edges_from(positive_edges_list, color='green', weight=1)
node_labels = {node:str(node) for node in G.nodes()}

# Draw graph
pos=nx.spring_layout(G)
nx.draw_networkx_nodes(G, pos)
nx.draw_networkx_edges(G, pos, 
                        edgelist=positive_edges_list, edge_color='g')
nx.draw_networkx_edges(G, pos, 
                        edgelist=negative_edges_list, edge_color='r')
nx.draw_networkx_labels(G,pos,node_labels)
```

We can create a scaled Jacobian $M$


```
sz = G.order()

# Interaction-sign matrix
Mq = np.zeros(shape=(sz,sz))
Mq[tuple(zip(*positive_edges_list))] = 1
# i.e. Mq[((4, 4, 2, 2, 3, 3), (2, 3, 3, 1, 1, 0))]
Mq[tuple(zip(*negative_edges_list))] = -1

# Assign random uniform values to M
M = np.multiply(np.random.random_sample((sz, sz)), Mq)
eigM = np.linalg.eigvals(M)

print 'M = \n', M, '\n'
print 'eigenvalues of M: \n', eigM
```

    M = 
    [[-0.58152272  0.          0.         -0.14127391  0.        ]
     [ 0.         -0.80860701 -0.57168714 -0.59986014  0.        ]
     [ 0.          0.5646323  -0.55274168  0.50760127 -0.59904304]
     [ 0.16858447  0.95959136 -0.85644885 -0.06868719 -0.24935855]
     [ 0.          0.          0.03009555  0.53615731 -0.07643393]] 
    
    eigenvalues of M: 
    [-0.49151243+1.23541285j -0.49151243-1.23541285j  0.18055314+0.j
     -0.61868079+0.j         -0.66684001+0.j        ]


This corresponds to a whole family of possible systems. For example, we could define them by their steady state values:


```
n1 = np.random.random_sample(sz)
n2 = np.random.random_sample(sz)
print n1, '\n'
print n2, '\n'
```

    [ 0.49459123  0.42365168  0.70069605  0.5194794   0.60422636] 
    
    [ 0.83482282  0.83135622  0.88581108  0.33283174  0.1965191 ] 
    



```
N1 = np.diag(n1)
N2 = np.diag(n2)
A1 = np.dot(M, np.linalg.inv(N1))
A2 = np.dot(M, np.linalg.inv(N2))

print A1, '\n'
print A2, '\n'
```

    [[-1.17576432  0.          0.         -0.27195286  0.        ]
     [ 0.         -1.90865998 -0.81588463 -1.15473326  0.        ]
     [ 0.          1.33277485 -0.78884658  0.97713455 -0.99142157]
     [ 0.34085617  2.26504792 -1.22228296 -0.13222313 -0.41269063]
     [ 0.          0.          0.04295093  1.03210506 -0.12649883]] 
    
    [[-0.6965822   0.          0.         -0.42446044  0.        ]
     [ 0.         -0.97263602 -0.64538269 -1.80229244  0.        ]
     [ 0.          0.67917012 -0.623995    1.52509872 -3.04826878]
     [ 0.20194042  1.15424812 -0.96685272 -0.20637211 -1.26887693]
     [ 0.          0.          0.03397513  1.61089596 -0.38893892]] 
    


They have equal eigenvalues at the steady state


```
J1 = np.dot(N1,A1)
J2 = np.dot(N2,A2)

eigM = np.linalg.eigvals(M)
eigJ1 = np.linalg.eigvals(J1)
eigJ2 = np.linalg.eigvals(J2)

print eigM, '\n'
print eigJ1, '\n'
print eigJ2, '\n'
```

    [-0.49151243+1.23541285j -0.49151243-1.23541285j  0.18055314+0.j
     -0.61868079+0.j         -0.66684001+0.j        ] 
    
    [-0.49151243+1.23541285j -0.49151243-1.23541285j  0.18055314+0.j
     -0.66684001+0.j         -0.61868079+0.j        ] 
    
    [-0.49151243+1.23541285j -0.49151243-1.23541285j -0.61868079+0.j
     -0.66684001+0.j          0.18055314+0.j        ] 
    


The sign structure of the sensitivity matrix is also shared


```
print np.linalg.inv(-M), '\n'
print np.linalg.inv(-J1), '\n'
print np.linalg.inv(-J2), '\n'
```

    [[ 1.74245556  0.01785382 -0.10828187  0.0787587   0.59170512]
     [-0.16038225  0.69565713  0.16539543 -0.55322963  0.50859035]
     [ 0.32546384  0.84236854 -0.70162248  1.12266937  1.83629001]
     [-0.09398404 -0.07349128  0.4457183  -0.3241927  -2.43562293]
     [-0.53111518 -0.1838365   2.85029731 -1.83205217 -3.27881743]] 
    
    [[ 1.74245556  0.0208434  -0.07643152  0.07498538  0.48434193]
     [-0.13737852  0.69565713  0.10000064 -0.45117605  0.35659675]
     [ 0.46109033  1.39323018 -0.70162248  1.5143045   2.12946878]
     [-0.09871338 -0.09011461  0.33044496 -0.3241927  -2.09400983]
     [-0.6488465  -0.26219384  2.45787708 -2.13092995 -3.27881743]] 
    
    [[ 1.74245556  0.01792826 -0.10204905  0.19754593  2.51359249]
     [-0.15971627  0.69565713  0.15522782 -1.38187209  2.1515453 ]
     [ 0.34534212  0.89754471 -0.70162248  2.987915    8.27708869]
     [-0.03747007 -0.02942208  0.16747273 -0.3241927  -4.12505766]
     [-0.12502566 -0.04345596  0.63234462 -1.08172749 -3.27881743]] 
    


The systems in the family also share their $\mathbf{r}$, so whether this family of systems make sense or not depends upon what kinds of assumptions we're willing to make about the intrinsic growth rates. Perhaps for some top predators we might like to require that $r_i \leq 0$


```
r1 = -np.dot(A1,n1)
r2 = -np.dot(A2,n2)

print r1
print r2
```

    [ 0.72279663  1.98015428  0.07955115  0.04631876 -0.48981893]
    [ 0.72279663  1.98015428  0.07955115  0.04631876 -0.48981893]


Here is another scaled system, this time assuming that some scaling factor has been applied to the $M$ (for example, to keep its values between zero and 1). The sign-structure of the eigenvalues, the sensitivity matrix, and the intrinsic growth rates $\mathbf{r}$ is preserved. Given that our concerns are usually about the signs of these attributes, we can say that our $M$s search the entire space of possible outcomes in this regard.


```
alpha = 42
Mf = alpha*M # Just apply some scaling
n3 = np.random.random_sample(sz)
N3 = np.diag(n3)
A3 = np.dot(Mf, np.linalg.inv(N3))
r3 = -np.dot(A3,n3)
J3 = np.dot(N3,A3)

print 'eiganvalues: \n', eigJ1
print np.linalg.eigvals(J3)/alpha, '\n'
print 'r vector: \n', r1
print r3/alpha
```

    eiganvalues: 
    [-0.49151243+1.23541285j -0.49151243-1.23541285j  0.18055314+0.j
     -0.66684001+0.j         -0.61868079+0.j        ]
    [-0.49151243+1.23541285j -0.49151243-1.23541285j  0.18055314+0.j
     -0.61868079+0.j         -0.66684001+0.j        ] 
    
    r vector: 
    [ 0.72279663  1.98015428  0.07955115  0.04631876 -0.48981893]
    [ 0.72279663  1.98015428  0.07955115  0.04631876 -0.48981893]


## Species deletion

### Effect on steady state

The steady state of the full system $-\mathbf{r} = A \mathbf{n}$
may also be expressed as

\begin{equation}
-\mathbf{r}_{i\neq k} = A_{i \neq k,j \neq k}
\mathbf{n}_{i\neq k} +
A_{i \neq k,k} n_k
\end{equation}

where e.g. $A_{i \neq k,j \neq k}$ denotes the original $A$ matrix with its $k^{\text{th}}$ rows and columns removed.

A vector of the steady-state population size for all but the $k^{\text{th}}$ species can then be found

\begin{equation}
\mathbf{n}_{i\neq k} = (A_{i \neq k,j\neq k})^{-1}
\left( 
    - \mathbf{r}_{i\neq k} -
    n_k A_{i \neq k, k}
\right)
\end{equation}

We are interested in removal of the $k^{\text{th}}$ species. The steady state of the subsystem formed by removing species
$k$ is

\begin{equation}
\mathbf{s} =
(A_{i \neq k,j \neq k})^{-1}
\left(
-\mathbf{r}_{i\neq k} \right)
\end{equation}

We wish to know
if the removal of species $k$ caused any other species to go
extinct. Species extinction occurs when

\begin{equation}
\begin{split}
s_i &< 0 \\
\delta_i = \frac{s_i - n_i}{n_i} &< -1 
\end{split}
\end{equation}

Taking the difference between the original full system and the system that results after removal of species $k$

\begin{equation}
\begin{split}
\mathbf{s} - \mathbf{n}_{i \neq k}
& = A^{-1}_{i \neq k,j \neq k} \left(
    n_k A_{i\neq k, k} \right) \\
& = A^{-1}_{i \neq k,j \neq k} 
    M_{i\neq k,k}
\end{split}
\end{equation}

Because $A^{-1} = N M^{-1}$ and because $N$ is a diagonal matrix then $A^{-1}_{i\neq k, j \neq k} = 
N_{i\neq k, j \neq k}(M_{i\neq k, j \neq k})^{-1}$. Therefore

\begin{equation}
\mathbf{s} - \mathbf{n}_{i \neq k}
= N_{i\neq k, j \neq k}(M_{i\neq k, j \neq k})^{-1}
    M_{i\neq k,k}
\end{equation}

and so the effect of the removal of species $k$ on the steady state of the remaining species in the system can be found for all systems conforming to $M$ from the elements of $M$ alone

\begin{equation}
\delta_i
= \left[ (M_{i\neq k, j \neq k})^{-1}
    M_{i\neq k,k} \right]_i
\end{equation}


#### Example of effect on steady state

First we'll verify the re-expression of the steady state


```
k = 2 # Our species to remove

# Removal of kth row and column from system 1 from above
A_notk = np.delete(A1, (k), axis=0)
A_notk_k = A_notk[:,k]
A_notk = np.delete(A_notk, (k), axis=1)
r_notk = np.delete(r1, (k), axis=0)

print 'A = \n', A1, '\n'
print 'A_notk = \n', A_notk, '\n'

# Recalculate n
n_notk = np.linalg.solve(A_notk,-r_notk - n1[k]*A_notk_k)

print 'n = \n', n1, '\n'
print 'n_notk = \n', n_notk, '\n'
```

    A = 
    [[-1.17576432  0.          0.         -0.27195286  0.        ]
     [ 0.         -1.90865998 -0.81588463 -1.15473326  0.        ]
     [ 0.          1.33277485 -0.78884658  0.97713455 -0.99142157]
     [ 0.34085617  2.26504792 -1.22228296 -0.13222313 -0.41269063]
     [ 0.          0.          0.04295093  1.03210506 -0.12649883]] 
    
    A_notk = 
    [[-1.17576432  0.         -0.27195286  0.        ]
     [ 0.         -1.90865998 -1.15473326  0.        ]
     [ 0.34085617  2.26504792 -0.13222313 -0.41269063]
     [ 0.          0.          1.03210506 -0.12649883]] 
    
    n = 
    [ 0.49459123  0.42365168  0.70069605  0.5194794   0.60422636] 
    
    n_notk = 
    [ 0.49459123  0.42365168  0.5194794   0.60422636] 
    


We can find the new steady state after removal of species $k$.


```
s = np.linalg.solve(A_notk, -r_notk)

print 's = \n', s, '\n'
```

    s = 
    [ 0.41826063  0.52352028  0.84948803  3.05885802] 
    


The scaled changes in the steady-state $\delta$ can be found directly


```
delta = (s - n_notk)/n_notk

print 'delta = \n', delta, '\n'
```

    delta = 
    [-0.15433067  0.2357328   0.63526799  4.06243725] 
    


and we see here that the species that have gone extinct have a $\delta_i < -1$.

These scaled changes in the steady-state population size can be found directly from our scaled $M$


```
Mf_notk = np.delete(Mf, (k), axis=0)
Mf_notk_k = Mf_notk[:,k]
Mf_notk = np.delete(Mf_notk, (k), axis=1)

delta_from_M = np.linalg.solve(Mf_notk, Mf_notk_k)

print 'delta calculated from M = \n', delta_from_M, '\n'
```

    delta calculated from M = 
    [-0.15433067  0.2357328   0.63526799  4.06243725] 
    


$M$ can also be used to find the effect of species deletion for the whole family of systems. The scaled changes $\delta$ are the same for the second system above that scaled to $M$. 


```
# Removal of kth row and column from system 2 above, to compare
A_notk = np.delete(A2, (k), axis=0)
A_notk = np.delete(A_notk, (k), axis=1)
s = np.linalg.solve(A_notk, -r_notk)
n_notk = np.delete(n2, (k), axis=0)
delta2 = (s - n_notk)/n_notk

print 'delta from another system corresponding to M = \n', delta2, '\n'

```

    delta from another system corresponding to M = 
    [-0.15433067  0.2357328   0.63526799  4.06243725] 
    


### Effect on stability

The scaled changes relate $\mathbf{s}$ and $\mathbf{n}$ by 

\begin{equation}
s_i = n_i (\delta_i + 1)
\end{equation}

In the original system

\begin{equation}
J_{ij} = n_i a_{ij}
\end{equation}

and in the subsystem

\begin{equation}
\begin{split}
J_{s,ij} &= s_i a_{ij} \\
&= J_{ij} (\delta_i + 1) 
\end{split}
\end{equation}
so
\begin{equation}
J_s = \text{diag}(\boldsymbol{\delta}+1) J
\end{equation}

In the scaled system
\begin{equation}
m_{ij} = a_{ij} n_j
\end{equation}

and in the scaled subsystem

\begin{equation}
b_{ij} = m_{ij} (\delta_i + 1) 
\end{equation}

so 

\begin{equation}
B = M \text{diag}(\boldsymbol{\delta}+1) 
\end{equation}

So
\begin{equation}
\begin{split}
B &= N^{-1} J N \text{diag}(\boldsymbol{\delta}+1) \\
 &= N^{-1} (\text{diag}(\boldsymbol{\delta}+1))^{-1} J_s
 N \text{diag}(\boldsymbol{\delta}+1) \\
 &= N^{-1} (\text{diag}(\boldsymbol{\delta}+1))^{-1} J_s
 \text{diag}(\boldsymbol{\delta}+1) N 
\end{split}
\end{equation}
which is a similarity transform and therefore
\begin{equation}
\text{eig}(B) = \text{eig}(J_s)
\end{equation}

Also if $M$ is further scaled by $\alpha$, then $B$ and the eigenvalues are scaled by $\alpha$.

#### Example of effect on stability


```
J_s = np.dot(np.diag(s), A_notk)
print 'eig(Js) = \n', np.linalg.eigvals(J_s), '\n'

B = np.dot(Mf_notk, np.diag(delta2+1))

print 'eig(B) = \n', np.linalg.eigvals(B)/alpha, '\n'
```

    eig(Js) = 
    [-0.39291970+1.46001967j -0.39291970-1.46001967j -0.71629101+0.j
     -0.48813165+0.j        ] 
    
    eig(B) = 
    [-0.39291970+1.46001967j -0.39291970-1.46001967j -0.71629101+0.j
     -0.48813165+0.j        ] 
    


# Species addition

## Effect on steady state

We wish to know if the addition of species $k$ caused any other species to go extinct. Extinction of a resident species occurs when
\begin{equation}
\varepsilon_i = \frac{n_{i \neq k} - s_i}{s_i} < -1 
\end{equation}
As before,
\begin{equation}
\mathbf{n}_{i\neq k} = (A_{i \neq k,j\neq k})^{-1}
\left( - \mathbf{r}_{i\neq k} - n_k A_{i \neq k, k}
\right)
\end{equation}
\begin{equation}
\mathbf{s} =
(A_{i \neq k,j \neq k})^{-1}
\left(
-\mathbf{r}_{i\neq k} \right)
\end{equation}
so
\begin{equation}
\begin{split}
\mathbf{n}_{i \neq k} - \mathbf{s} 
 & = -(A_{i \neq k,j\neq k})^{-1} (n_k A_{i \neq k, k}) \\
 & = -(A_{i \neq k,j\neq k})^{-1} M_{i\neq k,k} \\
\end{split}
\end{equation}

Let the Jacobian of the subsystem be *** FIX THIS
$B = A_{i \neq k,j\neq k} S$ where $S$ is the diagonal matrix formed of
$\mathbf{s}$. Then, because $S$ is diagonal,
$(A_{i \neq k,j\neq k})^{-1} = S B^{-1}$. So

\begin{equation}
\mathbf{n}_{i \neq k} -
\mathbf{s} 
= SB^{-1} M_{i\neq k,k}
\end{equation}
so
\begin{equation}
\varepsilon_i = [- B^{-1} M_{i\neq k,k}]_i
\end{equation}

## Invasion criteria for species $k$

At the new steady state
\begin{equation}
r_k  = -a_{kk} n_k - \sum_{j \neq k} a_{kj} n_j 
\end{equation}

The new species introduced at low numbers can invade if
\begin{equation}
\begin{split}
\gamma_k 
 & = r_k + \sum_{i \neq k} a_{ki} s_i > 0 \\
\gamma_k
 & = -m_{kk} - \sum_{j \neq k} m_{kj}
     \left( \frac{\varepsilon_j}{\varepsilon_j + 1} \right) > 0\\
\end{split}
\end{equation}


## Filling in the remainder of the new $M$

The old scaled Jacobian

\begin{equation}
b_{ij} = a_{ij} s_j
\end{equation}

At the new steady state, the scaled Jacobian

\begin{equation}
\begin{split}
m_{ij} 
 & = a_{ij} n_j \\
 & = a_{ij} s_j(\varepsilon_j + 1) \\
 & = b_{ij} (\varepsilon_j + 1) \\
\end{split}
\end{equation}
so
\begin{equation}
M_{i \neq k, j \neq k} = B \; \text{diag}(\varepsilon + 1)
\end{equation}

## Example of successful species addition

We'll take one of the scaled Jacobians above as our starting subystem


```
B = M # Turn previous system into our new subsystem
```

We'll add one species $k$ to the final row and column, such that it feeds upon species 2 and 4, and is fed upon by species 3.


```
k = sz
Mq_notk_k = np.array([0, 0, +1, -1, +1])
```

First, we'll find values for the new final column of $M$ that will not lead to the extinction of any of the resident species


```
# Find an invasion that will not cause the extinction of 
# resident species

mineps = -2

while mineps <= -1:

    M_notk_k = np.random.random_sample(sz) * Mq_notk_k
    eps = np.dot(-np.linalg.inv(B), M_notk_k)
    mineps = min(eps)
    
print 'eps = \n', eps
```

    eps = 
    [-0.03633213  0.40438499 -0.72889433  0.14955316  1.18732908]


The new species' ability to invade the system is determined by the final row of the new $M$. We search for values that will satisfy the invasion criteria $\gamma_k > 0$


```
# Find an invasion that is successful

Mq_k_notk = -Mq_notk_k

gamma = -1

while gamma <= 0:
    m_kk = - np.random.random_sample()
    M_k_notk = np.random.random_sample(sz) * Mq_k_notk
    
    gamma = - m_kk - sum(M_k_notk * (eps/(1+eps)))

print 'gamma = ', gamma

```

    gamma =  0.842156209921


We can now construct the complete $M$, with which we may wish to test for stability or further invasions / removals.


```
# Construct the new M
M = np.dot(B, np.diag(eps+1))

M = np.concatenate((M, np.array([M_notk_k]).T), axis=1)
M_k_all = np.concatenate( (np.array([M_k_notk]), np.array([[m_kk]])), axis=1)
M = np.concatenate( (M, M_k_all), axis=0)

print 'M = \n', M, '\n'
print 'Stability test: eig(M) = \n', np.linalg.eigvals(M), '\n'
```

    M = 
    [[-0.56039476  0.          0.         -0.16240187  0.          0.        ]
     [ 0.         -1.13559554 -0.15498763 -0.68957112  0.          0.        ]
     [ 0.          0.79296113 -0.1498514   0.58351464 -1.31030426  0.00412874]
     [ 0.16245944  1.34763569 -0.23218814 -0.07895958 -0.54542922 -0.69983696]
     [ 0.          0.          0.00815907  0.61634134 -0.16718615  0.03250467]
     [ 0.          0.         -0.05252907  0.85987774 -0.72924792 -0.69940172]] 
    
    Stability test: eig(M) = 
    [-0.52642204+1.43510357j -0.52642204-1.43510357j  0.18233554+0.j
     -0.49341171+0.j         -0.56296611+0.j         -0.86450279+0.j        ] 
    


To show that this is consistent with the original system, we can choose one of the original unscaled systems in the same way as in the previous section. First we can look at $r$ (the value of $r_k$ might be a constraint on our invader)


```
# Find an original system to compare it to

r1 = -np.sum(B, axis=1)
r2 = -np.sum(M, axis=1)

print 'r from B \n', r1, '\n'
print 'r from M \n', r2, '\n'
```

    r from B 
    [ 0.72279663  1.98015428  0.07955115  0.04631876 -0.48981893] 
    
    r from M 
    [ 0.72279663  1.98015428  0.07955115  0.04631876 -0.48981893  0.62130097] 
    


Here is one of the possible unscaled systems, derived from the subsystem, a $s$ and $A$ pairing


```
s = np.random.random_sample(sz)
S = np.diag(s)

A_notk = np.dot(B, np.linalg.inv(S))

print 's = ', s, '\n'
print 'Old A = \n', A_notk
```

    s =  [ 0.77594037  0.80386546  0.94726495  0.54612625  0.10132065] 
    
    Old A = 
    [[-0.74944254  0.          0.         -0.25868361  0.        ]
     [ 0.         -1.00589844 -0.60351345 -1.09839097  0.        ]
     [ 0.          0.70239653 -0.58351328  0.92945774 -5.91234908]
     [ 0.21726472  1.19372135 -0.90412808 -0.12577163 -2.46108328]
     [ 0.          0.          0.03177099  0.9817461  -0.75437661]]


We'll use this information and the change in the steady state to get our new $n$ and from that the extended $A$


```
# Choose a_kk such that n_k positive, find new A
n_k = np.random.random_sample()
n = np.append(s*(eps+1), n_k)
N = np.diag(n)
A = np.dot(M, np.linalg.inv(N))

print 'New A = \n', A, '\n'
print 'Stability test: eig(J) = \n', np.linalg.eigvals(np.dot(N,A)), 
```

    New A = 
    [[-0.74944254  0.          0.         -0.25868361  0.          0.        ]
     [ 0.         -1.00589844 -0.60351345 -1.09839097  0.          0.        ]
     [ 0.          0.70239653 -0.58351328  0.92945774 -5.91234908  0.04953504]
     [ 0.21726472  1.19372135 -0.90412808 -0.12577163 -2.46108328 -8.39636731]
     [ 0.          0.          0.03177099  0.9817461  -0.75437661  0.38997823]
     [ 0.          0.         -0.20454535  1.36966575 -3.29050922 -8.39114545]] 
    
    Stability test: eig(J) = 
    [-0.52642204+1.43510357j -0.52642204-1.43510357j  0.18233554+0.j
     -0.49341171+0.j         -0.56296611+0.j         -0.86450279+0.j        ]



```

```

44b5c8e9701b7f387c5cab14e75c6a268d8b16dd
ipynb
Jupyter Notebook
docs/species_additions.ipynb
nadiahpk/qualitative-modelling
e80a2655d3067232ac5b65a6d618f79aa7846c3f
Unlicense
2021-09-08T17:38:28.000Z
2021-09-23T02:30:52.000Z
docs/species_additions.ipynb
nadiahpk/qualitative-modelling
e80a2655d3067232ac5b65a6d618f79aa7846c3f
Unlicense
docs/species_additions.ipynb
nadiahpk/qualitative-modelling
e80a2655d3067232ac5b65a6d618f79aa7846c3f
Unlicense
Qwen/Qwen-72B
1. YES
2. YES