#!/usr/bin/env python
# coding: utf-8

# ***Note:*** 
# 
# - This example is discussed in detail by <cite data-cite="GawPan20X">Gawthrop and Pan (2020)</cite> available [here](https://arxiv.org/abs/2009.02217).
# 
# - This is the Redox.ipynb notebook. The PDF version is available [here](https://github.com/gawthrop/GawPan20/blob/main/Redox.pdf).
# 
# 

# # Introduction
# Redox reactions and proton pumps play a fundamental role in Biology. This note illustrates this using a bond graph model of complex I of the mitochondrial electron transport chain.
# 

# In[1]:


## Some useful imports
import BondGraphTools as bgt
print('BGT version is',bgt.version)
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt

## Stoichiometric analysis
import stoich as st

## SVG
import svgBondGraph as sbg

## Display (eg disp.SVG(), disp.
import IPython.display as disp

## Potential data
import phiData

quiet = True
TranslateSVG = False


# # Redox reaction
# A key energy-generating redox reaction that within the mitochondrial respiratory chain is
# \begin{equation}\label{eq:NAD}
# \mathrm{NADH + Q + H^+ \rightleftharpoons NAD+ +  QH2 }
# \end{equation}
# This reaction can be divided into the half reactions:
# \begin{align}
# \mathrm{NADH} &\rightleftharpoons \mathrm{2 e_1^- + H^+ + NAD^+}\label{eq:NAD_1}\\
# \mathrm{2 e_2^- + 2 H^+ + Q} &\rightleftharpoons \mathrm{QH_2}\label{eq:NAD_2}
# \end{align}
# A bond graph representation of this decomposition is given  below.

# ## Bond graph
# C:E1 and C:E2 represent the electron potentials and the other The C components represent the species; the two Re components the two half reactions. 

# In[2]:


## Redox reaction
if TranslateSVG:
    sbg.model('Redox_abg.svg',quiet=False)
import Redox_abg
disp.SVG('Redox_abg.svg')


# ## Stoichiometry

# In[3]:


## Stoichiometry
linear = ['E1','E2']
s = st.stoich(Redox_abg.model(),linear=linear,quiet=quiet)
chemostats = ['NADH','NAD','Q','QH2','H']
sc = st.statify(s,chemostats=chemostats)


# In[4]:


## Stoichiometric matrix
disp.Latex(st.sprintl(s,'N'))


# ## Reactions
# These are automatically generated from the bond graph

# In[5]:


## Reactions
disp.Latex(st.sprintrl(s,chemformula=True,all=True))


# ## Potentials
# The reaction (Faraday Equivalent) potentials are computed from tables gleaned from the literature

# In[6]:


## Standard potetials
phi_Std = phiData.phi_Std()

## Typical concentrations
conc = phiData.ParRubXu16_conc()

## From BazBeaVin16
conc['Q'] = conc['QH2'] = 1e-2

## pH 7
conc['H'] = 1e-7

## Table for paper and put values in to phi_NADH etc.
print('%% Table')
ch='\ch'
l='{'
r='}'
eol = r'\\'
phi_std = {}
for spec in ['NAD','NADH','Q','QH2','H']:
    phi0 = phi_Std[spec]
    con = conc[spec]
    phi_std_spec = phi0 + st.V_N()*np.log(con)
    phi_std[spec] = phi_std_spec
    #print(f'phi_Std_{spec} = {1000*phi0:0.0f}, phi_{spec} = {1000*phi_std[spec]:.0f}, conc_{spec}={conc[spec]}')
    print(f'{ch}{l}{spec}{r} & {1000*phi0:.0f} & {con:1.2e} & {1000*phi_std_spec:.0f}{eol}')
    exec(f'phi_{spec} = {phi_std_spec}')
  

## Print the worked example for the paper.
print('\n%% Equations')
E1 = 0.5*(phi_NADH - phi_NAD - phi_H)
E2 = 0.5*(phi_QH2 - phi_Q - 2*phi_H)
print(f'E1 = 0.5({1000*phi_NADH:.0f} - {1000*phi_NAD:.0f} - {1000*phi_H:.0f}) = {1000*E1:.0f} mV')
print(f'E2 = 0.5({1000*phi_QH2:.0f} - {1000*phi_Q:.0f} - 2x{1000*phi_H:.0f}) = {1000*E2:.0f} mV')
print(f'E1-E2 = {1000*(E1-E2):.0f} mV')
print(f'PMF = {1000*(E1-E2)/2:.0f} mV')


# 
# # Proton pump
# The redox reaction of complex I drives a proton pump.

# ## Bond graph
# C:E1 and C:E2 correspond to the redox reaction and provide the potential to drive protons in the interior $H^+_i$ to the exterior $H^+_e$ of the mitochondrial membrane. The protons have both electrical and chemical potential.

# In[7]:


## Proton pump
if TranslateSVG:
    sbg.model('ProtonPump_abg.svg')
import ProtonPump_abg
disp.SVG('ProtonPump_abg.svg')


# ## Stoichiometry

# In[8]:


## Stoichiometry
linear = ['E1','E2','Ei','Ee']
s = st.stoich(ProtonPump_abg.model(),linear=linear,quiet=quiet)
chemostats = ['E1','E2','Ei','Ee','Hi','He']
sc = st.statify(s,chemostats=chemostats)


# ## Reactions
# These are automatically generated from the bond graph

# In[9]:


## Reactions
disp.Latex(st.sprintrl(s,chemformula=True,all=True))


# In[10]:


## Flows
disp.Latex(st.sprintvl(s))


# # Complex I
# The model of mitochondrial complex I comprises two modules: the redox reaction and the proton pump. 

# In[11]:


## Complex I
if TranslateSVG:
    sbg.model('ComplexI_abg.svg')
import ComplexI_abg
disp.SVG('ComplexI_abg.svg')


# ## Stoichiometry

# In[12]:


## Stoichiometry
linear = ['E1','E2','pp_Ei','pp_Ee']
s = st.stoich(ComplexI_abg.model(),linear=linear,quiet=quiet)
print(s['species'])
chemostats = ['pp_Ee', 'pp_Ei', 'pp_He', 'pp_Hi', 'rr_H', 'rr_NAD', 'rr_NADH', 'rr_Q', 'rr_QH2']
sc = st.statify(s,chemostats=chemostats)
#print(s['species'])
#disp.Latex(st.sprint(s0,'K'))
#print(st.sprints(s))


# ## Reactions
# These are automatically generated from the bond graph

# In[13]:


## Reactions
disp.Latex(st.sprintrl(s,chemformula=True,all=True))


# In[14]:


## Flows
disp.Latex(st.sprintvl(s))


# In[15]:


## Path
sp = st.path(s,sc)
## Reactions
disp.Latex(st.sprintrl(sp,chemformula=True,all=True))