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 = 21
sizeStates = 20
sizeConstants = 50
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 (second)"
legend_algebraic[5] = "LRG in component beta_adrenergic_receptor_module (micromolar)"
legend_algebraic[7] = "RG in component beta_adrenergic_receptor_module (micromolar)"
legend_algebraic[0] = "LR in component beta_adrenergic_receptor_module (micromolar)"
legend_states[0] = "beta1_AR in component beta_adrenergic_receptor_module (micromolar)"
legend_states[1] = "Gs in component beta_adrenergic_receptor_module (micromolar)"
legend_states[2] = "beta1_ARact in component beta_adrenergic_receptor_module (micromolar)"
legend_states[3] = "beta1_AR_S464 in component beta_adrenergic_receptor_module (micromolar)"
legend_states[4] = "beta1_AR_S301 in component beta_adrenergic_receptor_module (micromolar)"
legend_states[5] = "L in component beta_adrenergic_receptor_module (micromolar)"
legend_constants[0] = "Ltotmax in component beta_adrenergic_receptor_module (micromolar)"
legend_constants[1] = "KL in component beta_adrenergic_receptor_module (micromolar)"
legend_constants[2] = "KR in component beta_adrenergic_receptor_module (micromolar)"
legend_constants[3] = "KC in component beta_adrenergic_receptor_module (micromolar)"
legend_constants[4] = "Gstot in component beta_adrenergic_receptor_module (micromolar)"
legend_constants[5] = "k_betaARK_plus in component beta_adrenergic_receptor_module (first_order_rate_constant)"
legend_constants[6] = "k_betaARK_minus in component beta_adrenergic_receptor_module (first_order_rate_constant)"
legend_constants[7] = "k_PKA_plus in component beta_adrenergic_receptor_module (second_order_rate_constant)"
legend_constants[8] = "k_PKA_minus in component beta_adrenergic_receptor_module (first_order_rate_constant)"
legend_states[6] = "Gs_beta_gamma in component Gs_activation_module (micromolar)"
legend_algebraic[9] = "PKACI in component PKA_activation_module (micromolar)"
legend_states[7] = "Gs_alpha_GTPtot in component Gs_activation_module (micromolar)"
legend_states[8] = "Gs_alpha_GDP in component Gs_activation_module (micromolar)"
legend_constants[9] = "k_gact in component Gs_activation_module (first_order_rate_constant)"
legend_constants[10] = "k_hyd in component Gs_activation_module (first_order_rate_constant)"
legend_constants[11] = "k_reassoc in component Gs_activation_module (second_order_rate_constant)"
legend_states[9] = "cAMPtot in component cyclic_AMP_metabolism_module (micromolar)"
legend_states[10] = "AC in component cyclic_AMP_metabolism_module (micromolar)"
legend_states[11] = "Gs_alpha_GTP in component cyclic_AMP_metabolism_module (micromolar)"
legend_algebraic[2] = "Gs_alpha_GTP_AC in component cyclic_AMP_metabolism_module (micromolar)"
legend_states[12] = "PDE in component cyclic_AMP_metabolism_module (micromolar)"
legend_constants[12] = "ACtot in component cyclic_AMP_metabolism_module (micromolar)"
legend_constants[13] = "ATP in component cyclic_AMP_metabolism_module (micromolar)"
legend_constants[14] = "PDEtot in component cyclic_AMP_metabolism_module (micromolar)"
legend_algebraic[1] = "PDEinhib in component cyclic_AMP_metabolism_module (micromolar)"
legend_constants[15] = "IBMXtot in component cyclic_AMP_metabolism_module (micromolar)"
legend_states[13] = "IBMX in component cyclic_AMP_metabolism_module (micromolar)"
legend_constants[16] = "fsktot in component cyclic_AMP_metabolism_module (micromolar)"
legend_states[14] = "fsk in component cyclic_AMP_metabolism_module (micromolar)"
legend_algebraic[6] = "fsk_AC in component cyclic_AMP_metabolism_module (micromolar)"
legend_constants[17] = "kAC_basal in component cyclic_AMP_metabolism_module (first_order_rate_constant)"
legend_constants[18] = "kAC_fsk in component cyclic_AMP_metabolism_module (first_order_rate_constant)"
legend_constants[19] = "k_PDE in component cyclic_AMP_metabolism_module (first_order_rate_constant)"
legend_constants[20] = "kAC_Gs_alpha in component cyclic_AMP_metabolism_module (first_order_rate_constant)"
legend_constants[21] = "Km_basal in component cyclic_AMP_metabolism_module (micromolar)"
legend_constants[22] = "Km_PDE in component cyclic_AMP_metabolism_module (micromolar)"
legend_constants[23] = "Km_fsk in component cyclic_AMP_metabolism_module (micromolar)"
legend_constants[24] = "K_fsk in component cyclic_AMP_metabolism_module (micromolar)"
legend_constants[25] = "Km_Gs_alpha_GTP in component cyclic_AMP_metabolism_module (micromolar)"
legend_constants[26] = "K_Gs_alpha in component cyclic_AMP_metabolism_module (micromolar)"
legend_constants[27] = "K_IBMX in component cyclic_AMP_metabolism_module (micromolar)"
legend_algebraic[10] = "cAMP in component PKA_activation_module (micromolar)"
legend_algebraic[11] = "PKACII in component PKA_activation_module (micromolar)"
legend_algebraic[12] = "PKAtemp in component PKA_activation_module (micromolar)"
legend_algebraic[13] = "ARCI in component PKA_activation_module (micromolar)"
legend_algebraic[14] = "A2RCI in component PKA_activation_module (micromolar)"
legend_algebraic[15] = "A2RI in component PKA_activation_module (micromolar)"
legend_algebraic[16] = "A2RCII in component PKA_activation_module (micromolar)"
legend_algebraic[17] = "A2RII in component PKA_activation_module (micromolar)"
legend_algebraic[18] = "ARCII in component PKA_activation_module (micromolar)"
legend_constants[28] = "Ki_pki in component PKA_activation_module (micromolar)"
legend_constants[29] = "PKAItot in component PKA_activation_module (micromolar)"
legend_constants[30] = "PKAIItot in component PKA_activation_module (micromolar)"
legend_constants[31] = "PKItot in component PKA_activation_module (micromolar)"
legend_constants[32] = "KA in component PKA_activation_module (micromolar)"
legend_constants[33] = "KB in component PKA_activation_module (micromolar)"
legend_constants[34] = "KD in component PKA_activation_module (micromolar)"
legend_constants[35] = "KPKI in component PKA_activation_module (micromolar)"
legend_algebraic[19] = "PKI in component PKA_activation_module (micromolar)"
legend_algebraic[3] = "fracPLBp in component phospholamban_regulation_module (dimensionless)"
legend_algebraic[20] = "fracPLB in component phospholamban_regulation_module (dimensionless)"
legend_constants[36] = "fracPLBo in component phospholamban_regulation_module (dimensionless)"
legend_states[15] = "PLBp in component phospholamban_regulation_module (micromolar)"
legend_algebraic[8] = "PLB in component phospholamban_regulation_module (micromolar)"
legend_states[16] = "Inhib1ptot in component phospholamban_regulation_module (micromolar)"
legend_states[17] = "Inhib1 in component phospholamban_regulation_module (micromolar)"
legend_states[18] = "Inhib1p in component phospholamban_regulation_module (micromolar)"
legend_algebraic[4] = "PP1_Inhib1p in component phospholamban_regulation_module (micromolar)"
legend_states[19] = "PP1 in component phospholamban_regulation_module (micromolar)"
legend_constants[37] = "PP1tot in component phospholamban_regulation_module (micromolar)"
legend_constants[38] = "PLBtot in component phospholamban_regulation_module (micromolar)"
legend_constants[39] = "Inhib1tot in component phospholamban_regulation_module (micromolar)"
legend_constants[40] = "epsilon in component phospholamban_regulation_module (dimensionless)"
legend_constants[41] = "kPKA_PLB in component phospholamban_regulation_module (first_order_rate_constant)"
legend_constants[42] = "KmPKA_PLB in component phospholamban_regulation_module (micromolar)"
legend_constants[43] = "kPKA_Inhib1 in component phospholamban_regulation_module (first_order_rate_constant)"
legend_constants[44] = "kPP1_PLB in component phospholamban_regulation_module (first_order_rate_constant)"
legend_constants[45] = "KmPP1_PLB in component phospholamban_regulation_module (micromolar)"
legend_constants[46] = "KmPKA_Inhib1 in component phospholamban_regulation_module (micromolar)"
legend_constants[47] = "VmaxPP2A_Inhib1 in component phospholamban_regulation_module (flux)"
legend_constants[48] = "KmPP2A_Inhib1 in component phospholamban_regulation_module (micromolar)"
legend_constants[49] = "KInhib1 in component phospholamban_regulation_module (micromolar)"
legend_rates[5] = "d/dt L in component beta_adrenergic_receptor_module (micromolar)"
legend_rates[1] = "d/dt Gs in component beta_adrenergic_receptor_module (micromolar)"
legend_rates[0] = "d/dt beta1_AR in component beta_adrenergic_receptor_module (micromolar)"
legend_rates[2] = "d/dt beta1_ARact in component beta_adrenergic_receptor_module (micromolar)"
legend_rates[3] = "d/dt beta1_AR_S464 in component beta_adrenergic_receptor_module (micromolar)"
legend_rates[4] = "d/dt beta1_AR_S301 in component beta_adrenergic_receptor_module (micromolar)"
legend_rates[7] = "d/dt Gs_alpha_GTPtot in component Gs_activation_module (micromolar)"
legend_rates[6] = "d/dt Gs_beta_gamma in component Gs_activation_module (micromolar)"
legend_rates[8] = "d/dt Gs_alpha_GDP in component Gs_activation_module (micromolar)"
legend_rates[11] = "d/dt Gs_alpha_GTP in component cyclic_AMP_metabolism_module (micromolar)"
legend_rates[14] = "d/dt fsk in component cyclic_AMP_metabolism_module (micromolar)"
legend_rates[10] = "d/dt AC in component cyclic_AMP_metabolism_module (micromolar)"
legend_rates[12] = "d/dt PDE in component cyclic_AMP_metabolism_module (micromolar)"
legend_rates[13] = "d/dt IBMX in component cyclic_AMP_metabolism_module (micromolar)"
legend_rates[9] = "d/dt cAMPtot in component cyclic_AMP_metabolism_module (micromolar)"
legend_rates[19] = "d/dt PP1 in component phospholamban_regulation_module (micromolar)"
legend_rates[18] = "d/dt Inhib1p in component phospholamban_regulation_module (micromolar)"
legend_rates[15] = "d/dt PLBp in component phospholamban_regulation_module (micromolar)"
legend_rates[16] = "d/dt Inhib1ptot in component phospholamban_regulation_module (micromolar)"
legend_rates[17] = "d/dt Inhib1 in component phospholamban_regulation_module (micromolar)"
return (legend_states, legend_algebraic, legend_voi, legend_constants)
def initConsts():
constants = [0.0] * sizeConstants; states = [0.0] * sizeStates;
states[0] = 0.0001
states[1] = 3.182
states[2] = 0.01205
states[3] = 0
states[4] = 1.154e-3
states[5] = 0.988
constants[0] = 1
constants[1] = 0.285
constants[2] = 0.062
constants[3] = 33
constants[4] = 3.83
constants[5] = 1.1e-3
constants[6] = 2.2e-3
constants[7] = 3.6e-3
constants[8] = 2.2e-3
states[6] = 0.02569
states[7] = 0.02505
states[8] = 6.44e-4
constants[9] = 16
constants[10] = 0.8
constants[11] = 1.21e3
states[9] = 0.8453
states[10] = 0.04706295
states[11] = 0.02241295
states[12] = 0.0389
constants[12] = 49.7e-3
constants[13] = 5e3
constants[14] = 38.9e-3
constants[15] = 0
states[13] = 0
constants[16] = 0
states[14] = 0
constants[17] = 0.2
constants[18] = 7.3
constants[19] = 5
constants[20] = 8.5
constants[21] = 1.03e3
constants[22] = 1.3
constants[23] = 860
constants[24] = 44
constants[25] = 0.4
constants[26] = 0.4
constants[27] = 30
constants[28] = 0.2e-3
constants[29] = 0.59
constants[30] = 0.025
constants[31] = 0.18
constants[32] = 9.14e-3
constants[33] = 1.64e-3
constants[34] = 4.375e-3
constants[35] = 2e-4
constants[36] = 0.9613
states[15] = 4.105
states[16] = 0.0526
states[17] = 0.2474
states[18] = 6.27339e-5
states[19] = 0.8374627
constants[37] = 0.89
constants[38] = 106
constants[39] = 0.3
constants[40] = 10
constants[41] = 54
constants[42] = 21
constants[43] = 60
constants[44] = 8.5
constants[45] = 7
constants[46] = 1
constants[47] = 14
constants[48] = 1
constants[49] = 1e-3
return (states, constants)
def computeRates(voi, states, constants):
rates = [0.0] * sizeStates; algebraic = [0.0] * sizeAlgebraic
rates[8] = constants[10]*states[7]-constants[11]*states[8]*states[6]
algebraic[2] = (states[11]*states[10])/constants[26]
rates[11] = states[7]-(algebraic[2]+states[11])
rates[10] = constants[12]-(algebraic[2]+states[10])
algebraic[1] = (states[12]*states[13])/constants[27]
rates[12] = constants[14]-(algebraic[1]+states[12])
rates[13] = constants[15]-(algebraic[1]+states[13])
algebraic[4] = (states[19]*states[18])/constants[49]
rates[19] = constants[37]-(algebraic[4]+states[19])
rates[18] = states[16]-(algebraic[4]+states[18])
algebraic[5] = (states[5]*states[0]*states[1])/(constants[1]*constants[2])
algebraic[0] = (states[5]*states[0])/constants[1]
rates[5] = constants[0]-(algebraic[0]+algebraic[5]+states[5])
rates[3] = constants[5]*(algebraic[0]+algebraic[5])-constants[6]*states[3]
algebraic[6] = (states[14]*states[10])/constants[24]
rates[14] = constants[16]-(algebraic[6]+states[14])
rates[9] = ((constants[17]*states[10]*constants[13])/(constants[21]+constants[13])+(constants[20]*algebraic[2]*constants[13])/(constants[25]+constants[13])+(constants[18]*algebraic[6]*constants[13])/(constants[23]+constants[13]))-(constants[19]*states[12]*states[9])/(constants[22]+states[9])
algebraic[7] = (states[0]*states[1])/constants[3]
rates[1] = constants[4]-(algebraic[7]+algebraic[5]+states[1])
rates[0] = states[2]-(algebraic[0]+algebraic[5]+algebraic[7]+states[0])
rates[7] = constants[9]*(algebraic[7]+algebraic[5])-constants[10]*states[7]
rates[6] = constants[9]*(algebraic[7]+algebraic[5])-constants[11]*states[8]*states[6]
rootfind_0(voi, constants, rates, states, algebraic)
rates[2] = (constants[6]*states[3]-constants[5]*(algebraic[0]+algebraic[5]))+(constants[8]*states[4]-constants[7]*algebraic[9]*states[2])
rates[4] = constants[7]*algebraic[9]*states[2]-constants[8]*states[4]
algebraic[8] = constants[38]-states[15]
rates[15] = (constants[41]*algebraic[9]*algebraic[8])/(constants[42]+algebraic[8])-(constants[44]*states[19]*states[15])/(constants[45]+states[15])
rates[16] = (constants[43]*algebraic[9]*states[16])/(constants[46]+states[16])-(constants[47]*states[17])/(constants[48]+states[17])
rates[17] = (constants[47]*states[17])/(constants[48]+states[17])-(constants[43]*algebraic[9]*states[16])/(constants[46]+states[16])
return(rates)
def computeAlgebraic(constants, states, voi):
algebraic = array([[0.0] * len(voi)] * sizeAlgebraic)
states = array(states)
voi = array(voi)
algebraic[2] = (states[11]*states[10])/constants[26]
algebraic[1] = (states[12]*states[13])/constants[27]
algebraic[4] = (states[19]*states[18])/constants[49]
algebraic[5] = (states[5]*states[0]*states[1])/(constants[1]*constants[2])
algebraic[0] = (states[5]*states[0])/constants[1]
algebraic[6] = (states[14]*states[10])/constants[24]
algebraic[7] = (states[0]*states[1])/constants[3]
algebraic[8] = constants[38]-states[15]
algebraic[3] = states[15]/constants[38]
algebraic[20] = algebraic[8]/constants[38]
return algebraic
initialGuess0 = None
def rootfind_0(voi, constants, rates, states, algebraic):
"""Calculate values of algebraic variables for DAE"""
from scipy.optimize import fsolve
global initialGuess0
if initialGuess0 is None: initialGuess0 = ones(11)*0.1
if not iterable(voi):
soln = fsolve(residualSN_0, initialGuess0, args=(algebraic, voi, constants, rates, states), xtol=1E-6)
initialGuess0 = soln
algebraic[9] = soln[0]
algebraic[10] = soln[1]
algebraic[11] = soln[2]
algebraic[12] = soln[3]
algebraic[13] = soln[4]
algebraic[14] = soln[5]
algebraic[15] = soln[6]
algebraic[16] = soln[7]
algebraic[17] = soln[8]
algebraic[18] = soln[9]
algebraic[19] = soln[10]
else:
for (i,t) in enumerate(voi):
soln = fsolve(residualSN_0, initialGuess0, args=(algebraic[:,i], voi[i], constants, rates[:i], states[:,i]), xtol=1E-6)
initialGuess0 = soln
algebraic[9][i] = soln[0]
algebraic[10][i] = soln[1]
algebraic[11][i] = soln[2]
algebraic[12][i] = soln[3]
algebraic[13][i] = soln[4]
algebraic[14][i] = soln[5]
algebraic[15][i] = soln[6]
algebraic[16][i] = soln[7]
algebraic[17][i] = soln[8]
algebraic[18][i] = soln[9]
algebraic[19][i] = soln[10]
def residualSN_0(algebraicCandidate, algebraic, voi, constants, rates, states):
resid = array([0.0] * 11)
algebraic[9] = algebraicCandidate[0]
algebraic[10] = algebraicCandidate[1]
algebraic[11] = algebraicCandidate[2]
algebraic[12] = algebraicCandidate[3]
algebraic[13] = algebraicCandidate[4]
algebraic[14] = algebraicCandidate[5]
algebraic[15] = algebraicCandidate[6]
algebraic[16] = algebraicCandidate[7]
algebraic[17] = algebraicCandidate[8]
algebraic[18] = algebraicCandidate[9]
algebraic[19] = algebraicCandidate[10]
resid[0] = (algebraic[19]-(constants[31]*constants[28])/(constants[28]+algebraic[9]+algebraic[11]))
resid[1] = (algebraic[14]-(algebraic[9]/constants[34])*algebraic[9]*(1.00000+algebraic[19]/constants[28]))
resid[2] = (algebraic[15]-algebraic[9]*(1.00000+algebraic[19]/constants[28]))
resid[3] = (algebraic[16]-(algebraic[11]/constants[34])*algebraic[11]*(1.00000+algebraic[19]/constants[28]))
resid[4] = (algebraic[17]-algebraic[11]*(1.00000+algebraic[19]/constants[28]))
resid[5] = (algebraic[13]-(constants[32]/algebraic[10])*algebraic[14])
resid[6] = (algebraic[18]-(constants[32]/algebraic[10])*algebraic[16])
resid[7] = (algebraic[10]-((states[9]-(algebraic[13]+2.00000*algebraic[14]+2.00000*algebraic[15]))-(algebraic[18]+2.00000*algebraic[16]+2.00000*algebraic[17])))
resid[8] = (algebraic[12]-((constants[32]*constants[33])/constants[34]+(constants[32]*algebraic[10])/constants[34]+(algebraic[10]*algebraic[10])/constants[34]))
resid[9] = (algebraic[9]-(2.00000*constants[29]*algebraic[10]*algebraic[10]-algebraic[9]*(1.00000+algebraic[19]/constants[28])*(algebraic[12]*algebraic[9]+algebraic[10]*algebraic[10])))
resid[10] = (algebraic[11]-(2.00000*constants[30]*algebraic[10]*algebraic[10]-algebraic[11]*(1.00000+algebraic[19]/constants[28])*(algebraic[12]*algebraic[11]+algebraic[10]*algebraic[10])))
return resid
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)