### Analysis of deterministc approximation to the mitochondrial graph

Functionality of this notebook is similar to the accompanying python [script](./src/odes.py).  
Please feel free to adjust any of them to your needs.

Time-dependent behaviour of the mitochondria network can be represented as a system of differential-algebraic equations for the number of nodes $x_i$ of degrees $i$ = 1 to 3:  
\begin{array} 
\dot{x_1} = -a_2x_1x_2 + (3/2)bx_3 - a_1x_1(x_1-1) + 2bx_2 \\
\dot{x_2} = a_2x_1x_2 - (3/2)bx_3 \\
x_1 + 2x_2 +3x_3 = 2h
\end{array}
where $a_1$, $a_2$ and $b$ are ratec onstants for bulk fusion, side fusion and fission respectively. 


```python
# uncomment to plot into separate windows
#%matplotlib   

import math

import numpy as np
import scipy.integrate

from src.odes import \
    eqs, \
    is_stable, \
    node_numbers_equil, \
    num_deg2, \
    plot_node_numbers, \
    plot_phase_equil, \
    plot_stability, \
    plot_time_evol
```


```python
# Initialize the parameters:

def a(m):
     """Initialize fusion rate constants.
     """
     init = -np.floor(2 * m / 3)
     return np.logspace(init, init + m - 1, num=m, base=2)

m1, m2 = 57, 57          # grid dimensions
b = 1                    # fission rate constant
a1, a2 = a(m1), a(m2)    # fusion rate constant
c1, c2 = a1/b, a2/b      # reduced rates

h = [10000, 30000]       # total number of edges in the graph
```

#### The Steady state

Find the equilibrium solution as $t\rightarrow+\infty$  
using reduced parameters $c_1 = a_1/b$, $c_2 = a_2/b$:
\begin{align}
& 0 = c_1c_2x_1^3 + c_1(1 - c_2)x_1^2 + (1 - c_1)x_1 - 2h \\
& x_2 = c_1x_1(x_1 - 1)/2 \\
& x_3 = 2c_2x_1x_2/3
\end{align}

and plot it.


```python
x = [[[node_numbers_equil(cc1, cc2, hh) 
       for cc1 in c1] 
      for cc2 in c2] 
     for hh in h]
```


```python
for xx, hh in zip(x, h):
    plot_node_numbers(c1, c2, hh, xx, figsize=(15, 5))
```

Plot the solution in phase coordinates:


```python
plot_phase_equil(x, c1, h, figsize=(10,10))
```

Examine the solution.
The equilibrium is asymptotically stable if real parts of all eigenvalues of the Jacobian are strictly negative.


```python
st = [is_stable(xx, b, a1, a2, hh) for xx, hh in zip(x, h)]
```

Plot the stability map indicating stable and unstable solutions with blue and red markers respectively:


```python
for s, hh in zip(st, h):
    plot_stability(b, c1, c2, hh, s)
```

#### Ttransient bechavior

Slove the ODEs directly for specific parameters and plot the results:


```python
ht = h[0]                                   # graph total mass (in edges)
bt = b                                      # fission rate constant
a1t = a1[20]                                # end-to-end fusion rate constant
a2t = a2[30]                                # end-to-side fusion rate constant
tspan = [0., 20.]                           # time interval
tsol = np.linspace(tspan[0], tspan[1], 100) # time points for plotting
# initial values:
x1_0 = np.linspace(0, ht, 10)
x3_0 = (ht - x1_0) / 2

x123 = []
for x10, x30 in zip(x1_0, x3_0):
    # new scipy ivp solver: requires scipy >= 1.4:
    sol = scipy.integrate.solve_ivp(
        eqs, 
        t_span=tspan, 
        y0=[x10, x30], 
        args=(b, a1t, a2t, ht), 
        dense_output=True
    )
    x13 = sol.sol(tsol)
    x123.append([x13[0,:], num_deg2(x13, ht), x13[1,:]])
```


```python
plot_time_evol(b, a1t, a2t, ht, x123, tsol, figsize=(16, 5))
```


```python

```

492f96d8b49f25e780b419e7009b3af2640f5a86
ipynb
Jupyter Notebook
odes.ipynb
vsukhor/mitoSim-postProcessor
a0d42a0b0a69cb8f8e48b0ff8baf1802c3c432f4
MIT
odes.ipynb
vsukhor/mitoSim-postProcessor
a0d42a0b0a69cb8f8e48b0ff8baf1802c3c432f4
MIT
odes.ipynb
vsukhor/mitoSim-postProcessor
a0d42a0b0a69cb8f8e48b0ff8baf1802c3c432f4
MIT
Qwen/Qwen-72B
1. YES
2. YES