# Simple Biochemical Reaction Networks

This tutorial has been adapted from the tutorials in https://github.com/mic-pan/BGT_BiochemicalNetworkTutorials


```julia
# Since BondGraphs is not yet in the package manager, we will need to include it directly from Github
# NOTE: You will need Julia >= 1.7
using Pkg; Pkg.add(url="https://github.com/jedforrest/BondGraphs.jl")
using BondGraphs
```


```julia
using Plots
using ModelingToolkit
using Catalyst
```

## A simple biochemical example
$$
A + B \rightleftharpoons C
$$
BondGraphs can convert Catalyst reaction networks into a BondGraph


```julia
abc = @reaction_network ABC begin
    1, A + B --> C
end
```




\begin{align}
\require{mhchem}
\ce{ A + B &->[1] C}
\end{align}





```julia
bg_abc = BondGraph(abc)
```




    BondGraph ABC (5 Nodes, 4 Bonds)




```julia
plot(bg_abc, fontsize=12)
```




    

    




```julia
re, one, A, B, C, = nodes(bg_abc)
```




    5-element Vector{AbstractNode}:
     Re:R1
     𝟏_2
     Ce:A
     Ce:B
     Ce:C




```julia
constitutive_relations(A)
```




\begin{align}
0 =&  - E_1\left( t \right) + R T \log\left( K q\left( t \right) \right) \\
\frac{dq(t)}{dt} =& F_1\left( t \right)
\end{align}





```julia
constitutive_relations(re)
```




\begin{align}
0 =& F_1\left( t \right) + F_2\left( t \right) \\
0 =&  - r \left(  - e^{\frac{E_2\left( t \right)}{R T}} + e^{\frac{E_1\left( t \right)}{R T}} \right) + F_1\left( t \right)
\end{align}





```julia
eqs = constitutive_relations(bg_abc)
```




\begin{align}
\frac{dA_{+}q(t)}{dt} =&  - R1_{+}r \left(  - C_{+}K C_{+q}\left( t \right) + A_{+}K B_{+}K A_{+q}\left( t \right) B_{+q}\left( t \right) \right) \\
\frac{dB_{+}q(t)}{dt} =&  - R1_{+}r \left(  - C_{+}K C_{+q}\left( t \right) + A_{+}K B_{+}K A_{+q}\left( t \right) B_{+q}\left( t \right) \right) \\
\frac{dC_{+}q(t)}{dt} =& R1_{+}r \left(  - C_{+}K C_{+q}\left( t \right) + A_{+}K B_{+}K A_{+q}\left( t \right) B_{+q}\left( t \right) \right)
\end{align}





```julia
eqs = constitutive_relations(bg_abc; sub_defaults=true)
```




\begin{align}
\frac{dA_{+}q(t)}{dt} =&  - \left(  - C_{+q}\left( t \right) + A_{+q}\left( t \right) B_{+q}\left( t \right) \right) \\
\frac{dB_{+}q(t)}{dt} =&  - \left(  - C_{+q}\left( t \right) + A_{+q}\left( t \right) B_{+q}\left( t \right) \right) \\
\frac{dC_{+}q(t)}{dt} =&  - C_{+q}\left( t \right) + A_{+q}\left( t \right) B_{+q}\left( t \right)
\end{align}





```julia
tspan = (0., 3.)
u0 = [1, 2, 3]
sol = simulate(bg_abc, tspan; u0)
plot(sol, lw=2)
```




    

    



## Stoichiometric Example
$$
A + 2B \rightleftharpoons 3C
$$


```julia
a2b2c = @reaction_network A2B2C begin
    1, A + 2B --> 3C
end
```




\begin{align}
\require{mhchem}
\ce{ A + 2 B &->[1] 3 C}
\end{align}





```julia
bg_a2b2c = BondGraph(a2b2c)
plot(bg_a2b2c, fontsize=12)
```




    

    




```julia
eqs = constitutive_relations(bg_a2b2c; sub_defaults=true)
```




\begin{align}
\frac{dA_{+}q(t)}{dt} =&  - \left(  - \left( C_{+q}\left( t \right) \right)^{3} + \left( B_{+q}\left( t \right) \right)^{2} A_{+q}\left( t \right) \right) \\
\frac{dB_{+}q(t)}{dt} =&  - 2.0 \left(  - \left( C_{+q}\left( t \right) \right)^{3} + \left( B_{+q}\left( t \right) \right)^{2} A_{+q}\left( t \right) \right) \\
\frac{dC_{+}q(t)}{dt} =& 3 \left(  - \left( C_{+q}\left( t \right) \right)^{3} + \left( B_{+q}\left( t \right) \right)^{2} A_{+q}\left( t \right) \right)
\end{align}





```julia
tspan = (0., 1.)
u0 = [2, 1, 2]
sol = simulate(bg_a2b2c, tspan; u0)
plot(sol, lw=2)
```




    

    



## The reversible Michaelis-Menten model
$$
E + S \rightleftharpoons C \\ 
C \rightleftharpoons E + P
$$


```julia
rn_mm = @reaction_network MM_reversible begin
    (1, 1), E + S <--> C
    (1, 1), C <--> E + P
end
```




\begin{align}
\require{mhchem}
\ce{ E + S &<=>[1][1] C}\\
\ce{ C &<=>[1][1] E + P}
\end{align}




Chemostats are chemical species that are held constant throught the reaction (`Se` in bond graph terms)


```julia
bg_mm = BondGraph(rn_mm; chemostats=["S", "P"])
```




    BondGraph MM_reversible (10 Nodes, 10 Bonds)




```julia
plot(bg_mm, fontsize=12)
```




    

    




```julia
# Update default concentration for the substrate
bg_mm.S.xs = t -> 2
```




    #1 (generic function with 1 method)




```julia
constitutive_relations(bg_mm)
```




\begin{align}
\frac{dE_{+}q(t)}{dt} =& R2_{+}r \left( C_{+}K C_{+q}\left( t \right) - E_{+}K P_{+}K E_{+q}\left( t \right) \right) - R1_{+}r \left(  - C_{+}K C_{+q}\left( t \right) + 2 E_{+}K S_{+}K E_{+q}\left( t \right) \right) \\
\frac{dC_{+}q(t)}{dt} =& R1_{+}r \left(  - C_{+}K C_{+q}\left( t \right) + 2 E_{+}K S_{+}K E_{+q}\left( t \right) \right) - R2_{+}r \left( C_{+}K C_{+q}\left( t \right) - E_{+}K P_{+}K E_{+q}\left( t \right) \right)
\end{align}





```julia
tspan = (0., 3.)
u0 = [1,2]
sol = simulate(bg_mm, tspan; u0)
plot(sol, lw=2)
```




    

    




```julia
# Rerun simulation with a substrate which increases over time
bg_mm.S.xs = t -> 1 + t

tspan = (0., 10.)
u0 = [1,2]
sol = simulate(bg_mm, tspan; u0)
plot(sol, lw=2)
```




    

    



0799c04498d85126c73633435c56979f223c4b6a
ipynb
Jupyter Notebook
Biochemical/SimpleReactions.ipynb
jedforrest/BondGraphsTutorials
d93b12d9c3e31c9c6b53d9997d92c81587ee299e
MIT
Biochemical/SimpleReactions.ipynb
jedforrest/BondGraphsTutorials
d93b12d9c3e31c9c6b53d9997d92c81587ee299e
MIT
Biochemical/SimpleReactions.ipynb
jedforrest/BondGraphsTutorials
d93b12d9c3e31c9c6b53d9997d92c81587ee299e
MIT
Qwen/Qwen-72B
1. YES
2. YES