# Analysis of Wilhelm - "The Smallest Chemical Reaction with Bistability"

# Preliminaries

## Imports


```python
import tellurium as te
import sympy
import matplotlib.pyplot as plt
import numpy as np
from common_python.sympy import sympyUtil as su
from common_python.ODEModel.ODEModel import ODEModel
```

## Constants


```python
su.addSymbols("k_1 x k_2 k_3 y k_4")
```

# Equations


```python
stateDct = {
    x: 2 * k_1 * y - k_2 * x**2 - k_3 * x * y - k_4 * x,
    y: k_2 * x **2 - k_1 * y
}
```

# Analysis of Fixed Points


```python
stateDct[x]
```




$\displaystyle 2 k_{1} y - k_{2} x^{2} - k_{3} x y - k_{4} x$




```python
stateDct[y]
```




$\displaystyle - k_{1} y + k_{2} x^{2}$




```python
model = ODEModel(stateDct)
```


```python
valueDcts = [f.valueDct for f in model.fixedPoints]
valueDcts
```




    [{x: 0.0, y: 0.0},
     {x: -(-1.0*sqrt(k_1)*sqrt(k_2)*sqrt(0.25*k_1*k_2 - k_3*k_4)/k_3**2 + 0.5*(k_1*k_2 - 2.0*k_3*k_4)/k_3**2)*(-k_1*k_2 + k_3**2*(-1.0*sqrt(k_1)*sqrt(k_2)*sqrt(0.25*k_1*k_2 - k_3*k_4)/k_3**2 + 0.5*(k_1*k_2 - 2.0*k_3*k_4)/k_3**2) + k_3*k_4)/(k_2*k_4),
      y: -1.0*sqrt(k_1)*sqrt(k_2)*sqrt(0.25*k_1*k_2 - k_3*k_4)/k_3**2 + 0.5*(k_1*k_2 - 2.0*k_3*k_4)/k_3**2},
     {x: -(1.0*sqrt(k_1)*sqrt(k_2)*sqrt(0.25*k_1*k_2 - k_3*k_4)/k_3**2 + 0.5*(k_1*k_2 - 2.0*k_3*k_4)/k_3**2)*(-k_1*k_2 + k_3**2*(1.0*sqrt(k_1)*sqrt(k_2)*sqrt(0.25*k_1*k_2 - k_3*k_4)/k_3**2 + 0.5*(k_1*k_2 - 2.0*k_3*k_4)/k_3**2) + k_3*k_4)/(k_2*k_4),
      y: 1.0*sqrt(k_1)*sqrt(k_2)*sqrt(0.25*k_1*k_2 - k_3*k_4)/k_3**2 + 0.5*(k_1*k_2 - 2.0*k_3*k_4)/k_3**2}]




```python
sympy.simplify(valueDcts[1][x] - valueDcts[2][x])
```




$\displaystyle - \frac{2.0 \sqrt{k_{1}} \sqrt{0.25 k_{1} k_{2} - k_{3} k_{4}}}{\sqrt{k_{2}} k_{3}}$



So, require that $0.25 k_1 k_2 \geq k_3 k_4$ for real fixed point.

## Non-zero fixed points


```python
fp = valueDcts[1]
fp[x].expand()
```




$\displaystyle - \frac{1.0 \sqrt{k_{1}} \sqrt{0.25 k_{1} k_{2} - k_{3} k_{4}}}{\sqrt{k_{2}} k_{3}} + \frac{0.5 k_{1}}{k_{3}}$




```python
fp = valueDcts[2]
fp[x].expand()
```




$\displaystyle \frac{1.0 \sqrt{k_{1}} \sqrt{0.25 k_{1} k_{2} - k_{3} k_{4}}}{\sqrt{k_{2}} k_{3}} + \frac{0.5 k_{1}}{k_{3}}$




```python
fp = valueDcts[1]
fp[y].expand()
```




$\displaystyle - \frac{1.0 \sqrt{k_{1}} \sqrt{k_{2}} \sqrt{0.25 k_{1} k_{2} - k_{3} k_{4}}}{k_{3}^{2}} + \frac{0.5 k_{1} k_{2}}{k_{3}^{2}} - \frac{1.0 k_{4}}{k_{3}}$




