################################################################################
#                                                                              #
#              STEPS model of Cerebellar Long-Term Depression                  #
#                                                                              #
#              Project supervision: Erik De Schutter (2009-2016)               #
#                                                                              #
#                           Scripts Authors:                                   #
#                                                                              #
#  Original development (2008-2012): Gabriela Antunes                          #
#  Modified for STEPS 1.3 (2013):    Iain Hepburn and Cory Simon               #
#  Addition of RKIP, improved calcium dynamics, and other                      #
#      modifications (2014-2016):    Anant Jain, Iain Hepburn, Himanshu Gangal #
#                                                                              #
#                                                                              #
################################################################################


# Example command-line usage:
# python LTD_STEPS_2.1.py ./data LTDsim 50 15 7 55 1 AMPAR Ca ERK MEK PKC

########################################################################     

import sys
import os

###################### Command line arguments ##########################

# The path to store data
path=sys.argv[1]

# ID of the data (e.g. 'LTDsim')
dataid=sys.argv[2]

# Ca pulse peak (uM)
ca_peak=sys.argv[3]

# Ca pulse width (s)
ca_width=sys.argv[4]

# Raf number
raf_number = sys.argv[5]

# PKC number
pkc_number = sys.argv[6]

# Index of the data (to be converted to int, e.g. '7'). 
# This will also control the random number seed
idx = sys.argv[7]

# Data to save. Species names to be given as a subset of:
# 'AMPAR', 'Ca', 'AA', 'PKC', 'ERK', 'MEK', 'PLA2', 'Raf', 'RKIP'
# Activites of those species will be recorded over time
data_record = sys.argv[8:]

if len(data_record) == 0: 
    print "WARNING: Will not record any data"

datapath = path+'/'+dataid+'_'+ca_peak+'_'+ca_width+'_'+raf_number+'_'+pkc_number

# Create storage locations if necessary.
try: os.mkdir(path)
except: pass
try: os.mkdir(datapath)
except: pass
try: os.mkdir(datapath +'/'+ str(idx))
except: pass


# Automatically plot some figures or not. Basically plots activity
# over time for any species recorded
plotfigs = True

########################################################################

import datetime
import math
import numpy
from numpy import mean
from pylab import *
import pickle

import steps.model as smodel
import steps.solver as swmdirect
import steps.geom as swm
import steps.rng as srng

# Simualation data recording time step (s)
DT = 0.05 

# Final time (s) simulated
INT = 3600

# Number of runs used to calculate the mean
# To see individual runs this should be set to 1
NITER = 1 

# Initial time that will be discharged (s). 
# The initial 10 minutes of all simulations should be discharged to allow the model to reach equilibrium.
tinit = 600 

# Avogadro constant
Na = 6.02214129e23

########################################################################

"""
# Ca2+ pulse. 
# The duration and amplitude of the pulse should be checked directly in the Ca2+ concentration
# The changes in [Ca2+] generated by the pulses are sensitive to stochasticity
# max_molar is the molar concentration per second injection (M/s)
# start defines the time the pulse starts (s)
# duration defines the duration of the pulse(s)

# start and duration are taken from command line arguments
"""

def square_pulse_ica(t,max_molar=float(ca_peak)*1e-6, duration=float(ca_width), start=1200.0):
    if t < start or t > start + duration: return 0
    else: return max_molar
    

######################################################################## 

