- Author:
- Shelley Fong <s.fong@auckland.ac.nz>
- Date:
- 2021-10-28 14:44:03+13:00
- Desc:
- Adding script to find BG parameters
- Permanent Source URI:
- https://staging.physiomeproject.org/workspace/7a5/rawfile/84946bde4e7490903b7cd9807abe9ff235101cf8/parameter_finder/find_bg_parameters_SERCA.py
# This script calculates the bond graph parameters for the SERCA model of
# Tran et al. (2009) by using the kinetic parameters and stoichiometric
# matrix.
# translated from MATLAB to Python by SF using Pan's original code
import os
import csv
import json
import math
import numpy as np
import sympy
from scipy.linalg import null_space
def read_IDs(path):
data = []
with open(path,'r') as f:
reader = csv.reader(f)
for row in reader:
data.append(row[0])
f.close()
return data
def load_matrix(stoich_path):
matrix = []
with open(stoich_path,'r') as f:
reader = csv.reader(f,delimiter=',')
for row in reader:
matrix.append([int(r) for r in row])
f.close()
return matrix
def kinetic_parameters(M, dims, V):
# R1-2
k_12p = 25900
# R2-4
K_dCai = 0.91
k_24m = 259e8
k_24p = k_24m/pow(K_dCai,2)
# R2-2a
K_dH1 = 1.09e-5
k_22am = 259e8
k_22ap = k_22am/K_dH1
# R4-5
K_dHi = 3.54e-3
k_45p = 259e8
k_45m = k_45p/K_dHi
# R5-6
k_56p = 2540
k_56m = 67200
# R6-8
K_dCasr = 2.24
k_68m = 259e8
k_68p = k_68m*pow(K_dCasr,2)
# R8-9
K_dHsr = 1.05e-8
k_89m = 259e8
k_89p = k_89m/K_dHsr
# R9-10
K_dH = 7.24e-5
k_910p = 259e8
k_910m = k_910p/K_dH
# R10-1
k_101p = 20.5
k_101m = 149
G_0_MgATP = 11900
R = 8.314
T = 310
K_MgATP = math.exp(-G_0_MgATP/(R*T))*pow(10,6)
# Calculate remaining parameter using detailed balance
k_12m = (k_12p*k_24p*k_45p*k_56p*k_68p*k_89p*k_910p*k_101p)/(k_24m*k_45m*k_56m*k_68m*k_89m*k_910m*k_101m*K_MgATP)
## Calculate bond graph constants from kinetic parameters
# Note: units of kappa are fmol/s, units of K are fmol^-1
K_C = [K_MgATP, 1]
k = [k_12p, k_24p, k_22ap, k_45p, k_56p, k_68p, k_89p, k_910p, k_101p, k_12m, k_24m, k_22am, k_45m, k_56m, k_68m, k_89m, k_910m, k_101m]
# constraints
N_cT = np.zeros([2,len(M[0])])
# Add constraint for MgATP hydrolysis
N_cT[0][21] = 1
N_cT[0][22] = -1
N_cT[0][23] = -1
N_cT[0][18] = -1
# Add constraint for calcium transport
N_cT[1][19] = 1
N_cT[1][20] = -1
W = [1]*18 + [ V['V_isr'], V['V_myo'], V['V_sr'], V['V_myo'], V['V_myo'], V['V_myo']]
return (k, N_cT, K_C, W)
if __name__ == "__main__":
## booleans
write_parameters_file = True
include_constraints = True
## Set directories
current_dir = os.getcwd()
data_dir = current_dir + '\data'
output_dir = current_dir + '\output'
modname = os.path.dirname(current_dir).split('\\')[-1].split('_')[-1]
## Define the volumes of each compartment (units pL)
W_i = 38.0 # Intracellular volume
W_sr = 2.28 # SR volume
W_isr = W_i + W_sr # Intracellular volume + SR volume
V = dict()
V['V_myo'] = W_i
V['V_sr'] = W_sr
V['V_isr'] = W_isr
## Load forward matrix
stoich_path = data_dir + '\\all_forward_matrix.txt'
N_f = load_matrix(stoich_path)
## Load reverse matrix
stoich_path = data_dir + '\\all_reverse_matrix.txt'
N_r = load_matrix(stoich_path)
N_fT = np.transpose(N_f)
N_rT = np.transpose(N_r)
## Calculate stoichiometric matrix
N = [[N_r[j][i] - N_f[j][i] for i in range(len(N_f[0]))] for j in range(len(N_f))]
N_T = [[N_rT[j][i] - N_fT[j][i] for i in range(len(N_fT[0]))] for j in range(len(N_fT))]
num_rows = len(N)
num_cols = len(N[0])
dims = dict()
dims['num_rows'] = num_rows
dims['num_cols'] = num_cols
I = np.identity(num_cols)
M = np.append(np.append(I, N_fT,1), np.append(I, N_rT,1),0)
[k_kinetic, N_cT, K_C, W] = kinetic_parameters(M, dims, V)
if not include_constraints:
N_cT = []
try:
M = np.append(M, N_cT,0)
k = np.append(k_kinetic, K_C, 0)
except:
k = k_kinetic
# FIND PARAMETERS HERE
lambda_expo = np.matmul(np.linalg.pinv(M), [math.log(ik) for ik in k])
lambdaW = [math.exp(l) for l in lambda_expo]
# Check that kinetic parameters are reproduced by bond graph parameters
k_est = np.matmul(M,[math.log(k) for k in lambdaW])
k_est = [math.exp(k) for k in k_est]
diff = [(k[i] - k_est[i])/k[i] for i in range(len(k))]
error = np.sum([abs(d) for d in diff])
# Check that there is a detailed balance constraint
Z = null_space(np.transpose(M))
N_rref = sympy.Matrix(N).rref()
R_mat = null_space(N)
kf = k_kinetic[:num_cols]
kr = k_kinetic[num_cols:]
K_eq = [kf[i]/kr[i] for i in range(len(kr))]
zero_est = np.matmul(np.transpose(R_mat),K_eq)
zero_est_log = np.matmul(np.transpose(R_mat),[math.log(k) for k in K_eq])
lambdak = [lambdaW[i]/W[i] for i in range(len(W))]
kappa = lambdak[:len(N[0])]
K = lambdak[len(N[0]):]
rxnID = read_IDs('data\\rxnID.txt')
Kname = read_IDs('data\\Kname.txt')
# ### print outputs ###
for ik in range(len(kappa)):
print('var kappa_%s: fmol_per_sec {init: %g, pub: out};' %(rxnID[ik],kappa[ik]))
for ik in range(len(Kname)):
print('var K_%s: per_fmol {init: %g, pub: out};' %(Kname[ik],K[ik]))
file = open(output_dir + '/all_parameters_out.json', 'w')
data = { "K": K, "kappa": kappa, "k_kinetic": k_kinetic }
json.dump(data, file)
