# 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)