```python
fp = valueDcts[2]
fp[y].expand()
```




$\displaystyle \frac{1.0 \sqrt{k_{1}} \sqrt{k_{2}} \sqrt{0.25 k_{1} k_{2} - k_{3} k_{4}}}{k_{3}^{2}} + \frac{0.5 k_{1} k_{2}}{k_{3}^{2}} - \frac{1.0 k_{4}}{k_{3}}$



## When are fixed points > 0?

### Value of $x$


```python
# Evaluation of x^* for fixed point 1
fp = valueDcts[1]
expandedX = fp[x].expand()
expandedX
```




$\displaystyle - \frac{1.0 \sqrt{k_{1}} \sqrt{0.25 k_{1} k_{2} - k_{3} k_{4}}}{\sqrt{k_{2}} k_{3}} + \frac{0.5 k_{1}}{k_{3}}$



This is equivalent to the following.


```python
inequality = (k_2 * k_3 * expandedX.args[0]) **2 / (k_1 * k_2) > (k_2 * k_3 * expandedX.args[1]) ** 2 / (k_1 * k_2)
inequality
```




$\displaystyle 0.25 k_{1} k_{2} > 0.25 k_{1} k_{2} - 1.0 k_{3} k_{4}$



Since all $k > 0 $, $x^{\star} > 0$.

### Value of $y$


```python
# Evaluation of y^* for fixed point 1. This should be greater than 0.
fp = valueDcts[1]
expandedY = fp[y].expand()
expandedY > 0
```




$\displaystyle - \frac{1.0 \sqrt{k_{1}} \sqrt{k_{2}} \sqrt{0.25 k_{1} k_{2} - k_{3} k_{4}}}{k_{3}^{2}} + \frac{0.5 k_{1} k_{2}}{k_{3}^{2}} - \frac{1.0 k_{4}}{k_{3}} > 0$




```python
su.addSymbols("k_5")  # k_5 = k_3 * k_4 /(k_1 * k_2)
term = k_3 **2 / (k_1 * k_2)
# lhs = term * expandedY.args[1]
rhs = term * expandedY.args[0] - term * expandedY.args[2]
rhs = rhs.subs(k_3 * k_4 /(k_1 * k_2), k_5)
term * expandedY.args[1] > - sympy.simplify(term * expandedY.args[0] - term * expandedY.args[2])
```




$\displaystyle 0.5 > \frac{1.0 k_{3} k_{4}}{k_{1} k_{2}} - \frac{1.0 \sqrt{0.25 k_{1} k_{2} - k_{3} k_{4}}}{\sqrt{k_{1}} \sqrt{k_{2}}}$




```python
# Re-written
0.5 > k_5 - sympy.sqrt(0.25  - k_5)
```




$\displaystyle k_{5} - \sqrt{0.25 - k_{5}} < 0.5$




```python
-sympy.sqrt(0.25  - k_5) < 0.5 - k_5
```




$\displaystyle - \sqrt{0.25 - k_{5}} < 0.5 - k_{5}$




```python
k_5 - 0.5 < sympy.sqrt(0.25  - k_5)
```




$\displaystyle k_{5} - 0.5 < \sqrt{0.25 - k_{5}}$



This is true for $k_5 \in [0, 0.25]$. If $k_5 > 0.25$, there is no real solution.

### Evaluation of fixed points

Plot for values of $k_5$ outside the operating region.

ec6857d344b9d3599a2b17bc78189aa6ab4d4b73
ipynb
Jupyter Notebook
notebooks/Wilhelm.ipynb
sys-bio/IntegralControl
45954551a7bf10b111b1f1b15f53a435839517c5
MIT
notebooks/Wilhelm.ipynb
sys-bio/IntegralControl
45954551a7bf10b111b1f1b15f53a435839517c5
MIT
notebooks/Wilhelm.ipynb
sys-bio/IntegralControl
45954551a7bf10b111b1f1b15f53a435839517c5
MIT
Qwen/Qwen-72B
1. YES
2. YES