Generated Code
The following is python code generated by the CellML API from this CellML file. (Back to language selection)
The raw code is available.
# Size of variable arrays:
sizeAlgebraic = 24
sizeStates = 24
sizeConstants = 58
from math import *
from numpy import *
def createLegends():
legend_states = [""] * sizeStates
legend_rates = [""] * sizeStates
legend_algebraic = [""] * sizeAlgebraic
legend_voi = ""
legend_constants = [""] * sizeConstants
legend_voi = "time in component environment (minute)"
legend_states[0] = "Glc_x_0 in component Glc_x_0 (millimolar)"
legend_algebraic[0] = "inGlc in component inGlc (flux)"
legend_states[1] = "Glc_x in component Glc_x (millimolar)"
legend_algebraic[1] = "GlcTrans in component GlcTrans (flux)"
legend_states[2] = "Glc in component Glc (millimolar)"
legend_algebraic[2] = "HK in component HK (flux)"
legend_states[3] = "G6P in component G6P (millimolar)"
legend_algebraic[3] = "PGI in component PGI (flux)"
legend_algebraic[10] = "storage in component storage (flux)"
legend_states[4] = "F6P in component F6P (millimolar)"
legend_algebraic[4] = "PFK in component PFK (flux)"
legend_states[5] = "FBP in component FBP (millimolar)"
legend_algebraic[5] = "ALD in component ALD (flux)"
legend_states[6] = "GAP in component GAP (millimolar)"
legend_algebraic[7] = "GAPDH in component GAPDH (flux)"
legend_algebraic[6] = "TIM in component TIM (flux)"
legend_states[7] = "DHAP in component DHAP (millimolar)"
legend_algebraic[14] = "lpGlyc in component lpGlyc (flux)"
legend_states[8] = "BPG in component BPG (millimolar)"
legend_algebraic[8] = "lpPEP in component lpPEP (flux)"
legend_states[9] = "PEP in component PEP (millimolar)"
legend_algebraic[9] = "PK in component PK (flux)"
legend_states[10] = "Pyr in component Pyr (millimolar)"
legend_algebraic[11] = "PDC in component PDC (flux)"
legend_states[11] = "ATP in component ATP (millimolar)"
legend_algebraic[13] = "consum in component consum (flux)"
legend_algebraic[15] = "AK in component AK (flux)"
legend_states[12] = "ADP in component ADP (millimolar)"
legend_states[13] = "AMP in component AMP (millimolar)"
legend_states[14] = "CN_x_0 in component CN_x_0 (millimolar)"
legend_algebraic[23] = "inCN in component inCN (flux)"
legend_states[15] = "CN_x in component CN_x (millimolar)"
legend_algebraic[21] = "lacto in component lacto (flux)"
legend_states[16] = "ACA in component ACA (millimolar)"
legend_algebraic[18] = "difACA in component difACA (flux)"
legend_algebraic[12] = "ADH in component ADH (flux)"
legend_states[17] = "ACA_x in component ACA_x (millimolar)"
legend_algebraic[16] = "outACA in component outACA (flux)"
legend_states[18] = "EtOH in component EtOH (millimolar)"
legend_algebraic[19] = "difEtOH in component difEtOH (flux)"
legend_states[19] = "EtOH_x in component EtOH_x (millimolar)"
legend_algebraic[17] = "outEtOH in component outEtOH (flux)"
legend_states[20] = "Glyc in component Glyc (millimolar)"
legend_algebraic[22] = "difGlyc in component difGlyc (flux)"
legend_states[21] = "Glyc_x in component Glyc_x (millimolar)"
legend_algebraic[20] = "outGlyc in component outGlyc (flux)"
legend_states[22] = "NADH in component NADH (millimolar)"
legend_states[23] = "NAD in component NAD (millimolar)"
legend_constants[0] = "k0 in component model_parameters (first_order_rate_constant)"
legend_constants[1] = "y_vol in component model_parameters (dimensionless)"
legend_constants[2] = "K2_Glc in component GlcTrans (millimolar)"
legend_constants[3] = "K2_IG6P in component GlcTrans (millimolar)"
legend_constants[4] = "K2_IIG6P in component GlcTrans (millimolar)"
legend_constants[5] = "V_2m in component GlcTrans (flux)"
legend_constants[6] = "P2 in component GlcTrans (dimensionless)"
legend_constants[7] = "K3_DGlc in component HK (millimolar)"
legend_constants[8] = "K3_Glc in component HK (millimolar)"
legend_constants[9] = "K3_ATP in component HK (millimolar)"
legend_constants[10] = "V_3m in component HK (flux)"
legend_constants[11] = "V_4m in component PGI (flux)"
legend_constants[12] = "K4_G6P in component PGI (millimolar)"
legend_constants[13] = "K4_F6P in component PGI (millimolar)"
legend_constants[14] = "K4_eq in component PGI (dimensionless)"
legend_constants[15] = "K5 in component PFK (millimolar2)"
legend_constants[16] = "kappa5 in component PFK (dimensionless)"
legend_constants[17] = "V_5m in component PFK (flux)"
legend_constants[18] = "K6_eq in component ALD (millimolar)"
legend_constants[19] = "K6_FBP in component ALD (millimolar)"
legend_constants[20] = "K6_DHAP in component ALD (millimolar)"
legend_constants[21] = "K6_GAP in component ALD (millimolar)"
legend_constants[22] = "K6_IGAP in component ALD (millimolar)"
legend_constants[23] = "V_6r in component ALD (flux)"
legend_constants[24] = "V_6f in component ALD (flux)"
legend_constants[25] = "K7_eq in component TIM (dimensionless)"
legend_constants[26] = "K7_DHAP in component TIM (millimolar)"
legend_constants[27] = "K7_GAP in component TIM (millimolar)"
legend_constants[28] = "V_7m in component TIM (flux)"
legend_constants[29] = "K8_NAD in component GAPDH (millimolar)"
legend_constants[30] = "K8_NADH in component GAPDH (millimolar)"
legend_constants[31] = "K8_GAP in component GAPDH (millimolar)"
legend_constants[32] = "K8_BPG in component GAPDH (millimolar)"
legend_constants[33] = "K8_eq in component GAPDH (dimensionless)"
legend_constants[34] = "V_8m in component GAPDH (flux)"
legend_constants[35] = "k9f in component lpPEP (second_order_rate_constant)"
legend_constants[36] = "k9r in component lpPEP (second_order_rate_constant)"
legend_constants[37] = "K10_PEP in component PK (millimolar)"
legend_constants[38] = "K10_ADP in component PK (millimolar)"
legend_constants[39] = "V_10m in component PK (flux)"
legend_constants[40] = "K11 in component PDC (millimolar)"
legend_constants[41] = "V_11m in component PDC (flux)"
legend_constants[42] = "K12_NADH in component ADH (millimolar)"
legend_constants[43] = "K12_ACA in component ADH (millimolar)"
legend_constants[44] = "V_12m in component ADH (flux)"
legend_constants[45] = "k13 in component difEtOH (first_order_rate_constant)"
legend_constants[46] = "K15_NADH in component lpGlyc (millimolar)"
legend_constants[47] = "K15_INADH in component lpGlyc (millimolar)"
legend_constants[48] = "K15_INAD in component lpGlyc (millimolar)"
legend_constants[49] = "K15_DHAP in component lpGlyc (millimolar)"
legend_constants[50] = "V_15m in component lpGlyc (flux)"
legend_constants[51] = "k16 in component difGlyc (first_order_rate_constant)"
legend_constants[52] = "k18 in component difACA (first_order_rate_constant)"
legend_constants[53] = "k20 in component lacto (second_order_rate_constant)"
legend_constants[54] = "k22 in component storage (second_order_rate_constant)"
legend_constants[55] = "k23 in component consum (first_order_rate_constant)"
legend_constants[56] = "k24f in component AK (second_order_rate_constant)"
legend_constants[57] = "k24r in component AK (second_order_rate_constant)"
legend_rates[0] = "d/dt Glc_x_0 in component Glc_x_0 (millimolar)"
legend_rates[1] = "d/dt Glc_x in component Glc_x (millimolar)"
legend_rates[2] = "d/dt Glc in component Glc (millimolar)"
legend_rates[3] = "d/dt G6P in component G6P (millimolar)"
legend_rates[4] = "d/dt F6P in component F6P (millimolar)"
legend_rates[5] = "d/dt FBP in component FBP (millimolar)"
legend_rates[6] = "d/dt GAP in component GAP (millimolar)"
legend_rates[7] = "d/dt DHAP in component DHAP (millimolar)"
legend_rates[8] = "d/dt BPG in component BPG (millimolar)"
legend_rates[9] = "d/dt PEP in component PEP (millimolar)"
legend_rates[10] = "d/dt Pyr in component Pyr (millimolar)"
legend_rates[11] = "d/dt ATP in component ATP (millimolar)"
legend_rates[12] = "d/dt ADP in component ADP (millimolar)"
legend_rates[13] = "d/dt AMP in component AMP (millimolar)"
legend_rates[14] = "d/dt CN_x_0 in component CN_x_0 (millimolar)"
legend_rates[15] = "d/dt CN_x in component CN_x (millimolar)"
legend_rates[16] = "d/dt ACA in component ACA (millimolar)"
legend_rates[17] = "d/dt ACA_x in component ACA_x (millimolar)"
legend_rates[18] = "d/dt EtOH in component EtOH (millimolar)"
legend_rates[19] = "d/dt EtOH_x in component EtOH_x (millimolar)"
legend_rates[20] = "d/dt Glyc in component Glyc (millimolar)"
legend_rates[21] = "d/dt Glyc_x in component Glyc_x (millimolar)"
legend_rates[22] = "d/dt NADH in component NADH (millimolar)"
legend_rates[23] = "d/dt NAD in component NAD (millimolar)"
return (legend_states, legend_algebraic, legend_voi, legend_constants)
def initConsts():
constants = [0.0] * sizeConstants; states = [0.0] * sizeStates;
states[0] = 24.0
states[1] = 6.7
states[2] = 0.573074
states[3] = 4.2
states[4] = 0.49
states[5] = 4.64
states[6] = 0.115
states[7] = 2.95
states[8] = 0.00027
states[9] = 0.04
states[10] = 8.7
states[11] = 2.1
states[12] = 1.5
states[13] = 0.33
states[14] = 5.60
states[15] = 5.20358
states[16] = 1.48153
states[17] = 1.28836
states[18] = 19.2379
states[19] = 16.4514
states[20] = 4.196
states[21] = 1.68478
states[22] = 0.33
states[23] = 0.65
constants[0] = 0.048
constants[1] = 59.0
constants[2] = 1.7
constants[3] = 1.2
constants[4] = 7.2
constants[5] = 1014.96
constants[6] = 1.0
constants[7] = 0.37
constants[8] = 0.0
constants[9] = 0.1
constants[10] = 5.7547
constants[11] = 496.042
constants[12] = 0.8
constants[13] = 0.15
constants[14] = 0.13
constants[15] = 0.021
constants[16] = 0.15
constants[17] = 45.4327
constants[18] = 0.081
constants[19] = 0.3
constants[20] = 2.0
constants[21] = 4.0
constants[22] = 10.0
constants[23] = 1.10391E4
constants[24] = 2.20782E3
constants[25] = 0.055
constants[26] = 1.23
constants[27] = 1.27
constants[28] = 1.16365E2
constants[29] = 0.1
constants[30] = 0.06
constants[31] = 0.6
constants[32] = 0.01
constants[33] = 0.0055
constants[34] = 8.33858E2
constants[35] = 4.43866E5
constants[36] = 1.52862E3
constants[37] = 0.2
constants[38] = 0.17
constants[39] = 3.43096E2
constants[40] = 0.3
constants[41] = 5.31328E1
constants[42] = 0.1
constants[43] = 0.71
constants[44] = 8.98023E1
constants[45] = 16.72
constants[46] = 0.13
constants[47] = 0.034
constants[48] = 0.13
constants[49] = 25.0
constants[50] = 8.14797E1
constants[51] = 1.9
constants[52] = 24.7
constants[53] = 2.83828E-3
constants[54] = 2.25932
constants[55] = 3.20760
constants[56] = 4.32900E2
constants[57] = 1.33333E2
return (states, constants)
def computeRates(voi, states, constants):
rates = [0.0] * sizeStates; algebraic = [0.0] * sizeAlgebraic
algebraic[0] = constants[0]*(states[0]-states[1])
rates[0] = algebraic[0]
algebraic[1] = (constants[5]/constants[1])*((states[1]/constants[2])/(1.00000+states[1]/constants[2]+((constants[6]*(states[1]/constants[2])+1.00000)/(constants[6]*(states[2]/constants[2])+1.00000))*(1.00000+states[2]/constants[2])+states[3]/constants[3]+(states[2]*states[3])/(constants[2]*constants[4])))-(constants[5]/constants[1])*((states[2]/constants[2])/(1.00000+states[2]/constants[2]+((constants[6]*(states[2]/constants[2])+1.00000)/(constants[6]*(states[1]/constants[2])+1.00000))*(1.00000+states[1]/constants[2])+states[3]/constants[3]+(states[2]*states[3])/(constants[2]*constants[4])))
rates[1] = algebraic[0]+algebraic[1]
algebraic[2] = (constants[10]*states[11]*states[2])/(constants[7]*constants[9]+constants[8]*states[11]+constants[9]*states[2]+states[2]*states[11])
rates[2] = algebraic[1]+algebraic[2]
algebraic[3] = (constants[11]*states[3])/(constants[12]+states[3]+(constants[12]/constants[13])*states[4])-(constants[11]*(states[4]/constants[14]))/(constants[12]+states[3]+(constants[12]/constants[13])*states[4])
algebraic[4] = (constants[17]*(power(states[4], 2.00000)))/(constants[15]*(1.00000+constants[16]*(states[11]/states[13]))+power(states[4], 2.00000))
rates[4] = algebraic[3]+algebraic[4]
algebraic[5] = (constants[24]*states[5])/(constants[19]+states[5]+(states[6]*constants[20]*constants[24])/(constants[18]*constants[23])+(states[7]*constants[21]*constants[24])/(constants[18]*constants[23])+(states[5]*states[6])/constants[22]+(states[6]*states[7]*constants[24])/(constants[18]*constants[23]))-(constants[24]*((states[6]*states[7])/constants[18]))/(constants[19]+states[5]+(states[6]*constants[20]*constants[24])/(constants[18]*constants[23])+(states[7]*constants[21]*constants[24])/(constants[18]*constants[23])+(states[5]*states[6])/constants[22]+(states[6]*states[7]*constants[24])/(constants[18]*constants[23]))
rates[5] = algebraic[4]+algebraic[5]
algebraic[7] = (constants[34]*states[6]*states[23])/(constants[31]*constants[29]*(1.00000+states[6]/constants[31]+states[8]/constants[32])*(1.00000+states[23]/constants[29]+states[22]/constants[30]))-(constants[34]*((states[8]*states[22])/constants[33]))/(constants[31]*constants[29]*(1.00000+states[6]/constants[31]+states[8]/constants[32])*(1.00000+states[23]/constants[29]+states[22]/constants[30]))
algebraic[6] = (constants[28]*states[7])/(constants[26]+states[7]+(constants[26]/constants[27])*states[6])-(constants[28]*(states[6]/constants[25]))/(constants[26]+states[7]+(constants[26]/constants[27])*states[6])
rates[6] = algebraic[5]+algebraic[6]+algebraic[7]
algebraic[8] = constants[35]*states[8]*states[12]-constants[36]*states[9]*states[11]
rates[8] = algebraic[7]+algebraic[8]
algebraic[9] = (constants[39]*states[12]*states[9])/((constants[37]+states[9])*(constants[38]+states[12]))
rates[9] = algebraic[8]+algebraic[9]
algebraic[10] = constants[54]*states[11]*states[3]
rates[3] = algebraic[2]+algebraic[3]+algebraic[10]
algebraic[11] = (constants[41]*states[10])/(constants[40]+states[10])
rates[10] = algebraic[9]+algebraic[11]
algebraic[14] = (constants[50]*states[7])/(constants[49]*(1.00000+(constants[47]/states[22])*(1.00000+states[23]/constants[48]))+states[7]*(1.00000+(constants[46]/states[22])*(1.00000+states[23]/constants[48])))
rates[7] = algebraic[5]+algebraic[6]+algebraic[14]
algebraic[13] = constants[55]*states[11]
algebraic[15] = constants[56]*states[13]*states[11]-constants[57]*(power(states[12], 2.00000))
rates[11] = algebraic[4]+algebraic[9]+algebraic[8]+algebraic[10]+algebraic[2]+algebraic[13]+algebraic[15]
rates[12] = algebraic[4]+algebraic[9]+algebraic[8]+algebraic[10]+algebraic[2]+algebraic[13]+algebraic[15]
rates[13] = algebraic[15]
algebraic[12] = (constants[44]*states[16]*states[22])/((constants[42]+states[22])*(constants[43]+states[16]))
rates[22] = algebraic[7]+algebraic[12]+algebraic[14]
rates[23] = algebraic[7]+algebraic[12]+algebraic[14]
algebraic[18] = (constants[52]/constants[1])*(states[16]-states[17])
rates[16] = algebraic[18]+algebraic[11]+algebraic[12]
algebraic[19] = (constants[45]/constants[1])*(states[18]-states[19])
rates[18] = algebraic[19]+algebraic[12]
algebraic[17] = constants[0]*states[19]
rates[19] = algebraic[19]+algebraic[17]
algebraic[21] = constants[53]*states[17]*states[15]
algebraic[16] = constants[0]*states[17]
rates[17] = algebraic[21]+algebraic[18]+algebraic[16]
algebraic[22] = (constants[51]/constants[1])*(states[20]-states[21])
rates[20] = algebraic[22]+algebraic[14]
algebraic[20] = constants[0]*states[21]
rates[21] = algebraic[22]+algebraic[20]
algebraic[23] = constants[0]*(states[14]-states[15])
rates[14] = algebraic[23]
rates[15] = algebraic[21]+algebraic[23]
return(rates)
def computeAlgebraic(constants, states, voi):
algebraic = array([[0.0] * len(voi)] * sizeAlgebraic)
states = array(states)
voi = array(voi)
algebraic[0] = constants[0]*(states[0]-states[1])
algebraic[1] = (constants[5]/constants[1])*((states[1]/constants[2])/(1.00000+states[1]/constants[2]+((constants[6]*(states[1]/constants[2])+1.00000)/(constants[6]*(states[2]/constants[2])+1.00000))*(1.00000+states[2]/constants[2])+states[3]/constants[3]+(states[2]*states[3])/(constants[2]*constants[4])))-(constants[5]/constants[1])*((states[2]/constants[2])/(1.00000+states[2]/constants[2]+((constants[6]*(states[2]/constants[2])+1.00000)/(constants[6]*(states[1]/constants[2])+1.00000))*(1.00000+states[1]/constants[2])+states[3]/constants[3]+(states[2]*states[3])/(constants[2]*constants[4])))
algebraic[2] = (constants[10]*states[11]*states[2])/(constants[7]*constants[9]+constants[8]*states[11]+constants[9]*states[2]+states[2]*states[11])
algebraic[3] = (constants[11]*states[3])/(constants[12]+states[3]+(constants[12]/constants[13])*states[4])-(constants[11]*(states[4]/constants[14]))/(constants[12]+states[3]+(constants[12]/constants[13])*states[4])
algebraic[4] = (constants[17]*(power(states[4], 2.00000)))/(constants[15]*(1.00000+constants[16]*(states[11]/states[13]))+power(states[4], 2.00000))
algebraic[5] = (constants[24]*states[5])/(constants[19]+states[5]+(states[6]*constants[20]*constants[24])/(constants[18]*constants[23])+(states[7]*constants[21]*constants[24])/(constants[18]*constants[23])+(states[5]*states[6])/constants[22]+(states[6]*states[7]*constants[24])/(constants[18]*constants[23]))-(constants[24]*((states[6]*states[7])/constants[18]))/(constants[19]+states[5]+(states[6]*constants[20]*constants[24])/(constants[18]*constants[23])+(states[7]*constants[21]*constants[24])/(constants[18]*constants[23])+(states[5]*states[6])/constants[22]+(states[6]*states[7]*constants[24])/(constants[18]*constants[23]))
algebraic[7] = (constants[34]*states[6]*states[23])/(constants[31]*constants[29]*(1.00000+states[6]/constants[31]+states[8]/constants[32])*(1.00000+states[23]/constants[29]+states[22]/constants[30]))-(constants[34]*((states[8]*states[22])/constants[33]))/(constants[31]*constants[29]*(1.00000+states[6]/constants[31]+states[8]/constants[32])*(1.00000+states[23]/constants[29]+states[22]/constants[30]))
algebraic[6] = (constants[28]*states[7])/(constants[26]+states[7]+(constants[26]/constants[27])*states[6])-(constants[28]*(states[6]/constants[25]))/(constants[26]+states[7]+(constants[26]/constants[27])*states[6])
algebraic[8] = constants[35]*states[8]*states[12]-constants[36]*states[9]*states[11]
algebraic[9] = (constants[39]*states[12]*states[9])/((constants[37]+states[9])*(constants[38]+states[12]))
algebraic[10] = constants[54]*states[11]*states[3]
algebraic[11] = (constants[41]*states[10])/(constants[40]+states[10])
algebraic[14] = (constants[50]*states[7])/(constants[49]*(1.00000+(constants[47]/states[22])*(1.00000+states[23]/constants[48]))+states[7]*(1.00000+(constants[46]/states[22])*(1.00000+states[23]/constants[48])))
algebraic[13] = constants[55]*states[11]
algebraic[15] = constants[56]*states[13]*states[11]-constants[57]*(power(states[12], 2.00000))
algebraic[12] = (constants[44]*states[16]*states[22])/((constants[42]+states[22])*(constants[43]+states[16]))
algebraic[18] = (constants[52]/constants[1])*(states[16]-states[17])
algebraic[19] = (constants[45]/constants[1])*(states[18]-states[19])
algebraic[17] = constants[0]*states[19]
algebraic[21] = constants[53]*states[17]*states[15]
algebraic[16] = constants[0]*states[17]
algebraic[22] = (constants[51]/constants[1])*(states[20]-states[21])
algebraic[20] = constants[0]*states[21]
algebraic[23] = constants[0]*(states[14]-states[15])
return algebraic
def solve_model():
"""Solve model with ODE solver"""
from scipy.integrate import ode
# Initialise constants and state variables
(init_states, constants) = initConsts()
# Set timespan to solve over
voi = linspace(0, 10, 500)
# Construct ODE object to solve
r = ode(computeRates)
r.set_integrator('vode', method='bdf', atol=1e-06, rtol=1e-06, max_step=1)
r.set_initial_value(init_states, voi[0])
r.set_f_params(constants)
# Solve model
states = array([[0.0] * len(voi)] * sizeStates)
states[:,0] = init_states
for (i,t) in enumerate(voi[1:]):
if r.successful():
r.integrate(t)
states[:,i+1] = r.y
else:
break
# Compute algebraic variables
algebraic = computeAlgebraic(constants, states, voi)
return (voi, states, algebraic)
def plot_model(voi, states, algebraic):
"""Plot variables against variable of integration"""
import pylab
(legend_states, legend_algebraic, legend_voi, legend_constants) = createLegends()
pylab.figure(1)
pylab.plot(voi,vstack((states,algebraic)).T)
pylab.xlabel(legend_voi)
pylab.legend(legend_states + legend_algebraic, loc='best')
pylab.show()
if __name__ == "__main__":
(voi, states, algebraic) = solve_model()
plot_model(voi, states, algebraic)