def gen_model(): 
    mdl=smodel.Model() 
    
    Ca = smodel.Spec('Ca', mdl)
    PMCA = smodel.Spec('PMCA', mdl)
    Ca1PMCA = smodel.Spec('Ca1PMCA', mdl) 
    NCX = smodel.Spec('NCX', mdl) 
    Ca1NCX = smodel.Spec('Ca1NCX', mdl)  
    Ca2NCX = smodel.Spec('Ca2NCX', mdl)
    SERCA = smodel.Spec('SERCA', mdl) 
    Ca1SERCA = smodel.Spec('Ca1SERCA', mdl)  
    Ca2SERCA = smodel.Spec('Ca2SERCA', mdl)            
    PV = smodel.Spec('PV', mdl)
    MgPV = smodel.Spec('MgPV', mdl)  
    CaPV = smodel.Spec('CaPV', mdl)
    Mg2PV = smodel.Spec('Mg2PV', mdl)  
    Ca2PV = smodel.Spec('Ca2PV', mdl)
    CBf = smodel.Spec('CBf', mdl)              
    CaCBf = smodel.Spec('CaCBf', mdl) 
    Ca2CBf = smodel.Spec('Ca2CBf', mdl) 
    CBs = smodel.Spec('CBs', mdl)              
    CaCBs = smodel.Spec('CaCBs', mdl)  
    Ca2CBs = smodel.Spec('Ca2CBs', mdl) 
    PKC = smodel.Spec('PKC', mdl)
    Ca1PKC = smodel.Spec('Ca1PKC', mdl)
    Ca3PKC = smodel.Spec('Ca3PKC', mdl)
    AAPKC = smodel.Spec('AAPKC', mdl)    
    AACa1PKC = smodel.Spec('AACa1PKC', mdl)  
    AACa3PKC = smodel.Spec('AACa3PKC', mdl)                      
    PKCstar = smodel.Spec('PKCstar', mdl) 
    PKCstar2 = smodel.Spec('PKCstar2', mdl)       
    PKCstar4 = smodel.Spec('PKCstar4', mdl)
    PKCstar3 = smodel.Spec('PKCstar3', mdl) 
    Raf = smodel.Spec('Raf', mdl)  
    Rafstar = smodel.Spec('Rafstar', mdl)           
    PP5 = smodel.Spec('PP5', mdl)       
    PP5Rafstar = smodel.Spec('PP5Rafstar', mdl)        
    MEK = smodel.Spec('MEK', mdl)                    
    RafstarMEK = smodel.Spec('RafstarMEK', mdl)        
    MEKp = smodel.Spec('MEKp', mdl)
    PP2A = smodel.Spec('PP2A', mdl) 
    PP2AMEKp = smodel.Spec('PP2AMEKp', mdl)             
    RafstarMEKp = smodel.Spec('RafstarMEKp', mdl)        
    MEKstar = smodel.Spec('MEKstar', mdl)                
    PP2AMEKstar = smodel.Spec('PP2AMEKstar', mdl)        
    ERK = smodel.Spec('ERK', mdl)                         
    MEKstarERK = smodel.Spec('MEKstarERK', mdl)           
    ERKp = smodel.Spec('ERKp', mdl)                       
    MKP = smodel.Spec('MKP', mdl)                       
    MKPERKp = smodel.Spec('MKPERKp', mdl)             
    MEKstarERKp = smodel.Spec('MEKstarERKp', mdl)        
    ERKstar = smodel.Spec('ERKstar', mdl)             
    MKPERKstar = smodel.Spec('MKPERKstar', mdl) 
    PLA2 = smodel.Spec('PLA2', mdl)            
    Ca1PLA2 = smodel.Spec('Ca1PLA2', mdl)
    Ca2PLA2 = smodel.Spec('Ca2PLA2', mdl)
    PLA2memb = smodel.Spec('PLA2memb', mdl)
    Ca1PLA2memb = smodel.Spec('Ca1PLA2memb', mdl) 
    PLA2star1 = smodel.Spec('PLA2star1', mdl)                     
    ERKstarPLA2 = smodel.Spec('ERKstarPLA2', mdl)         
    PLA2star2 = smodel.Spec('PLA2star2', mdl)                              
    Ca1PLA2star2 = smodel.Spec('Ca1PLA2star2', mdl)
    Ca2PLA2star2 = smodel.Spec('Ca2PLA2star2', mdl)
    PLA2star2memb = smodel.Spec('PLA2star2memb', mdl)
    Ca1PLA2star2memb = smodel.Spec('Ca1PLA2star2memb', mdl)
    Ca2PLA2star2memb = smodel.Spec('Ca2PLA2star2memb', mdl)
    PLA2star1APC = smodel.Spec('PLA2star1APC', mdl)
    ERKstarCa1PLA2 = smodel.Spec('ERKstarCa1PLA2', mdl)
    ERKstarCa2PLA2 = smodel.Spec('ERKstarCa2PLA2', mdl)
    PP2ACa1PLA2star2 = smodel.Spec('PP2ACa1PLA2star2', mdl)
    PP2ACa2PLA2star2 = smodel.Spec('PP2ACa2PLA2star2', mdl)
    PP1Ca1PLA2star2 = smodel.Spec('PP1Ca1PLA2star2', mdl)
    PP1Ca2PLA2star2 = smodel.Spec('PP1Ca2PLA2star2', mdl)
    PLA2star2membAPC = smodel.Spec('PLA2star2membAPC', mdl)
    Ca1PLA2star2membAPC = smodel.Spec('Ca1PLA2star2membAPC', mdl)
    Ca2PLA2star2membAPC = smodel.Spec('Ca2PLA2star2membAPC', mdl)
    PP1 = smodel.Spec('PP1', mdl)
    PP1PLA2star2 = smodel.Spec('PP1PLA2star2', mdl) 
    PP2APLA2star2 = smodel.Spec('PP2APLA2star2', mdl)           
    AA = smodel.Spec('AA', mdl)              
    AMPAR = smodel.Spec('AMPAR', mdl)    
    PKCstarAMPAR = smodel.Spec('PKCstarAMPAR', mdl) 
    PKCstar2AMPAR = smodel.Spec('PKCstar2AMPAR', mdl)       
    PKCstar4AMPAR = smodel.Spec('PKCstar4AMPAR', mdl)
    PKCstar3AMPAR = smodel.Spec('PKCstar3AMPAR', mdl) 
    PKCstarGRIPAMPAR = smodel.Spec('PKCstarGRIPAMPAR', mdl) 
    PKCstar2GRIPAMPAR = smodel.Spec('PKCstar2GRIPAMPAR', mdl)       
    PKCstar4GRIPAMPAR = smodel.Spec('PKCstar4GRIPAMPAR', mdl)
    PKCstar3GRIPAMPAR = smodel.Spec('PKCstar3GRIPAMPAR', mdl)   
    AMPAR_P = smodel.Spec('AMPAR_P', mdl)                  
    PP2AAMPAR_P = smodel.Spec('PP2AAMPAR_P', mdl)   
    PP2AGRIPAMPAR_P = smodel.Spec('PP2AGRIPAMPAR_P', mdl)   
    AMPARextra = smodel.Spec('AMPARextra', mdl)  
    AMPARextra_P = smodel.Spec('AMPARextra_P', mdl)
    PP2AAMPARextra_P = smodel.Spec('PP2AAMPARextra_P', mdl)
    GRIP = smodel.Spec('GRIP', mdl)
    GRIPAMPAR = smodel.Spec('GRIPAMPAR', mdl)
    GRIPAMPAR_P = smodel.Spec('GRIPAMPAR_P', mdl)   
    AMPARdend = smodel.Spec('AMPARdend', mdl)  
    AMPARdend_P = smodel.Spec('AMPARdend_P', mdl)  
    PP2AAMPARdend_P = smodel.Spec('PP2AAMPARdend_P', mdl)  
    AMPARcyt = smodel.Spec('AMPARcyt', mdl)  
    AMPARcyt_P = smodel.Spec('AMPARcyt_P', mdl)  
    PP2AAMPARcyt_P = smodel.Spec('PP2AAMPARcyt_P', mdl)
    RKIP = smodel.Spec('RKIP', mdl)
    RafRKIP = smodel.Spec('RafRKIP', mdl)
    RKIPstar = smodel.Spec('RKIPstar', mdl)
    RP = smodel.Spec('RP', mdl)
    RKIPstarRP = smodel.Spec ('RKIPstarRP', mdl)
    RafRKIPPKCstar = smodel.Spec('RafRKIPPKCstar', mdl)
    RafRKIPPKCstar2= smodel.Spec ('RafRKIPPKCstar2', mdl)
    RafRKIPPKCstar3= smodel.Spec ('RafRKIPPKCstar3', mdl)
    RafRKIPPKCstar4= smodel.Spec ('RafRKIPPKCstar4', mdl)
    MEKRKIP= smodel.Spec('MEKRKIP', mdl)
    MEKRKIPPKCstar = smodel.Spec('MEKRKIPPKCstar', mdl)
    MEKRKIPPKCstar2= smodel.Spec ('MEKRKIPPKCstar2', mdl)
    MEKRKIPPKCstar3= smodel.Spec ('MEKRKIPPKCstar3', mdl)
    MEKRKIPPKCstar4= smodel.Spec ('MEKRKIPPKCstar4', mdl)
    
    extra = smodel.Volsys('extra', mdl)
    vsys = smodel.Volsys('vsys', mdl)
    s = smodel.Surfsys('memb', mdl)
    ERs = smodel.Surfsys('ERmemb', mdl)
    cytER = smodel.Volsys('cytER', mdl)

    ########################################################################
    ########################################################################

    # Reactions
    
    # Ca influx (Ca = calcium)
    CainfluxR = smodel.Reac('Cainflux', vsys, rhs=[Ca])
    # reaction constant for zero-order reaction units M/s in STEPS 1.2 and above
    CainfluxR.kcst = 0.0 
    
    # Ca + PMCA <->  Ca1PMCA -> PMCA
    Reac1 = smodel.SReac('pump2_f', s, ilhs = [Ca], slhs =  [PMCA], srhs = [Ca1PMCA])
    Reac1.kcst = 2500e6
    Reac2 = smodel.SReac('Reac2', s, slhs = [Ca1PMCA], irhs = [Ca], srhs = [PMCA])
    Reac2.kcst = 200 
    Reac3 = smodel.SReac('Reac3', s, slhs = [Ca1PMCA], srhs = [PMCA])	
    Reac3.kcst = 50
    
    # Ca + NCX <->  Ca1NCX + Ca <->  Ca2NCX -> NCX 
    Reac4 = smodel.SReac('Reac4', s, ilhs = [Ca,Ca], slhs = [NCX], srhs = [Ca2NCX])
    Reac4.kcst = 93.827e12  
    Reac5 = smodel.SReac('Reac5', s, slhs = [Ca2NCX], irhs = [Ca,Ca], srhs = [NCX])
    Reac5.kcst = 4000.0
    Reac6 = smodel.SReac('Reac6', s, slhs = [Ca2NCX], srhs = [NCX])
    Reac6.kcst = 1000
    
    # Ca + SERCA <->  Ca1SERCA +Ca <->  Ca2SERCA  ->  SERCA
    Reac7 = smodel.SReac('Reac7', ERs, olhs = [Ca,Ca], slhs = [SERCA], srhs = [Ca2SERCA])
    Reac7.kcst = 3428.7485714376226e12 
    Reac8 = smodel.SReac('Reac8', ERs, slhs = [Ca2SERCA], orhs = [Ca,Ca], srhs = [SERCA])
    Reac8.kcst = 199.95577085780272
    Reac9 = smodel.SReac('Reac9', ERs, slhs = [Ca2SERCA], srhs = [SERCA], irhs = [Ca,Ca])
    Reac9.kcst = 50
    
    # Leak
    Reac10 = smodel.Reac('Reac10', vsys, rhs = [Ca])	
    # reaction constant for zero-order reaction units M/s in STEPS 1.2 and above
    Reac10.kcst = 150/(Na*0.08e-15)
    
    
    
    # PV + Ca <-> CaPV + Ca <-> Ca2PV
    Reac15 = smodel.Reac('Reac15', vsys, lhs = [PV, Ca], rhs = [CaPV])
    Reac15.kcst = 107e6 
    Reac16 = smodel.Reac('Reac16', vsys, lhs = [CaPV], rhs = [PV, Ca])
    Reac16.kcst = 0.95
    Reac17 = smodel.Reac('Reac17', vsys, lhs = [CaPV, Ca], rhs = [Ca2PV])
    Reac17.kcst = 107e6 
    Reac18 = smodel.Reac('Reac18', vsys, lhs = [Ca2PV], rhs = [Ca, CaPV])
    Reac18.kcst = 0.95
    
    
    
    # PV + Mg <-> MgPV <-> Mg2PV, apparent rate constants, [Mg] = 590uM 
    Reac19 = smodel.Reac('Reac19', vsys, lhs = [PV], rhs = [MgPV])
    Reac19.kcst = 472
    Reac20 = smodel.Reac('Reac20', vsys, lhs = [MgPV], rhs = [PV])
    Reac20.kcst = 25
    Reac21 = smodel.Reac('Reac21', vsys, lhs = [MgPV], rhs = [Mg2PV])
    Reac21.kcst = 472
    Reac22 = smodel.Reac('Reac22', vsys, lhs = [Mg2PV], rhs = [MgPV])
    Reac22.kcst = 25
    
    
    
    # CBs + Ca <-> CaCBs + Ca <-> Ca2CBs, CBs - high affinity site
    Reac23 = smodel.Reac('Reac23', vsys, lhs = [CBs, Ca], rhs = [CaCBs])
    Reac23.kcst = 5.5e6
    Reac24 = smodel.Reac('Reac24', vsys, lhs = [CaCBs], rhs = [CBs, Ca])
    Reac24.kcst = 2.6
    Reac25 = smodel.Reac('Reac25', vsys, lhs = [CaCBs, Ca], rhs = [Ca2CBs])
    Reac25.kcst = 5.5e6 
    Reac26 = smodel.Reac('Reac26', vsys, lhs = [Ca2CBs], rhs = [CaCBs, Ca])
    Reac26.kcst = 2.6
    
    
    
    # CBf + Ca <-> CaCBf + Ca <-> Ca2CBf, CBf - medium affinity site
    Reac27 = smodel.Reac('Reac27', vsys, lhs = [CBf, Ca], rhs = [CaCBf])
    Reac27.kcst = 43.5e6 
    Reac28 = smodel.Reac('Reac28', vsys, lhs = [CaCBf], rhs = [CBf, Ca])
    Reac28.kcst = 35.8
    Reac29 = smodel.Reac('Reac29', vsys, lhs = [CaCBf, Ca], rhs = [Ca2CBf])
    Reac29.kcst = 43.5e6 
    Reac30 = smodel.Reac('Reac30', vsys, lhs = [Ca2CBf], rhs = [CaCBf, Ca])
    Reac30.kcst = 35.8
    
    
    
    # PKC + Ca <-> CaPKC + 2Ca <-> Ca3PKC <-> Ca3PKC* (PKCstar)
    Reac31 = smodel.Reac('Reac31', vsys, lhs = [PKC, Ca], rhs = [Ca1PKC])
    Reac31.kcst = 13.3e6
    Reac32 = smodel.Reac('Reac32', vsys, lhs = [Ca1PKC], rhs = [PKC, Ca])
    Reac32.kcst = 12
    Reac33 = smodel.Reac('Reac33', vsys, lhs = [Ca1PKC, Ca, Ca], rhs = [Ca3PKC])
    Reac33.kcst = 1.e12 
    Reac34 = smodel.Reac('Reac34', vsys, lhs = [Ca3PKC], rhs = [Ca1PKC, Ca, Ca])
    Reac34.kcst = 12
    Reac35 = smodel.SReac('Reac35', s, ilhs = [Ca3PKC], srhs = [PKCstar])
    Reac35.kcst = 11.3  
    Reac36 = smodel.SReac('Reac36', s, slhs = [PKCstar], irhs = [Ca3PKC])
    Reac36.kcst = 0.23 
    
    
    
    # PKC + AA <-> AAPKC <-> AAPKC* (PKCstar2)
    Reac37 = smodel.Reac('Reac37', vsys, lhs = [PKC, AA], rhs = [AAPKC])
    Reac37.kcst = 1e6 
    Reac38 = smodel.Reac('Reac38', vsys, lhs = [AAPKC], rhs = [PKC, AA])
    Reac38.kcst = 10
    Reac39 = smodel.SReac('Reac39', s, ilhs = [AAPKC], srhs = [PKCstar2])
    Reac39.kcst = 0.017 
    Reac40 = smodel.SReac('Reac40', s, slhs = [PKCstar2], irhs = [AAPKC])
    Reac40.kcst = 0.0055
    
    
    
    # Ca1PKC + AA <-> AACa1PKC <-> AACa1PKC* (PKCstar3) + 2Ca <-> AACa3PKC* (PKCstar4)
    Reac41 = smodel.Reac('Reac41', vsys, lhs = [Ca1PKC, AA], rhs = [AACa1PKC])
    Reac41.kcst = 1e6 
    Reac42 = smodel.Reac('Reac42', vsys, lhs = [AACa1PKC], rhs = [Ca1PKC, AA])
    Reac42.kcst = 10
    Reac43 = smodel.SReac('Reac43', s, ilhs = [AACa1PKC], srhs = [PKCstar3])
    Reac43.kcst = 0.017 
    Reac44 = smodel.SReac('Reac44', s, slhs = [PKCstar3], irhs = [AACa1PKC])
    Reac44.kcst = 0.0055
    Reac45 = smodel.SReac('Reac45', s, slhs = [PKCstar3], ilhs = [Ca, Ca], srhs = [PKCstar4])
    Reac45.kcst =  1.0e12
    Reac46 = smodel.SReac('Reac46', s, slhs = [PKCstar4], srhs = [PKCstar3], irhs = [Ca, Ca])
    Reac46.kcst = 12
    
    
    
    # AAPKC + Ca <-> AACa1PKC <-> AACa1PKC + 2Ca <-> AACa3PKC  <-> AACa3PKC* (PKCstar4)
    Reac47 = smodel.Reac('Reac47', vsys, lhs = [AAPKC, Ca], rhs = [AACa1PKC])
    Reac47.kcst = 13.3e6
    Reac48 = smodel.Reac('Reac48', vsys, lhs = [AACa1PKC], rhs = [AAPKC, Ca])
    Reac48.kcst = 12
    Reac49 = smodel.Reac('Reac49', vsys, lhs = [AACa1PKC, Ca, Ca], rhs = [AACa3PKC])
    Reac49.kcst = 1.0e12 
    Reac50 = smodel.Reac('Reac50', vsys, lhs = [AACa3PKC], rhs = [AACa1PKC, Ca, Ca])
    Reac50.kcst = 12
    Reac51 = smodel.SReac('Reac51', s, ilhs = [AACa3PKC], srhs = [PKCstar4])
    Reac51.kcst = 11.3 
    Reac52 = smodel.SReac('Reac52', s, slhs = [PKCstar4], irhs = [AACa3PKC])
    Reac52.kcst = 0.23 
    
    
    
    # Ca3PKC + AA <-> AACa3PKC
    Reac53 = smodel.Reac('Reac53', vsys, lhs = [Ca3PKC, AA], rhs = [AACa3PKC])
    Reac53.kcst = 1e6 
    Reac54 = smodel.Reac('Reac54', vsys, lhs = [AACa3PKC], rhs = [Ca3PKC, AA])
    Reac54.kcst = 10
    
    
    
    # AAPKC* + Ca <-> AACa1PKC*
    Reac55 = smodel.SReac('Reac55', s, slhs = [PKCstar2], ilhs = [Ca], srhs = [PKCstar3])
    Reac55.kcst =  13.3e6
    Reac56 = smodel.SReac('Reac56', s, slhs = [PKCstar3], srhs = [PKCstar2], irhs = [Ca])
    Reac56.kcst = 12
    
          
    
    # Ca3PKC* + AA <-> AACa3PKC*
    Reac57 = smodel.SReac('Reac57', s, ilhs = [AA], slhs = [PKCstar], srhs = [PKCstar4])
    Reac57.kcst = 1e6 
    Reac58 = smodel.SReac('Reac58', s, slhs = [PKCstar4], irhs = [AA], srhs = [PKCstar])
    Reac58.kcst = 10
        
    
    
    # Ca3PKC* + RafRKIP <-> Ca3PKC*.RafRKIP -> Ca3PKC* + RafRKIP
    Reac59 = smodel.SReac('Reac59', s, ilhs = [RafRKIP], slhs = [PKCstar], srhs = [RafRKIPPKCstar])
    Reac59.kcst = 0.625e6
    Reac60 = smodel.SReac('Reac60', s, slhs = [RafRKIPPKCstar], irhs = [RafRKIP], srhs = [PKCstar])
    Reac60.kcst = 0.00245
    Reac61 = smodel.SReac('Reac61', s, slhs = [RafRKIPPKCstar], srhs = [PKCstar], irhs = [RKIPstar, Rafstar])
    Reac61.kcst = 0.0315
    
  
    
    # AAPKC* + RafRKIP <-> AAPKC*.RafRKIP -> AAPKC* + RafRKIP
    Reac62 = smodel.SReac('Reac62', s, ilhs = [RafRKIP], slhs = [PKCstar2], srhs = [RafRKIPPKCstar2])
    Reac62.kcst =  0.625e6  
    Reac63 = smodel.SReac('Reac63', s, slhs = [RafRKIPPKCstar2], irhs = [RafRKIP], srhs = [PKCstar2])
    Reac63.kcst = 0.00245
    Reac64 = smodel.SReac('Reac64', s, slhs = [RafRKIPPKCstar2], srhs = [PKCstar2], irhs = [RKIPstar, Rafstar])
    Reac64.kcst =  0.0315
    
        
    
    # AACa1PKC* + RafRKIP -> AACa1PKC*.RafRKIP -> AACa1PKC* + RafRKIP 
    Reac65 = smodel.SReac('Reac65', s, ilhs = [RafRKIP], slhs = [PKCstar3], srhs = [RafRKIPPKCstar3])
    Reac65.kcst =  0.625e6  
    Reac66 = smodel.SReac('Reac66', s, slhs = [RafRKIPPKCstar3], irhs = [RafRKIP], srhs = [PKCstar3])
    Reac66.kcst = 0.00245
    Reac67 = smodel.SReac('Reac67', s, slhs = [RafRKIPPKCstar3], srhs = [PKCstar3], irhs = [RKIPstar, Rafstar])
    Reac67.kcst =  0.0315

        
    
    # AACa3PKC* + RafRKIP <-> AACa3PKC*.RafRKIP -> AACa3PKC* + RafRKIP
    Reac68 = smodel.SReac('Reac68', s, ilhs = [RafRKIP], slhs = [PKCstar4], srhs = [RafRKIPPKCstar4])
    Reac68.kcst =  0.625e6  
    Reac69 = smodel.SReac('Reac69', s, slhs = [RafRKIPPKCstar4], irhs = [RafRKIP], srhs = [PKCstar4])
    Reac69.kcst = 0.00245
    Reac70 = smodel.SReac('Reac70', s, slhs = [RafRKIPPKCstar4], srhs = [PKCstar4], irhs = [RKIPstar, Rafstar])
    Reac70.kcst = 0.0315
    
        
    
    # Raf + RKIP --k1-> RafRKIP --k2-> Raf + RKIP
    Reac71 = smodel.Reac('Reac71',vsys, lhs = [Raf , RKIP], rhs = [RafRKIP])
    Reac71.kcst = 0.53e6
    Reac72 = smodel.Reac('Reac72', vsys, lhs = [RafRKIP], rhs =[Raf , RKIP])
    Reac72.kcst = 0.072
    
    
    
    # Ca3PKC* + MEKRKIP <-> Ca3PKC*.MEKRKIP -> Ca3PKC* + MEKRKIP
    Reac73 = smodel.SReac('Reac73', s, ilhs = [MEKRKIP], slhs = [PKCstar], srhs = [MEKRKIPPKCstar])
    Reac73.kcst = 0.625e6
    Reac74 = smodel.SReac('Reac74', s, slhs = [MEKRKIPPKCstar], irhs = [MEKRKIP], srhs = [PKCstar])
    Reac74.kcst = 0.00245
    Reac75 = smodel.SReac('Reac75', s, slhs = [MEKRKIPPKCstar], srhs = [PKCstar], irhs = [RKIPstar, MEKp])
    Reac75.kcst = 0.0315
    
    
    
    # AAPKC* + MEKRKIP <-> AAPKC*.MEKRKIP -> AAPKC* + MEKRKIP
    Reac76 = smodel.SReac('Reac76', s, ilhs = [MEKRKIP], slhs = [PKCstar2], srhs = [MEKRKIPPKCstar2])
    Reac76.kcst =  0.625e6  
    Reac77 = smodel.SReac('Reac77', s, slhs = [MEKRKIPPKCstar2], irhs = [MEKRKIP], srhs = [PKCstar2])
    Reac77.kcst = 0.00245
    Reac78 = smodel.SReac('Reac78', s, slhs = [MEKRKIPPKCstar2], srhs = [PKCstar2], irhs = [MEKp, RKIPstar])
    Reac78.kcst =  0.0315
    
        
    
    # AACa1PKC* + MEKRKIP -> AACa1PKC*.MEKRKIP -> AACa1PKC* + MEKRKIP 
    Reac79 = smodel.SReac('Reac79', s, ilhs = [MEKRKIP], slhs = [PKCstar3], srhs = [MEKRKIPPKCstar3])
    Reac79.kcst =  0.625e6  
    Reac80 = smodel.SReac('Reac80', s, slhs = [MEKRKIPPKCstar3], irhs = [MEKRKIP], srhs = [PKCstar3])
    Reac80.kcst = 0.00245
    Reac81 = smodel.SReac('Reac81', s, slhs = [MEKRKIPPKCstar3], srhs = [PKCstar3], irhs = [MEKp, RKIPstar])
    Reac81.kcst =  0.0315

    
    
    # AACa3PKC* + MEKRKIP <-> AACa3PKC*.MEKRKIP -> AACa3PKC* + MEKRKIP
    Reac82 = smodel.SReac('Reac82', s, ilhs = [MEKRKIP], slhs = [PKCstar4], srhs = [MEKRKIPPKCstar4])
    Reac82.kcst =  0.625e6  
    Reac83 = smodel.SReac('Reac83', s, slhs = [MEKRKIPPKCstar4], irhs = [MEKRKIP], srhs = [PKCstar4])
    Reac83.kcst = 0.00245
    Reac84 = smodel.SReac('Reac84', s, slhs = [MEKRKIPPKCstar4], srhs = [PKCstar4], irhs = [RKIPstar, MEKp])
    Reac84.kcst = 0.0315
    
    
    
    # MEK + RKIP --k1-> MEKRKIP --k2-> MEK + RKIP
    Reac85 = smodel.Reac('Reac85',vsys, lhs = [MEK , RKIP], rhs = [MEKRKIP])
    Reac85.kcst = 0.53e6
    Reac86 = smodel.Reac('Reac86', vsys, lhs = [MEKRKIP], rhs =[MEK , RKIP])
    Reac86.kcst = 0.072
    
    
    
    # RKIPstar + RP --k9-> RKIPstarRP --k10-> RKIPstar + RP
    Reac87 = smodel.Reac('Reac87' , vsys, lhs =[RKIPstar, RP], rhs =[RKIPstarRP])
    Reac87.kcst = 0.92e6
    Reac88 = smodel.Reac('Reac88', vsys, lhs =[RKIPstarRP], rhs=[RKIPstar, RP])
    Reac88.kcst = 0.00122
    # RKIPstarRP --k11-> RKIP + RP
    Reac89 = smodel.Reac('Reac89' , vsys, lhs =[RKIPstarRP], rhs = [RKIP, RP])
    Reac89.kcst = 0.87
    
    
    
    # PP5 (PP) + Raf* <-> PP5.Raf* -> PP5 + Raf
    Reac90 = smodel.Reac('Reac90', vsys, lhs = [PP5, Rafstar], rhs = [PP5Rafstar])
    Reac90.kcst = .55e6 
    Reac91 = smodel.Reac('Reac91', vsys, lhs = [PP5Rafstar], rhs = [PP5, Rafstar])
    Reac91.kcst = 2
    Reac92 = smodel.Reac('Reac92', vsys, lhs = [PP5Rafstar], rhs = [PP5, Raf])
    Reac92.kcst = 0.5
    
    
    
    # Raf* + MEK <-> Raf*.MEK -> Raf* + MEKp <->Raf*.MEKp -> Raf* + MEK* (MEKstar)
    Reac93 = smodel.Reac('Reac93', vsys, lhs = [Rafstar, MEK], rhs = [RafstarMEK])
    Reac93.kcst = 0.65e6 
    Reac94 = smodel.Reac('Reac94', vsys, lhs = [RafstarMEK], rhs = [Rafstar, MEK])
    Reac94.kcst = 0.065
    Reac95 = smodel.Reac('Reac95', vsys, lhs = [RafstarMEK], rhs = [Rafstar, MEKp])
    Reac95.kcst = 1.0
    Reac96 = smodel.Reac('Reac96', vsys, lhs = [Rafstar, MEKp], rhs = [RafstarMEKp])
    Reac96.kcst = 0.65e6 
    Reac97 = smodel.Reac('Reac97', vsys, lhs = [RafstarMEKp], rhs = [Rafstar, MEKp])
    Reac97.kcst = 0.065
    Reac98 = smodel.Reac('Reac98', vsys, lhs = [RafstarMEKp], rhs = [Rafstar, MEKstar])
    Reac98.kcst = 1.0
    
    
    
    # PP2A + MEK* <-> PP2A.MEK* -> PP2A + MEKp <-> PP2A.MEKp -> PP2A + MEK
    Reac99 = smodel.Reac('Reac99', vsys, lhs = [PP2A, MEKstar], rhs = [PP2AMEKstar])
    Reac99.kcst = 0.75e6 
    Reac100 = smodel.Reac('Reac100', vsys, lhs = [PP2AMEKstar], rhs = [PP2A, MEKstar])
    Reac100.kcst = 2
    Reac101 = smodel.Reac('Reac101', vsys, lhs = [PP2AMEKstar], rhs = [PP2A, MEKp])
    Reac101.kcst = 0.5
    Reac102 = smodel.Reac('Reac102', vsys, lhs = [PP2A, MEKp], rhs = [PP2AMEKp])
    Reac102.kcst = 0.75e6 
    Reac103 = smodel.Reac('Reac103', vsys, lhs = [PP2AMEKp], rhs = [PP2A, MEKp])
    Reac103.kcst = 2
    Reac104 = smodel.Reac('Reac104', vsys, lhs = [PP2AMEKp], rhs = [PP2A, MEK])
    Reac104.kcst = 0.5
    
    
    
    # MEK* + ERK <-> MEK*.ERK -> MEK* + ERKp <-> MEK*.ERKp -> MEK* + ERK* (ERKstar)
    Reac105 = smodel.Reac('Reac105', vsys, lhs = [MEKstar, ERK], rhs = [MEKstarERK])
    Reac105.kcst = 16.2e6 
    Reac106 = smodel.Reac('Reac106', vsys, lhs = [MEKstarERK], rhs = [MEKstar, ERK])
    Reac106.kcst = 0.6
    Reac107 = smodel.Reac('Reac107', vsys, lhs = [MEKstarERK], rhs = [MEKstar, ERKp])
    Reac107.kcst = 0.15
    Reac108 = smodel.Reac('Reac108', vsys, lhs = [MEKstar, ERKp], rhs = [MEKstarERKp])
    Reac108.kcst = 16.2e6 
    Reac109 = smodel.Reac('Reac109', vsys, lhs = [MEKstarERKp], rhs = [MEKstar, ERKp])
    Reac109.kcst = 0.6
    Reac110 = smodel.Reac('Reac110', vsys, lhs = [MEKstarERKp], rhs = [MEKstar, ERKstar])
    Reac110.kcst = 0.3 
    
    
    
    # MKP + ERK* <-> MKP.ERK* -> MKP + ERKp <-> MKP.ERKp -> MKP + ERK
    Reac111 = smodel.Reac('Reac111', vsys, lhs = [MKP, ERKstar], rhs = [MKPERKstar])
    Reac111.kcst = 13e6 
    Reac112 = smodel.Reac('Reac112', vsys, lhs = [MKPERKstar], rhs = [MKP, ERKstar])
    Reac112.kcst = 0.396
    Reac113 = smodel.Reac('Reac113', vsys, lhs = [MKPERKstar], rhs = [MKP, ERKp])
    Reac113.kcst = 0.099 
    Reac114 = smodel.Reac('Reac114', vsys, lhs = [MKP, ERKp], rhs = [MKPERKp])
    Reac114.kcst = 28e6 
    Reac115 = smodel.Reac('Reac115', vsys, lhs = [MKPERKp], rhs = [MKP, ERKp])
    Reac115.kcst = 0.56
    Reac116 = smodel.Reac('Reac116', vsys, lhs = [MKPERKp], rhs = [MKP, ERK])
    Reac116.kcst = 0.14 
    
    
    
    # Ca + PLA2 <-> Ca1PLA2 + Ca <-> Ca2PLA2 <-> Ca2PLA2* (PLA2star1)
    Reac117 = smodel.Reac('Reac117', vsys, lhs = [PLA2, Ca], rhs = [Ca1PLA2])
    Reac117.kcst = 1.93e6
    Reac118 = smodel.Reac('Reac118', vsys, lhs = [Ca1PLA2], rhs = [PLA2, Ca])
    Reac118.kcst = 108
    Reac119 = smodel.Reac('Reac119', vsys, lhs = [Ca, Ca1PLA2], rhs = [Ca2PLA2])
    Reac119.kcst = 10.8e6 
    Reac120 = smodel.Reac('Reac120', vsys, lhs = [Ca2PLA2], rhs = [Ca1PLA2, Ca])
    Reac120.kcst = 108
    Reac121 = smodel.SReac('Reac121', s, ilhs = [Ca2PLA2], srhs = [PLA2star1])
    Reac121.kcst = 300 
    Reac122 = smodel.SReac('Reac122', s, slhs = [PLA2star1], irhs = [Ca2PLA2])
    Reac122.kcst = 15
    
    
    
    # Ca1PLA2 <-> Ca1PLA2memb
    Reac123 = smodel.SReac('Reac123', s, ilhs = [Ca1PLA2],  srhs = [Ca1PLA2memb])
    Reac123.kcst = 30
    Reac124 = smodel.SReac('Reac124', s, slhs = [Ca1PLA2memb], irhs = [Ca1PLA2])
    Reac124.kcst = 15
    
    
    
    # PLA2 <-> PLA2memb 
    Reac125 = smodel.SReac('Reac125', s, ilhs = [PLA2], srhs = [PLA2memb])
    Reac125.kcst = 3
    Reac126 = smodel.SReac('Reac126', s, slhs = [PLA2memb], irhs = [PLA2])
    Reac126.kcst = 15
    
    
    
    # PLA2memb + Ca <-> Ca1PLA2memb + Ca <-> Ca2PLA2* (PLA2star1)
    Reac127 = smodel.SReac('Reac127', s, slhs = [PLA2memb], ilhs = [Ca], srhs = [Ca1PLA2memb])
    Reac127.kcst = 1.93e6
    Reac128 = smodel.SReac('Reac128', s, ilhs = [Ca1PLA2memb], srhs = [PLA2memb], irhs = [Ca])
    Reac128.kcst = 0.41
    Reac129 = smodel.SReac('Reac129', s, slhs = [Ca1PLA2memb], ilhs = [Ca], srhs = [PLA2star1])
    Reac129.kcst = 10.8e6 
    Reac130 = smodel.SReac('Reac130', s, ilhs = [PLA2star1], srhs = [Ca1PLA2memb], irhs = [Ca])
    Reac130.kcst = 2.5
    
    
    
    # PLA2star1 <-> PLA2star1APC -> PLA2star1 + AA
    Reac131 = smodel.SReac('Reac131', s, slhs = [PLA2star1],  srhs = [PLA2star1APC])
    Reac131.kcst = 43  
    Reac132 = smodel.SReac('Reac132', s, slhs = [PLA2star1APC],  srhs = [PLA2star1])
    Reac132.kcst = 600
    Reac133 = smodel.SReac('Reac133', s, slhs = [PLA2star1APC], irhs = [AA], srhs = [PLA2star1])
    Reac133.kcst = 450
    
    
    
    # ERK* + PLA2 <-> ERK*.PLA2 -> ERK* + PLA2p (PLA2star2)
    Reac134 = smodel.Reac('Reac134', vsys, lhs = [ERKstar, PLA2], rhs = [ERKstarPLA2])
    Reac134.kcst = 4e6 
    Reac135 = smodel.Reac('Reac135', vsys, lhs = [ERKstarPLA2], rhs = [ERKstar, PLA2])
    Reac135.kcst = 1
    Reac136 = smodel.Reac('Reac136', vsys, lhs = [ERKstarPLA2], rhs = [ERKstar, PLA2star2])
    Reac136.kcst = 14
    
    
    
    # ERK* + Ca1PLA2 <-> ERK*.Ca1PLA2 -> ERK* + Ca1PLA2p (Ca1PLA2star2)
    Reac137 = smodel.Reac('Reac137', vsys, lhs = [ERKstar, Ca1PLA2], rhs = [ERKstarCa1PLA2])
    Reac137.kcst = 4e6 
    Reac138 = smodel.Reac('Reac138', vsys, lhs = [ERKstarCa1PLA2], rhs = [ERKstar, Ca1PLA2])
    Reac138.kcst = 1
    Reac139 = smodel.Reac('Reac139', vsys, lhs = [ERKstarCa1PLA2], rhs = [ERKstar, Ca1PLA2star2])
    Reac139.kcst = 14
    
    
    
    # ERK* + Ca2PLA2 <-> ERK*.Ca2PLA2 -> ERK* + Ca2PLA2p (Ca2PLA2star2)
    Reac140 = smodel.Reac('Reac140', vsys, lhs = [ERKstar, Ca2PLA2], rhs = [ERKstarCa2PLA2])
    Reac140.kcst = 4e6 
    Reac141 = smodel.Reac('Reac141', vsys, lhs = [ERKstarCa2PLA2], rhs = [ERKstar, Ca2PLA2])
    Reac141.kcst = 1
    Reac142 = smodel.Reac('Reac142', vsys, lhs = [ERKstarCa2PLA2], rhs = [ERKstar, Ca2PLA2star2])
    Reac142.kcst = 14
    
    
    # PP2A + PLA2p <-> PP2A.PLA2p -> PP2A + PLA2
    Reac143 = smodel.Reac('Reac143', vsys, lhs = [PP2A, PLA2star2], rhs = [PP2APLA2star2])
    Reac143.kcst = 1.4e6
    Reac144 = smodel.Reac('Reac144', vsys, lhs = [PP2APLA2star2], rhs = [PP2A, PLA2star2])
    Reac144.kcst = 1.5
    Reac145 = smodel.Reac('Reac145', vsys, lhs = [PP2APLA2star2], rhs = [PLA2, PP2A])
    Reac145.kcst = 2.5
    
    
    
    # PP2A + Ca1PLA2p <-> PP2A.Ca1PLA2p -> PP2A + Ca1PLA2
    Reac146 = smodel.Reac('Reac146', vsys, lhs = [PP2A, Ca1PLA2star2], rhs = [PP2ACa1PLA2star2])
    Reac146.kcst = 1.4e6
    Reac147 = smodel.Reac('Reac147', vsys, lhs = [PP2ACa1PLA2star2], rhs = [PP2A, Ca1PLA2star2])
    Reac147.kcst = 1.5 
    Reac148 = smodel.Reac('Reac148', vsys, lhs = [PP2ACa1PLA2star2], rhs = [PP2A, Ca1PLA2])
    Reac148.kcst = 2.5
    
    
    
    # PP2A + Ca2PLA2p <-> PP2A.Ca2PLA2p -> PP2A + Ca2PLA2
    Reac149 = smodel.Reac('Reac149', vsys, lhs = [PP2A, Ca2PLA2star2], rhs = [PP2ACa2PLA2star2])
    Reac149.kcst = 1.4e6
    Reac150 = smodel.Reac('Reac150', vsys, lhs = [PP2ACa2PLA2star2], rhs = [PP2A, Ca2PLA2star2])
    Reac150.kcst = 1.5 
    Reac151 = smodel.Reac('Reac151', vsys, lhs = [PP2ACa2PLA2star2], rhs = [PP2A, Ca2PLA2])
    Reac151.kcst = 2.5    
    
    
    
    # PP1 + PLA2p <-> PP1.PLA2p -> PP1 + PLA2
    Reac152 = smodel.Reac('Reac152', vsys, lhs = [PP1, PLA2star2], rhs = [PP1PLA2star2])
    Reac152.kcst = 1.4e6
    Reac153 = smodel.Reac('Reac153', vsys, lhs = [PP1PLA2star2], rhs = [PP1, PLA2star2])
    Reac153.kcst = 1.5
    Reac154 = smodel.Reac('Reac154', vsys, lhs = [PP1PLA2star2], rhs = [PLA2, PP1])
    Reac154.kcst = 2.5
    
    
    
    # PP1 + Ca1PLA2p <-> PP1.Ca1PLA2p -> PP1 + Ca1PLA2
    Reac155 = smodel.Reac('Reac155', vsys, lhs = [PP1, Ca1PLA2star2], rhs = [PP1Ca1PLA2star2])
    Reac155.kcst = 1.4e6
    Reac156 = smodel.Reac('Reac156', vsys, lhs = [PP1Ca1PLA2star2], rhs = [PP1, Ca1PLA2star2])
    Reac156.kcst = 1.5
    Reac157 = smodel.Reac('Reac157', vsys, lhs = [PP1Ca1PLA2star2], rhs = [PP1, Ca1PLA2])
    Reac157.kcst = 2.5
    
    
    
    # PP1 + Ca2PLA2p <-> PP1.Ca2PLA2p -> PP1 + Ca2PLA2
    Reac158 = smodel.Reac('Reac158', vsys, lhs = [PP1, Ca2PLA2star2], rhs = [PP1Ca2PLA2star2])
    Reac158.kcst = 1.4e6
    Reac159 = smodel.Reac('Reac159', vsys, lhs = [PP1Ca2PLA2star2], rhs = [PP1, Ca2PLA2star2])
    Reac159.kcst = 1.5
    Reac160 = smodel.Reac('Reac160', vsys, lhs = [PP1Ca2PLA2star2], rhs = [PP1, Ca2PLA2])
    Reac160.kcst = 2.5
    
    
    
    # PLA2p <-> PLA2** (PLA2star2memb) <-> PLA2**APC -> PLA2** + AA 
    Reac161 = smodel.SReac('Reac161', s, ilhs = [PLA2star2],  srhs = [PLA2star2memb])
    Reac161.kcst = 50
    Reac162 = smodel.SReac('Reac162', s, slhs = [PLA2star2memb], irhs = [PLA2star2])
    Reac162.kcst = 15
    Reac163 = smodel.SReac('Reac163', s, slhs = [PLA2star2memb],  srhs = [PLA2star2membAPC])
    Reac163.kcst = 43 
    Reac164 = smodel.SReac('Reac164', s, slhs = [PLA2star2membAPC],  srhs = [PLA2star2memb])
    Reac164.kcst = 600
    Reac165 = smodel.SReac('Reac165', s, slhs = [PLA2star2membAPC], srhs = [PLA2star2memb], irhs = [AA])
    Reac165.kcst = 3600 
    
    
    
    # Ca1PLA2p <-> Ca1PLA2** (Ca1PLA2star2memb) <-> Ca1PLA2**APC -> Ca1PLA2** + AA
    Reac166 = smodel.SReac('Reac166', s, ilhs = [Ca1PLA2star2], srhs = [Ca1PLA2star2memb])
    Reac166.kcst = 30
    Reac167 = smodel.SReac('Reac167', s, slhs = [Ca1PLA2star2memb], irhs = [Ca1PLA2star2])
    Reac167.kcst = 15
    Reac168 = smodel.SReac('Reac168', s, slhs = [Ca1PLA2star2memb],  srhs = [Ca1PLA2star2membAPC])
    Reac168.kcst = 43 
    Reac169 = smodel.SReac('Reac169', s, slhs = [Ca1PLA2star2membAPC],  srhs = [Ca1PLA2star2memb])
    Reac169.kcst = 600
    Reac170 = smodel.SReac('Reac170', s, slhs = [Ca1PLA2star2membAPC], srhs = [Ca1PLA2star2memb], irhs = [AA])
    Reac170.kcst = 3600
    
    
    
    # Ca2PLA2p <-> Ca2PLA2** (Ca2PLA2star2memb) <-> Ca2PLA2**APC -> Ca2PLA2** + AA
    Reac171 = smodel.SReac('Reac171', s, ilhs = [Ca2PLA2star2], srhs = [Ca2PLA2star2memb])
    Reac171.kcst = 300
    Reac172 = smodel.SReac('Reac172', s, slhs = [Ca2PLA2star2memb], irhs = [Ca2PLA2star2])
    Reac172.kcst = 15   
    Reac173 = smodel.SReac('Reac173', s, slhs = [Ca2PLA2star2memb],  srhs = [Ca2PLA2star2membAPC])
    Reac173.kcst = 43 
    Reac174 = smodel.SReac('Reac174', s, slhs = [Ca2PLA2star2membAPC],  srhs = [Ca2PLA2star2memb])
    Reac174.kcst = 600
    Reac175 = smodel.SReac('Reac175', s, slhs = [Ca2PLA2star2membAPC], srhs = [Ca2PLA2star2memb], irhs = [AA])
    Reac175.kcst = 3600
    
    
    
    # Ca + PLA2p <-> Ca1PLA2p
    Reac176 = smodel.Reac('Reac176', vsys, lhs = [PLA2star2, Ca], rhs = [Ca1PLA2star2])
    Reac176.kcst = 1.93e6 
    Reac177 = smodel.Reac('Reac177', vsys, lhs = [Ca1PLA2star2], rhs = [PLA2star2, Ca])
    Reac177.kcst = 108
    
    
    
    # Ca + PLA2** <-> Ca1PLA2**
    Reac178 = smodel.SReac('Reac178', s, slhs = [PLA2star2memb], ilhs = [Ca], srhs = [Ca1PLA2star2memb])
    Reac178.kcst = 1.93e6 
    Reac179 = smodel.SReac('Reac179', s, ilhs = [Ca1PLA2star2memb], srhs = [PLA2star2memb], irhs = [Ca])
    Reac179.kcst = 0.41
    
    
    
    # Ca + Ca1PLA2p <-> Ca2PLA2p
    Reac180 = smodel.Reac('Reac180', vsys, lhs = [Ca, Ca1PLA2star2], rhs = [Ca2PLA2star2])
    Reac180.kcst = 10.8e6 
    Reac181 = smodel.Reac('Reac181', vsys, lhs = [Ca2PLA2star2], rhs = [Ca1PLA2star2, Ca])
    Reac181.kcst = 108
    
    
    
    # Ca + Ca1PLA2** <-> Ca2PLA2**
    Reac182 = smodel.SReac('Reac182', s, slhs = [Ca1PLA2star2memb], ilhs = [Ca], srhs = [Ca2PLA2star2memb])
    Reac182.kcst = 10.8e6
    Reac183 = smodel.SReac('Reac183', s, ilhs = [Ca2PLA2star2memb], srhs = [Ca1PLA2star2memb], irhs = [Ca])
    Reac183.kcst = 2.5
    
    
    
    # AA -> 0 ,degradation
    Reac184 = smodel.SReac('Reac184', s, ilhs = [AA])
    Reac184.kcst = .4
    
    
    
    # Ca3PKC* + AMPAR <-> Ca3PKC*.AMPAR -> Ca3PKC* + AMPARp 
    Reac185 = smodel.SReac('Reac185', s, slhs = [PKCstar, AMPAR], srhs = [PKCstarAMPAR])
    # units for 2D reaction based on square meters in STEPS 1.2 and above
    # units for 2nd-order reaction are m^2/mol.s
    Reac185.kcst = 0.4e6*(1.02e-12/1.87e-19)
    Reac186 = smodel.SReac('Reac186', s, slhs = [PKCstarAMPAR], srhs = [PKCstar, AMPAR])
    Reac186.kcst = 0.8
    Reac187 = smodel.SReac('Reac187', s, slhs = [PKCstarAMPAR], srhs = [PKCstar, AMPAR_P])
    Reac187.kcst = 0.3
    
    
    
    # AAPKC* + AMPAR <-> AAPKC*.AMPAR -> AAPKC* + AMPARp
    Reac188 = smodel.SReac('Reac188', s, slhs = [PKCstar2, AMPAR], srhs = [PKCstar2AMPAR])
    # units for 2D reaction based on square meters in STEPS 1.2 and above
    # units for 2nd-order reaction are m^2/mol.s
    Reac188.kcst = 0.4e6*(1.02e-12/1.87e-19)
    Reac189 = smodel.SReac('Reac189', s, slhs = [PKCstar2AMPAR], srhs = [PKCstar2, AMPAR])
    Reac189.kcst = 0.8
    Reac190 = smodel.SReac('Reac190', s, slhs = [PKCstar2AMPAR], srhs = [PKCstar2, AMPAR_P])
    Reac190.kcst = 0.3
    
    
    
    # AACa1PKC* + AMPAR <-> AACa1PKC*.AMPAR -> AACa1PKC* + AMPARp
    Reac191 = smodel.SReac('Reac191', s, slhs = [PKCstar3, AMPAR], srhs = [PKCstar3AMPAR])
    # units for 2D reaction based on square meters in STEPS 1.2 and above
    # units for 2nd-order reaction are m^2/mol.s
    Reac191.kcst = 0.4e6*(1.02e-12/1.87e-19)   
    Reac192 = smodel.SReac('Reac192', s, slhs = [PKCstar3AMPAR], srhs = [PKCstar3, AMPAR])
    Reac192.kcst = 0.8
    Reac193 = smodel.SReac('Reac193', s, slhs = [PKCstar3AMPAR], srhs = [PKCstar3, AMPAR_P])
    Reac193.kcst =  0.3
    
    
    
    # AACa3PKC* + AMPAR <-> AACa3PKC*.AMPAR -> AACa3PKC* + AMPARp 
    Reac194 = smodel.SReac('Reac194', s, slhs = [PKCstar4, AMPAR], srhs = [PKCstar4AMPAR])
    # units for 2D reaction based on square meters in STEPS 1.2 and above
    # units for 2nd-order reaction are m^2/mol.s
    Reac194.kcst = 0.4e6*(1.02e-12/1.87e-19)  
    Reac195 = smodel.SReac('Reac195', s, slhs = [PKCstar4AMPAR], srhs = [PKCstar4, AMPAR])
    Reac195.kcst = 0.8 
    Reac196 = smodel.SReac('Reac196', s, slhs = [PKCstar4AMPAR], srhs = [PKCstar4, AMPAR_P])
    Reac196.kcst = 0.3
    
    
    
    # PP2A + AMPARp <-> PP2A.AMPARp -> PP2A + AMPAR
    Reac197 = smodel.SReac('Reac197', s, ilhs = [PP2A], slhs = [AMPAR_P], srhs = [PP2AAMPAR_P])
    Reac197.kcst = 0.6e6
    Reac198 = smodel.SReac('Reac198', s, slhs = [PP2AAMPAR_P], irhs = [PP2A], srhs = [AMPAR_P])
    Reac198.kcst = 0.17
    Reac199 = smodel.SReac('Reac199', s, slhs = [PP2AAMPAR_P], srhs = [AMPAR], irhs = [PP2A])
    Reac199.kcst = 0.25 
    
    
    
    # Ca3PKC* + GRIP.AMPA <-> Ca3PKC*.GRIP.AMPA -> Ca3PKC* + GRIP.AMPAp 
    Reac200 = smodel.SReac('Reac200', s, slhs = [PKCstar, GRIPAMPAR], srhs = [PKCstarGRIPAMPAR])
    # units for 2D reaction based on square meters in STEPS 1.2 and above
    # units for 2nd-order reaction are m^2/mol.s
    Reac200.kcst = 0.4e6*(1.02e-12/1.87e-19)
    Reac201 = smodel.SReac('Reac201', s, slhs = [PKCstarGRIPAMPAR], srhs = [PKCstar, GRIPAMPAR])
    Reac201.kcst = 0.8
    Reac202 = smodel.SReac('Reac202', s, slhs = [PKCstarGRIPAMPAR], srhs = [PKCstar, GRIPAMPAR_P])
    Reac202.kcst = 0.3
    
    
    
    # AAPKC* + GRIP.AMPA <-> AAPKC*.GRIP.AMPA -> AAPKC* + GRIP.AMPAp
    Reac203 = smodel.SReac('Reac203', s, slhs = [PKCstar2, GRIPAMPAR], srhs = [PKCstar2GRIPAMPAR])
    # units for 2D reaction based on square meters in STEPS 1.2 and above
    # units for 2nd-order reaction are m^2/mol.s
    Reac203.kcst = 0.4e6*(1.02e-12/1.87e-19)
    Reac204 = smodel.SReac('Reac204', s, slhs = [PKCstar2GRIPAMPAR], srhs = [PKCstar2, GRIPAMPAR])
    Reac204.kcst = 0.8
    Reac205 = smodel.SReac('Reac205', s, slhs = [PKCstar2GRIPAMPAR], srhs = [PKCstar2, GRIPAMPAR_P])
    Reac205.kcst = 0.3
    
    
    
    # AACa1PKC* + GRIP.AMPA <-> AACa1PKC*.GRIP.AMPA -> AACa1PKC* + GRIP.AMPAp
    Reac206 = smodel.SReac('Reac206', s, slhs = [PKCstar3, GRIPAMPAR], srhs = [PKCstar3GRIPAMPAR])
    # units for 2D reaction based on square meters in STEPS 1.2 and above
    # units for 2nd-order reaction are m^2/mol.s
    Reac206.kcst = 0.4e6*(1.02e-12/1.87e-19)   
    Reac207 = smodel.SReac('Reac207', s, slhs = [PKCstar3GRIPAMPAR], srhs = [PKCstar3, GRIPAMPAR])
    Reac207.kcst = 0.8
    Reac208 = smodel.SReac('Reac208', s, slhs = [PKCstar3GRIPAMPAR], srhs = [PKCstar3, GRIPAMPAR_P])
    Reac208.kcst =  0.3
    
    
    
    # AACa3PKC* + GRIP.AMPA <-> AACa3PKC*.GRIP.AMPA -> AACa3PKC* + GRIP.AMPAp 
    Reac209 = smodel.SReac('Reac209', s, slhs = [PKCstar4, GRIPAMPAR], srhs = [PKCstar4GRIPAMPAR])
    # units for 2D reaction based on square meters in STEPS 1.2 and above
    # units for 2nd-order reaction are m^2/mol.s
    Reac209.kcst = 0.4e6*(1.02e-12/1.87e-19)  
    Reac210 = smodel.SReac('Reac210', s, slhs = [PKCstar4GRIPAMPAR], srhs = [PKCstar4, GRIPAMPAR])
    Reac210.kcst = 0.8 
    Reac211 = smodel.SReac('Reac211', s, slhs = [PKCstar4GRIPAMPAR], srhs = [PKCstar4, GRIPAMPAR_P])
    Reac211.kcst = 0.3
    
    
    
    # PP2A + GRIP.AMPAp <-> PP2A.GRIP.AMPAp -> PP2A + GRIP.AMPA
    Reac212 = smodel.SReac('Reac212', s, ilhs = [PP2A], slhs = [GRIPAMPAR_P], srhs = [PP2AGRIPAMPAR_P])
    Reac212.kcst = 0.6e6
    Reac213 = smodel.SReac('Reac213', s, slhs = [PP2AGRIPAMPAR_P], irhs = [PP2A], srhs = [GRIPAMPAR_P])
    Reac213.kcst = 0.17
    Reac214 = smodel.SReac('Reac214', s, slhs = [PP2AGRIPAMPAR_P], srhs = [GRIPAMPAR], irhs = [PP2A])
    Reac214.kcst = 0.25 
    
    
    
    # AMPAR + GRIP <-> GRIP.AMPA
    Reac215 = smodel.SReac('Reac215', s, slhs = [AMPAR, GRIP], srhs = [GRIPAMPAR])
    # units for 2D reaction based on square meters in STEPS 1.2 and above
    # units for 2nd-order reaction are m^2/mol.s
    Reac215.kcst = 1e6*(1.02e-12/1.87e-19)
    Reac216 = smodel.SReac('Reac216', s, slhs = [GRIPAMPAR], srhs = [GRIP, AMPAR])
    Reac216.kcst = 7
    
    
    
    # AMPARp + GRIP <-> GRIP.AMPAp (AMPAR = synaptic AMPAR)
    Reac217 = smodel.SReac('Reac217', s, slhs = [AMPAR_P, GRIP], srhs = [GRIPAMPAR_P])
    # units for 2D reaction based on square meters in STEPS 1.2 and above
    # units for 2nd-order reaction are m^2/mol.s
    Reac217.kcst = 1e6*(1.02e-12/1.87e-19)
    Reac218 = smodel.SReac('Reac218', s, slhs = [GRIPAMPAR_P], srhs = [GRIP, AMPAR_P])
    Reac218.kcst = 70
    
    
    
    # AMPAR <-> AMPARextra (AMPARextra = extra-synaptic AMPAR)  
    Reac219 = smodel.SReac('Reac219', s, slhs = [AMPAR], srhs = [AMPARextra])
    Reac219.kcst = 0.1
    Reac220 = smodel.SReac('Reac220', s, slhs = [AMPARextra], srhs = [AMPAR])
    Reac220.kcst = 0.02
    
    
    
    # AMPARp <-> AMPARextra_P
    Reac221 = smodel.SReac('Reac221', s, slhs = [AMPAR_P], srhs = [AMPARextra_P])
    Reac221.kcst =  0.1
    Reac222 = smodel.SReac('Reac222', s, slhs = [AMPARextra_P], srhs = [AMPAR_P])
    Reac222.kcst = 0.02
    
    
    
    # PP2A + AMPARextra_P  <-> PP2A.AMPARextra_P  -> PP2A + AMPARextra
    Reac223 = smodel.SReac('Reac223', s, ilhs = [PP2A], slhs = [AMPARextra_P], srhs = [PP2AAMPARextra_P])
    Reac223.kcst = 0.6e6
    Reac224 = smodel.SReac('Reac224', s, slhs = [PP2AAMPARextra_P], irhs = [PP2A], srhs = [AMPARextra_P])
    Reac224.kcst = 0.17
    Reac225 = smodel.SReac('Reac225', s, slhs = [PP2AAMPARextra_P], srhs = [AMPARextra], irhs = [PP2A])
    Reac225.kcst = 0.25
    
    
    
    # AMPARextra <-> AMPARdend (AMPARdend = dendritic AMPAR)  
    Reac226 = smodel.SReac('Reac226', s, slhs = [AMPARextra], srhs = [AMPARdend])
    Reac226.kcst = 0.02
    Reac227 = smodel.SReac('Reac227', s, slhs = [AMPARdend], srhs = [AMPARextra])
    Reac227.kcst = 0.00025
    
    
    
    # AMPARextra_P <-> AMPARdend_P 
    Reac228 = smodel.SReac('Reac228', s, slhs = [AMPARextra_P], srhs = [AMPARdend_P])
    Reac228.kcst =  0.02
    Reac229 = smodel.SReac('Reac229', s, slhs = [AMPARdend_P], srhs = [AMPARextra_P])
    Reac229.kcst = 0.00025
    
    
    
    # PP2A + AMPARp_dend  <-> PP2A.AMPARp_dend  -> PP2A + AMPAR_dend
    Reac230 = smodel.SReac('Reac230', s, ilhs = [PP2A], slhs = [AMPARdend_P], srhs = [PP2AAMPARdend_P])
    Reac230.kcst = 0.6e6
    Reac231 = smodel.SReac('Reac231', s, slhs = [PP2AAMPARdend_P], irhs = [PP2A], srhs = [AMPARdend_P])
    Reac231.kcst = 0.17
    Reac232 = smodel.SReac('Reac232', s, slhs = [PP2AAMPARdend_P], srhs = [AMPARdend], irhs = [PP2A])
    Reac232.kcst = 0.25
    
    
    
    # AMPARdend_P <-> AMPARcyt_P (AMPARcyt = cytosolic AMPAR) 
    Reac233 = smodel.SReac('Reac233', s, slhs = [AMPARdend_P], irhs = [AMPARcyt_P])
    Reac233.kcst =  0.003
    Reac234 = smodel.SReac('Reac234', s, ilhs = [AMPARcyt_P], srhs = [AMPARdend_P])
    Reac234.kcst =  0.002
    
    
    
    # PP2A + AMPARcyt_P  <-> PP2A.AMPARcyt_P  -> PP2A + AMPARcyt 
    Reac235 = smodel.Reac('Reac235', vsys, lhs = [PP2A, AMPARcyt_P], rhs = [PP2AAMPARcyt_P])
    Reac235.kcst = 0.6e6
    Reac236 = smodel.Reac('Reac236', vsys, lhs = [PP2AAMPARcyt_P], rhs = [PP2A, AMPARcyt_P])
    Reac236.kcst = 0.17
    Reac237 = smodel.Reac('Reac237', vsys, lhs = [PP2AAMPARcyt_P], rhs = [AMPARcyt, PP2A])
    Reac237.kcst = 0.25
    
    return mdl

########################################################################

# Geometric properties of the model: To change the size of the model, 
# all volumetric compartments must be altered by the same ratio
# and all the areas must be scaled considering a spherical shape.
# Alterations in the size of the compartments will change automatically the population
# of species given in concentration (initial condition), but the species 
# given in number of molecules must be altered manually by the same ratio
# to keep the balance among all the components of the model. 

def gen_geom():    
    g = swm.Geom()
    c = swm.Comp('vsys', g)
    c.addVolsys('vsys')
    c.vol = 0.08e-18
    
    extra = swm.Comp('extra', g)
    extra.addVolsys('extra')
    extra.vol = 1.87e-22
    
    s = swm.Patch('memb', g, icomp = c, ocomp = extra)
    s.addSurfsys('memb')
    s.area = 1.02e-12
    
    cytER = swm.Comp('cytER', g)
    cytER.addVolsys('cytER')
    cytER.vol = 0.017e-18
    
    ERs = swm.Patch('ERmemb', g, icomp = cytER, ocomp = c)
    ERs.addSurfsys('ERmemb')
    ERs.area = 0.32e-12
    
    return g

######################################################################

def run_simulation():
    rng=srng.create('mt19937',512)
    m=gen_model() 
    g=gen_geom() 
    tpnts=numpy.arange(0.0,INT,DT)
    ntpnts=tpnts.shape[0]
    
    # initialize arrays for storage at each time point and at each run.
    Ca = numpy.zeros((NITER, ntpnts))
    
    PMCA = numpy.zeros((NITER, ntpnts))
    Ca1PMCA = numpy.zeros((NITER, ntpnts))
    NCX = numpy.zeros((NITER, ntpnts))
    Ca1NCX = numpy.zeros((NITER, ntpnts))
    Ca2NCX = numpy.zeros((NITER, ntpnts))
    SERCA = numpy.zeros((NITER, ntpnts))
    Ca1SERCA = numpy.zeros((NITER, ntpnts))
    Ca2SERCA = numpy.zeros((NITER, ntpnts))
    CBf = numpy.zeros((NITER, ntpnts))
    CaCBf = numpy.zeros((NITER, ntpnts))
    Ca2CBf = numpy.zeros((NITER, ntpnts))
    CBs = numpy.zeros((NITER, ntpnts))
    CaCBs = numpy.zeros((NITER, ntpnts))
    Ca2CBs = numpy.zeros((NITER, ntpnts))
    PV = numpy.zeros((NITER, ntpnts))
    CaPV = numpy.zeros((NITER, ntpnts))
    MgPV = numpy.zeros((NITER, ntpnts))
    Ca2PV = numpy.zeros((NITER, ntpnts))
    Mg2PV = numpy.zeros((NITER, ntpnts))
    PP2A = numpy.zeros((NITER, ntpnts))
    PKC = numpy.zeros((NITER, ntpnts))
    Ca1PKC = numpy.zeros((NITER, ntpnts))
    Ca3PKC = numpy.zeros((NITER, ntpnts))
    AAPKC = numpy.zeros((NITER, ntpnts))
    AACa1PKC = numpy.zeros((NITER, ntpnts))
    AACa3PKC = numpy.zeros((NITER, ntpnts))
    PKCstar = numpy.zeros((NITER, ntpnts))
    PKCstar2 = numpy.zeros((NITER, ntpnts))
    PKCstar4 = numpy.zeros((NITER, ntpnts))
    PKCstar3 = numpy.zeros((NITER, ntpnts))
    RKIPstarRP = numpy.zeros((NITER, ntpnts))
    RafRKIPPKCstar = numpy.zeros((NITER, ntpnts))
    RafRKIPPKCstar2 = numpy.zeros((NITER, ntpnts))
    RafRKIPPKCstar3 = numpy.zeros((NITER, ntpnts))
    RafRKIPPKCstar4 = numpy.zeros((NITER, ntpnts))
    RafRKIP = numpy.zeros((NITER, ntpnts))
    MEKRKIPPKCstar = numpy.zeros((NITER, ntpnts))
    MEKRKIPPKCstar2 = numpy.zeros((NITER, ntpnts))
    MEKRKIPPKCstar3 = numpy.zeros((NITER, ntpnts))
    MEKRKIPPKCstar4 = numpy.zeros((NITER, ntpnts))
    MEKRKIP = numpy.zeros((NITER, ntpnts))
    Raf = numpy.zeros((NITER, ntpnts))
    RKIPstar = numpy.zeros((NITER, ntpnts))
    Rafstar = numpy.zeros((NITER, ntpnts))
    PP5Rafstar = numpy.zeros((NITER, ntpnts))
    PP5 = numpy.zeros((NITER, ntpnts))
    MEK = numpy.zeros((NITER, ntpnts))
    RafstarMEK = numpy.zeros((NITER, ntpnts))
    MEKp = numpy.zeros((NITER, ntpnts))
    PP2AMEKp = numpy.zeros((NITER, ntpnts))
    RafstarMEKp = numpy.zeros((NITER, ntpnts))
    MEKstar = numpy.zeros((NITER, ntpnts))
    PP2AMEKstar = numpy.zeros((NITER, ntpnts))
    ERK = numpy.zeros((NITER, ntpnts))
    MEKstarERK = numpy.zeros((NITER, ntpnts))
    ERKp = numpy.zeros((NITER, ntpnts))
    MKP = numpy.zeros((NITER, ntpnts))
    MKPERKp = numpy.zeros((NITER, ntpnts))
    MEKstarERKp = numpy.zeros((NITER, ntpnts))
    ERKstar = numpy.zeros((NITER, ntpnts))
    MKPERKstar = numpy.zeros((NITER, ntpnts))
    Ca1PLA2 = numpy.zeros((NITER, ntpnts))
    Ca2PLA2 = numpy.zeros((NITER, ntpnts))
    PLA2 = numpy.zeros((NITER, ntpnts))
    PLA2memb = numpy.zeros((NITER, ntpnts))
    Ca1PLA2memb = numpy.zeros((NITER, ntpnts))
    PLA2star1 = numpy.zeros((NITER, ntpnts))
    ERKstarPLA2 = numpy.zeros((NITER, ntpnts))
    PLA2star2 = numpy.zeros((NITER, ntpnts))
    AA = numpy.zeros((NITER, ntpnts))
    Ca1PLA2star2 = numpy.zeros((NITER, ntpnts))
    Ca2PLA2star2 = numpy.zeros((NITER, ntpnts))
    PLA2star2memb = numpy.zeros((NITER, ntpnts))
    Ca1PLA2star2memb = numpy.zeros((NITER, ntpnts))
    Ca2PLA2star2memb = numpy.zeros((NITER, ntpnts))
    PLA2star1APC = numpy.zeros((NITER, ntpnts))
    ERKstarCa1PLA2 = numpy.zeros((NITER, ntpnts))
    ERKstarCa2PLA2 = numpy.zeros((NITER, ntpnts))
    PP2ACa1PLA2star2 = numpy.zeros((NITER, ntpnts))
    PP2ACa2PLA2star2 = numpy.zeros((NITER, ntpnts))
    PP1Ca1PLA2star2 = numpy.zeros((NITER, ntpnts))
    PP1Ca2PLA2star2 = numpy.zeros((NITER, ntpnts))
    PLA2star2membAPC = numpy.zeros((NITER, ntpnts))
    Ca1PLA2star2membAPC = numpy.zeros((NITER, ntpnts))
    Ca2PLA2star2membAPC = numpy.zeros((NITER, ntpnts))
    PP1PLA2star2 = numpy.zeros((NITER, ntpnts))
    PP1 = numpy.zeros((NITER, ntpnts)) 
    PP2APLA2star2 = numpy.zeros((NITER, ntpnts))
    
    PKCstarAMPAR = numpy.zeros((NITER, ntpnts))
    PKCstar2AMPAR = numpy.zeros((NITER, ntpnts))
    PKCstar4AMPAR = numpy.zeros((NITER, ntpnts))
    PKCstar3AMPAR = numpy.zeros((NITER, ntpnts))

    PKCstarGRIPAMPAR = numpy.zeros((NITER, ntpnts))
    PKCstar2GRIPAMPAR = numpy.zeros((NITER, ntpnts))
    PKCstar4GRIPAMPAR = numpy.zeros((NITER, ntpnts))
    PKCstar3GRIPAMPAR = numpy.zeros((NITER, ntpnts))
    PP2AGRIPAMPAR_P = numpy.zeros((NITER, ntpnts))
    
    AMPAR = numpy.zeros((NITER, ntpnts))
    AMPAR_P = numpy.zeros((NITER, ntpnts))
    PP2AAMPAR_P = numpy.zeros((NITER, ntpnts))
    AMPARextra = numpy.zeros((NITER, ntpnts))
    AMPARextra_P = numpy.zeros((NITER, ntpnts))
    PP2AAMPARextra_P = numpy.zeros((NITER, ntpnts))
    GRIP = numpy.zeros((NITER, ntpnts))
    GRIPAMPAR_P = numpy.zeros((NITER, ntpnts))
    GRIPAMPAR = numpy.zeros((NITER, ntpnts))
    AMPARdend = numpy.zeros((NITER, ntpnts))
    AMPARdend_P = numpy.zeros((NITER, ntpnts))
    PP2AAMPARdend_P = numpy.zeros((NITER, ntpnts))
    AMPARcyt = numpy.zeros((NITER, ntpnts))
    AMPARcyt_P = numpy.zeros((NITER, ntpnts))
    PP2AAMPARcyt_P = numpy.zeros((NITER, ntpnts))
    
    sim = swmdirect.Wmdirect(m, g, rng)
    
    rng.initialize(int(idx)*100)
    
    for j in arange(NITER): 
        print "Run number {0}".format(j)
        sim.reset()
        
        # initial conditions
        sim.setCompConc('vsys', 'Ca', 0.045e-6)
        sim.setCompConc('vsys', 'CBf', 37.775e-6)
        sim.setCompConc('vsys', 'CaCBf', 2.1e-6)
        sim.setCompConc('vsys', 'Ca2CBf', 0.125e-6)
        sim.setCompConc('vsys', 'CBs', 36.25e-6)
        sim.setCompConc('vsys', 'CaCBs', 3.4e-6)
        sim.setCompConc('vsys', 'Ca2CBs', 0.125e-6)
        sim.setCompConc('vsys', 'PV', 1.15e-6)
        sim.setCompConc('vsys', 'CaPV', 8.4e-6)
        sim.setCompConc('vsys', 'MgPV', 30.45e-6)

        sim.setCompCount('vsys', 'PKC', int(pkc_number))
        
        sim.setCompCount('vsys', 'Raf',int(raf_number))
        
        sim.setCompConc('vsys', 'RKIP', 1.0e-6)
        sim.setCompConc('vsys', 'MEK',1.5e-6) 
        sim.setCompConc('vsys', 'RP' ,3e-6)
        sim.setCompConc('vsys', 'ERK',1.0e-6)
        sim.setCompConc('vsys', 'MKP',0.26e-6) 
        sim.setCompConc('vsys', 'PP5', 1.0e-6) 
        
        sim.setCompCount('vsys', 'PP2A', 35)

        sim.setCompConc('vsys', 'PLA2',0.4e-6)
        
        sim.setCompCount('vsys', 'PP1', 30)

        sim.setPatchCount('memb', 'PMCA', 10)
        sim.setPatchCount('memb', 'NCX', 3)
        sim.setPatchCount('ERmemb', 'SERCA', 5)

        sim.setCompConc('cytER', 'Ca', 150e-6)
        sim.setCompClamped('cytER', 'Ca', True) # clamped means the conc won't change as simulation runs.
        sim.setPatchCount('memb', 'AMPAR', 3) 
        sim.setPatchCount('memb', 'AMPARextra', 16) 
        sim.setPatchCount('memb', 'GRIP', 22) 
        sim.setPatchCount('memb', 'GRIPAMPAR', 72) 
        sim.setPatchCount('memb', 'AMPARdend', 1600)
        
        for i in range(ntpnts):
            
            sim.run(tpnts[i]) 
            if not i%1000: print "idx:", idx, " percent done", (tpnts[i]/INT*100)
            
            Ca[j,i] = sim.getCompConc('vsys', 'Ca')
            
            CBf[j,i] = sim.getCompConc('vsys', 'CBf')
            CaCBf[j,i] = sim.getCompConc('vsys', 'CaCBf')
            Ca2CBf[j,i] = sim.getCompConc('vsys', 'Ca2CBf')
            CaCBs[j,i] = sim.getCompConc('vsys', 'CaCBs')
            Ca2CBs[j,i] = sim.getCompConc('vsys', 'Ca2CBs')
            PV[j,i] = sim.getCompConc('vsys', 'PV')
            CaPV[j,i] = sim.getCompConc('vsys', 'CaPV')
            MgPV[j,i] = sim.getCompConc('vsys', 'MgPV')
            Ca2PV[j,i] = sim.getCompConc('vsys', 'Ca2PV')
            Mg2PV[j,i] = sim.getCompConc('vsys', 'Mg2PV')
            PP2A[j,i] = sim.getCompCount('vsys', 'PP2A')
            PP5[j,i] = sim.getCompCount('vsys', 'PP5')
            PKC[j,i] = sim.getCompCount('vsys', 'PKC')
            Ca3PKC[j,i] = sim.getCompCount('vsys', 'Ca3PKC')
            Ca1PKC[j,i] = sim.getCompCount('vsys', 'Ca1PKC')
            AAPKC[j,i] = sim.getCompCount('vsys', 'AAPKC')
            AACa1PKC[j,i] = sim.getCompCount('vsys', 'AACa1PKC')
            AACa3PKC[j,i] = sim.getCompCount('vsys', 'AACa3PKC')
            PKCstar[j,i] = sim.getPatchCount('memb', 'PKCstar')
            PKCstar2[j,i] = sim.getPatchCount('memb', 'PKCstar2')
            PKCstar3[j,i] = sim.getPatchCount('memb', 'PKCstar3')
            PKCstar4[j,i] = sim.getPatchCount('memb', 'PKCstar4')
            RafRKIP[j,i] = sim.getCompConc('vsys', 'RafRKIP')
            RafRKIPPKCstar[j,i] = sim.getPatchCount('memb', 'RafRKIPPKCstar')
            RafRKIPPKCstar2[j,i] = sim.getPatchCount('memb', 'RafRKIPPKCstar2')
            RafRKIPPKCstar3[j,i] = sim.getPatchCount('memb', 'RafRKIPPKCstar3')
            RafRKIPPKCstar4[j,i] = sim.getPatchCount('memb', 'RafRKIPPKCstar4')
            
            MEKRKIP[j,i] = sim.getCompConc('vsys', 'MEKRKIP')
            MEKRKIPPKCstar[j,i] = sim.getPatchCount('memb', 'MEKRKIPPKCstar')
            MEKRKIPPKCstar2[j,i] = sim.getPatchCount('memb', 'MEKRKIPPKCstar2')
            MEKRKIPPKCstar3[j,i] = sim.getPatchCount('memb', 'MEKRKIPPKCstar3')
            MEKRKIPPKCstar4[j,i] = sim.getPatchCount('memb', 'MEKRKIPPKCstar4')
            
            RKIPstar[j,i] = sim.getCompCount('vsys', 'RKIPstar')
            RKIPstarRP[j,i] = sim.getCompCount('vsys', 'RKIPstarRP')
            Raf[j,i] = sim.getCompCount('vsys', 'Raf')
            Rafstar[j,i] = sim.getCompCount('vsys', 'Rafstar')
            PP5Rafstar[j,i] = sim.getCompCount('vsys', 'PP5Rafstar')
            MEK[j,i] = sim.getCompCount('vsys', 'MEK')
            RafstarMEK[j,i] = sim.getCompCount('vsys', 'RafstarMEK')
            MEKp[j,i] = sim.getCompCount('vsys', 'MEKp')
            PP2AMEKp[j,i] = sim.getCompCount('vsys', 'PP2AMEKp')
            RafstarMEKp[j,i] = sim.getCompCount('vsys', 'RafstarMEKp')
            MEKstar[j,i] = sim.getCompCount('vsys', 'MEKstar')
            ERK[j,i] = sim.getCompCount('vsys', 'ERK')
            MEKstarERK[j,i] = sim.getCompCount('vsys', 'MEKstarERK')
            ERKp[j,i] = sim.getCompCount('vsys', 'ERKp')
            MKP[j,i] = sim.getCompCount('vsys', 'MKP')
            MKPERKp[j,i] = sim.getCompCount('vsys', 'MKPERKp')
            MEKstarERKp[j,i] = sim.getCompCount('vsys', 'MEKstarERKp')
            ERKstar[j,i] = sim.getCompCount('vsys', 'ERKstar')
            MKPERKstar[j,i] = sim.getCompCount('vsys', 'MKPERKstar')
            PP2A[j,i] = sim.getCompCount('vsys', 'PP2A')
            PLA2[j,i] = sim.getCompCount('vsys', 'PLA2')
            Ca1PLA2[j,i] = sim.getCompCount('vsys', 'Ca1PLA2')
            Ca2PLA2[j,i] = sim.getCompCount('vsys', 'Ca2PLA2')
            PLA2memb[j,i] = sim.getPatchCount('memb', 'PLA2memb')
            Ca1PLA2memb[j,i] = sim.getPatchCount('memb', 'Ca1PLA2memb')
            PLA2star1[j,i] = sim.getPatchCount('memb', 'PLA2star1')
            ERKstarPLA2[j,i] = sim.getCompCount('vsys', 'ERKstarPLA2')
            PLA2star2[j,i] = sim.getCompCount('vsys', 'PLA2star2')
            Ca1PLA2star2[j,i] = sim.getCompCount('vsys', 'Ca1PLA2star2')
            Ca2PLA2star2[j,i] = sim.getCompCount('vsys', 'Ca2PLA2star2')
            PLA2star2memb[j,i] = sim.getPatchCount('memb', 'PLA2star2memb')
            Ca1PLA2star2memb[j,i] = sim.getPatchCount('memb', 'Ca1PLA2star2memb')
            Ca2PLA2star2memb[j,i] = sim.getPatchCount('memb', 'Ca2PLA2star2memb')  
            AA[j,i] = sim.getCompCount('vsys', 'AA')
            PP2APLA2star2[j,i] = sim.getCompCount('vsys', 'PP2APLA2star2')
            PP1[j,i] = sim.getCompCount('vsys','PP1')
            PP1PLA2star2[j,i] = sim.getCompCount('vsys', 'PP1PLA2star2')
            PLA2star1APC[j,i] = sim.getPatchCount('memb', 'PLA2star1APC')
            ERKstarCa1PLA2[j,i] = sim.getCompConc('vsys', 'ERKstarCa1PLA2')
            ERKstarCa2PLA2[j,i] = sim.getCompConc('vsys', 'ERKstarCa2PLA2')             
            PP2ACa1PLA2star2[j,i] = sim.getCompConc('vsys', 'PP2ACa1PLA2star2')
            PP2ACa2PLA2star2[j,i] = sim.getCompConc('vsys', 'PP2ACa2PLA2star2')
            PP1Ca1PLA2star2[j,i] = sim.getCompConc('vsys', 'PP1Ca1PLA2star2')
            PP1Ca2PLA2star2[j,i] = sim.getCompConc('vsys', 'PP1Ca2PLA2star2')
            PLA2star2membAPC[j,i] = sim.getPatchCount('memb', 'PLA2star2membAPC')
            Ca1PLA2star2membAPC[j,i] = sim.getPatchCount('memb', 'Ca1PLA2star2membAPC')
            Ca2PLA2star2membAPC[j,i] = sim.getPatchCount('memb', 'Ca2PLA2star2membAPC')
            Ca1PMCA[j,i] = sim.getPatchCount('memb', 'Ca1PMCA')
            PMCA[j,i] = sim.getPatchCount('memb', 'PMCA')
            NCX[j,i] = sim.getPatchCount('memb', 'NCX')
            Ca2NCX[j,i] = sim.getPatchCount('memb', 'Ca2NCX')
            SERCA[j,i] = sim.getPatchCount('ERmemb', 'SERCA')
            Ca2SERCA[j,i] = sim.getPatchCount('ERmemb', 'Ca2SERCA')
            PP1PLA2star2[j,i] = sim.getCompConc('vsys', 'PP1PLA2star2')
            
            AMPAR[j,i] = sim.getPatchCount('memb', 'AMPAR')  
            AMPAR_P[j,i] = sim.getPatchCount('memb', 'AMPAR_P') 
            PP2AAMPAR_P[j,i] = sim.getPatchCount('memb', 'PP2AAMPAR_P') 
            AMPARextra[j,i] = sim.getPatchCount('memb', 'AMPARextra')
            AMPARextra_P[j,i] = sim.getPatchCount('memb', 'AMPARextra_P')
            PP2AAMPARextra_P[j,i] = sim.getPatchCount('memb', 'PP2AAMPARextra_P')
            GRIP[j,i] = sim.getPatchCount('memb', 'GRIP')
            GRIPAMPAR[j,i] = sim.getPatchCount('memb', 'GRIPAMPAR')
            GRIPAMPAR_P[j,i] = sim.getPatchCount('memb', 'GRIPAMPAR_P')
            PKCstarAMPAR[j,i] = sim.getPatchCount('memb', 'PKCstarAMPAR')  
            PKCstar2AMPAR[j,i] = sim.getPatchCount('memb', 'PKCstar2AMPAR') 
            PKCstar3AMPAR[j,i] = sim.getPatchCount('memb', 'PKCstar3AMPAR') 
            PKCstar4AMPAR[j,i] = sim.getPatchCount('memb', 'PKCstar4AMPAR')   
            AMPARdend[j,i] = sim.getPatchCount('memb', 'AMPARdend')
            AMPARdend_P[j,i] = sim.getPatchCount('memb', 'AMPARdend_P')
            PP2AAMPARdend_P[j,i] = sim.getPatchCount('memb', 'PP2AAMPARdend_P')
            AMPARcyt[j,i] = sim.getCompCount('vsys','AMPARcyt')
            AMPARcyt_P[j,i] = sim.getCompCount('vsys', 'AMPARcyt_P')
            PP2AAMPARcyt_P[j,i] = sim.getCompCount('vsys', 'PP2AAMPARcyt_P')   
            
            PKCstarGRIPAMPAR[j,i] =  sim.getPatchCount('memb', 'PKCstarGRIPAMPAR')
            PKCstar2GRIPAMPAR[j,i] =  sim.getPatchCount('memb', 'PKCstar2GRIPAMPAR')
            PKCstar4GRIPAMPAR[j,i]  = sim.getPatchCount('memb', 'PKCstar4GRIPAMPAR')
            PKCstar3GRIPAMPAR[j,i] = sim.getPatchCount('memb', 'PKCstar3GRIPAMPAR')
            PP2AGRIPAMPAR_P[j,i] = sim.getPatchCount('memb', 'PP2AGRIPAMPAR_P')
            
            sim.setCompReacK('vsys', 'Cainflux', square_pulse_ica(tpnts[i])) # changes reaction constant to be ...
    
#############################################################################
    
    if 'AMPAR' in data_record:
    
        AMPAR_mean = mean(AMPAR, axis=0)
        AMPAR_P_mean = mean(AMPAR_P, axis=0)
        PP2AAMPAR_P_mean = mean(PP2AAMPAR_P, axis=0)
        GRIPAMPAR_mean = mean(GRIPAMPAR, axis=0)
        GRIPAMPAR_P_mean = mean(GRIPAMPAR_P, axis=0)
        PKCstarAMPAR_mean = mean(PKCstarAMPAR, axis=0)
        PKCstar2AMPAR_mean = mean(PKCstar2AMPAR, axis=0)
        PKCstar3AMPAR_mean = mean(PKCstar3AMPAR, axis=0)
        PKCstar4AMPAR_mean = mean(PKCstar4AMPAR, axis=0)

        PKCstarGRIPAMPAR_mean = mean(PKCstarGRIPAMPAR, axis=0)
        PKCstar2GRIPAMPAR_mean = mean(PKCstar2GRIPAMPAR, axis=0)
        PKCstar3GRIPAMPAR_mean = mean(PKCstar3GRIPAMPAR, axis=0)
        PKCstar4GRIPAMPAR_mean = mean(PKCstar4GRIPAMPAR, axis=0)
        PP2AGRIPAMPAR_P_mean = mean(PP2AGRIPAMPAR_P, axis=0)

        AMPAR_syn_mean = AMPAR_mean + AMPAR_P_mean + PP2AAMPAR_P_mean + \
                         GRIPAMPAR_mean + GRIPAMPAR_P_mean + PKCstarAMPAR_mean + \
                         PKCstar2AMPAR_mean + PKCstar3AMPAR_mean + PKCstar4AMPAR_mean +\
                         PKCstarGRIPAMPAR_mean + PKCstar2GRIPAMPAR_mean + PKCstar3GRIPAMPAR_mean +\
                         PKCstar4GRIPAMPAR_mean + PP2AGRIPAMPAR_P_mean;
                         
        ofile_ampar=open(datapath +'/'+ str(idx) +'/AMPAR_data.txt','w')
        pickle.dump(AMPAR_syn_mean, ofile_ampar)
        
        ofile_ampar.close()     
        
        if plotfigs:
            plot(linspace(0,INT,INT/DT), AMPAR_syn_mean, 'g')
            ylim(ymin=0.0)
            xlabel('Time (s)')
            ylabel('Synaptic AMPAR (#)')
            savefig(datapath +'/'+ str(idx) +'/Syn_AMPAR.png')
            close()
        
    if 'Ca' in data_record:
    
        Ca_mean = mean(Ca, axis=0)
                         
        ofile_ca=open(datapath +'/'+ str(idx) +'/Ca_data.txt','w')
        pickle.dump(Ca_mean, ofile_ca)
        ofile_ca.close()
        
        if plotfigs:                        
            plot(linspace(0,INT,INT/DT),Ca_mean[:]*1e6,'k')
            xlabel('Time (s)')
            ylabel('Calcium (uM)')
            xlim(1190, 1260)
            savefig(datapath +'/'+ str(idx) +'/Ca_mean.png')
            close()        

    if 'AA' in data_record:
    
        AA_mean = mean(AA, axis=0)
        AAPKC_mean = mean(AAPKC, axis=0)
        AACa1PKC_mean = mean(AACa1PKC, axis=0)
        AACa3PKC_mean = mean(AACa3PKC, axis=0)
        
        AAA_mean = AA_mean+AAPKC_mean+AACa1PKC_mean+AACa3PKC_mean

        ofile_aa=open(datapath +'/'+ str(idx) +'/AA_data.txt','w')
        pickle.dump(AAA_mean, ofile_aa)
        ofile_aa.close()
 
        if plotfigs:                        
            plot(linspace(0,INT,INT/DT),AAA_mean[:],label = 'AA')
            xlabel('Time (s)')
            ylabel('Arachidonic acid #')
            savefig(datapath +'/'+ str(idx) +'/AA_mean.png')
            close()  

    if 'PKC' in data_record:
            
        PKCstar_mean = mean(PKCstar, axis=0)
        PKCstar2_mean = mean(PKCstar2, axis=0)
        PKCstar3_mean = mean(PKCstar3, axis=0)
        PKCstar4_mean = mean(PKCstar4, axis=0)
        RafRKIPPKCstar_mean = mean(RafRKIPPKCstar, axis=0)
        RafRKIPPKCstar2_mean = mean(RafRKIPPKCstar2, axis=0)
        RafRKIPPKCstar3_mean = mean(RafRKIPPKCstar3, axis=0)
        RafRKIPPKCstar4_mean = mean(RafRKIPPKCstar4, axis=0)
        PKCstarAMPAR_mean = mean(PKCstarAMPAR, axis=0)
        PKCstar2AMPAR_mean = mean(PKCstar2AMPAR, axis=0)
        PKCstar3AMPAR_mean = mean(PKCstar3AMPAR, axis=0)
        PKCstar4AMPAR_mean = mean(PKCstar4AMPAR, axis=0)
        PKCstarGRIPAMPAR_mean = mean(PKCstarGRIPAMPAR, axis=0)
        PKCstar2GRIPAMPAR_mean = mean(PKCstar2GRIPAMPAR, axis=0)
        PKCstar3GRIPAMPAR_mean = mean(PKCstar3GRIPAMPAR, axis=0)
        PKCstar4GRIPAMPAR_mean = mean(PKCstar4GRIPAMPAR, axis=0)
        
        MEKRKIPPKCstar_mean = mean(MEKRKIPPKCstar, axis=0)
        MEKRKIPPKCstar2_mean = mean(MEKRKIPPKCstar2, axis=0)
        MEKRKIPPKCstar3_mean = mean(MEKRKIPPKCstar3, axis=0)
        MEKRKIPPKCstar4_mean = mean(MEKRKIPPKCstar4, axis=0)
                
        PKC_mean = PKCstar_mean + PKCstar2_mean + PKCstar3_mean + PKCstar4_mean + \
                   RafRKIPPKCstar_mean + RafRKIPPKCstar2_mean + RafRKIPPKCstar3_mean + \
                   RafRKIPPKCstar4_mean + PKCstarAMPAR_mean + PKCstar2AMPAR_mean + \
                   PKCstar3AMPAR_mean + PKCstar4AMPAR_mean + \
                   MEKRKIPPKCstar_mean + MEKRKIPPKCstar2_mean + MEKRKIPPKCstar3_mean + \
                   MEKRKIPPKCstar4_mean + PKCstarGRIPAMPAR_mean + PKCstar2GRIPAMPAR_mean + \
                   PKCstar3GRIPAMPAR_mean + PKCstar4GRIPAMPAR_mean
                
        ofile_pkc=open(datapath +'/'+ str(idx) +'/PKC_data.txt','w')
        pickle.dump(PKC_mean, ofile_pkc)
        ofile_pkc.close()     
        
        if plotfigs:
            plot(linspace(0,INT,INT/DT), PKC_mean, 'r')
            ylim(ymin=0.0)
            xlabel('Time (s)')
            ylabel('PKC (#)')
            savefig(datapath +'/'+ str(idx) +'/PKC.png')
            close()
            
    if 'ERK' in data_record:
                    
        ERKstar_mean = mean(ERKstar, axis=0)
        MKPERKstar_mean = mean(MKPERKstar, axis=0)
        ERKstarPLA2_mean = mean(ERKstarPLA2, axis=0)
        ERKstarCa1PLA2_mean = mean(ERKstarCa1PLA2, axis=0)
        ERKstarCa2PLA2_mean = mean(ERKstarCa2PLA2, axis=0)
        
        ERKp_mean = mean(ERKp, axis=0)
        MKPERKp_mean = mean(MKPERKp, axis=0)
        MEKstarERKp_mean = mean(MEKstarERKp, axis=0)
        
        ERKstar_total = ERKstar_mean + MKPERKstar_mean + ERKstarPLA2_mean + ERKstarCa1PLA2_mean + ERKstarCa2PLA2_mean
        
        ERKp_total = ERKp_mean + MKPERKp_mean + MEKstarERKp_mean
        
        ERK_total = ERKstar_mean + MKPERKstar_mean + ERKstarPLA2_mean + ERKstarCa1PLA2_mean + ERKstarCa2PLA2_mean + \
                    ERKp_mean + MKPERKp_mean + MEKstarERKp_mean
        
        ofile_erktotal=open(datapath +'/'+ str(idx) +'/ERKtotal_data.txt','w')
        pickle.dump(ERK_total, ofile_erktotal)
        ofile_erktotal.close()  
        
        
        if plotfigs:
            plot(linspace(0,INT,INT/DT), ERK_total, 'k', label ='total')
            ylim(ymin=0.0)
            xlabel('Time (s)')
            ylabel('ERK (#)')
            savefig(datapath +'/'+ str(idx) +'/ERK.png')
            close()
    
                    
    if 'MEK' in data_record:
    
        MEKp_mean= mean(MEKp, axis=0)
        PP2AMEKp_mean= mean(PP2AMEKp, axis=0)
        RafstarMEKp_mean= mean(RafstarMEKp, axis=0)
        
        MEKstar_mean= mean(MEKstar, axis=0)
        MEKstarERK_mean= mean(MEKstarERK, axis=0)
        MEKstarERKp_mean= mean(MEKstarERKp, axis=0)
        PP2AMEKstar_mean = mean(PP2AMEKstar, axis=0)
        
        MEK_total = MEKp_mean + PP2AMEKp_mean + RafstarMEKp_mean + \
                   MEKstar_mean + MEKstarERK_mean + MEKstarERKp_mean + PP2AMEKstar_mean
                   
        MEKstar_total = MEKstar_mean + MEKstarERK_mean + MEKstarERKp_mean + PP2AMEKstar_mean
        
        MEKp_total = MEKp_mean + PP2AMEKp_mean + RafstarMEKp_mean

        ofile_mektotal=open(datapath +'/'+ str(idx) +'/MEKtotal_data.txt','w')
        pickle.dump(MEK_total, ofile_mektotal)
        ofile_mektotal.close()  
        
        if plotfigs:
            plot(linspace(0,INT,INT/DT), MEK_total, 'k', label='total')
            
            ylim(ymin=0.0)
            xlabel('Time (s)')
            ylabel('MEK (#)')
            savefig(datapath +'/'+ str(idx) +'/MEK.png')
            close()        
          
    
    if 'PLA2' in data_record:
        
        PLA2star1_mean = mean(PLA2star1, axis=0)
        PLA2star1APC_mean = mean(PLA2star1APC, axis=0)            
        PLA2star2_mean = mean(PLA2star2, axis=0)
        PLA2star2memb_mean = mean(PLA2star2memb, axis=0)
        PLA2star2membAPC_mean = mean(PLA2star2membAPC, axis=0)
        Ca1PLA2star2_mean = mean(Ca1PLA2star2, axis=0)
        Ca2PLA2star2_mean = mean(Ca2PLA2star2, axis=0)
        Ca1PLA2star2memb_mean = mean(Ca1PLA2star2memb, axis=0)
        Ca2PLA2star2memb_mean = mean(Ca2PLA2star2memb, axis=0)
        Ca1PLA2star2membAPC_mean = mean(Ca1PLA2star2membAPC, axis=0)
        Ca2PLA2star2membAPC_mean = mean(Ca2PLA2star2membAPC, axis=0)
        PP1PLA2star2_mean = mean(PP1PLA2star2, axis=0)
        PP1Ca1PLA2star2_mean = mean(PP1Ca1PLA2star2, axis=0)
        PP1Ca2PLA2star2_mean = mean(PP1Ca2PLA2star2, axis=0)
        PP2APLA2star2_mean = mean(PP1PLA2star2, axis=0)
        PP2ACa1PLA2star2_mean = mean(PP2ACa1PLA2star2, axis=0)
        PP2ACa2PLA2star2_mean = mean(PP2ACa2PLA2star2, axis=0)
                    
        PLA2_mean = PLA2star1_mean + PLA2star1APC_mean + PLA2star2_mean + \
                    PLA2star2memb_mean + PLA2star2membAPC_mean + \
                    Ca1PLA2star2_mean + Ca2PLA2star2_mean + Ca1PLA2star2memb_mean + \
                    Ca2PLA2star2memb_mean + Ca1PLA2star2membAPC_mean + Ca2PLA2star2membAPC_mean + \
                    PP1PLA2star2_mean + PP1Ca1PLA2star2_mean + PP1Ca2PLA2star2_mean + \
                    PP2APLA2star2_mean + PP2ACa1PLA2star2_mean + PP2ACa2PLA2star2_mean        
            
        ofile_pla2=open(datapath +'/'+ str(idx) +'/PLA2_data.txt','w')
        pickle.dump(PLA2_mean, ofile_pla2)
        ofile_pla2.close()             
    
        if plotfigs:
            plot(linspace(0,INT,INT/DT), PLA2_mean, 'y')
            ylim(ymin=0.0)
            xlabel('Time (s)')
            ylabel('cPLA2 (#)')
            savefig(datapath +'/'+ str(idx) +'/PLA2.png')
            close()    
    
    
    if 'Raf' in data_record:

        Rafstar_mean= mean(Rafstar, axis=0)
        PP5Rafstar_mean= mean(PP5Rafstar, axis=0)
        RafstarMEK_mean= mean(RafstarMEK, axis=0)
        RafstarMEKp_mean= mean(RafstarMEKp, axis=0)            
        
        Raf_mean = Rafstar_mean + PP5Rafstar_mean + RafstarMEK_mean + RafstarMEKp_mean;
        
        ofile_raf=open(datapath +'/'+ str(idx) +'/Raf_data.txt','w')
        pickle.dump(Raf_mean, ofile_raf)
        ofile_raf.close()                     
    
        if plotfigs:
            plot(linspace(0,INT,INT/DT), Raf_mean, 'm')
            ylim(ymin=0.0)
            xlabel('Time (s)')
            ylabel('Raf (#)')
            savefig(datapath +'/'+ str(idx) +'/Raf.png')
            close()            

    if 'RKIP' in data_record:

        RKIPstar_mean= mean(RKIPstar, axis=0)
        RKIPstarRP_mean= mean(RKIPstarRP, axis=0)

        RKIP_mean = RKIPstar_mean + RKIPstarRP_mean

        ofile_rkip=open(datapath +'/'+ str(idx) +'/RKIP_data.txt','w')
        pickle.dump(RKIP_mean, ofile_rkip)
        ofile_rkip.close()

        if plotfigs:
            plot(linspace(0,INT,INT/DT), RKIP_mean, 'r')
            ylim(ymin=0.0)
            xlabel('Time (s)')
            ylabel('RKIP (#)')
            savefig(datapath +'/'+ str(idx) +'/RKIP.png')
            close()


run_simulation()