#############################################################################
#                                                                           #
#                Model of Cerebellar Long-Term Depression                   #
#                                                                           #
#                     developed by Gabriela Antunes                         #  
#          updated for STEPS 1.3 by Iain Hepburn and Cory Simon             #
#                                                                           #
#############################################################################


# WARNING: THIS SCRIPT RUNS ON STEPS 1.2 and higher.
# Please use file ltd.py if you are using earlier versions of STEPS


# The reactions described in this script have the same identities and are presented in the same order
# as the reactions described in Supplementary Table I.
# Note, however, that some species might be represented by slightly different names.

import datetime
import steps.model as smodel
import math
import sys
import numpy 
import numpy as np
from numpy import mean
import steps.solver as swmdirect
import steps.geom as swm
import steps.rng as srng
from pylab import*

DT = 0.05 # time step
INT = 3600 # final time (s) simulated
NITER = 2 # number of runs used to calculate the mean
tinit = 600 # initial time that will be discharged. The initial 10 minutes of all simulations should be discharged to allow the model to reach equilibrium.

# Avogadro constant
Na = 6.02214129e23

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

# Ca2+ pulse. IMPORTANT: the duration and amplitude of the pulse should be checked directly in the Ca2+ concentration
# The changes in [Ca2+] generated by the pulses are highly stochastic
# altering max_molar alters the maximum molar concentration per second injection (M/s) 
# cntr defines the time the pulse is at maximum (s)
# FW100M defines the full-width 100th maximum of the Gaussian (s)

def gaussian_ica(t,max_molar=40.0e-6/0.05,cntr=1200.0,FW100M=3):
    c = FW100M/(2*math.sqrt(2*math.log(100)))
    return max_molar*np.exp(-(((t-cntr)**2)/(2*c**2)))

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

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) 
    Rafact = smodel.Spec('Rafact', mdl)     
    PKCstarRafact = smodel.Spec('PKCstarRafact', mdl)      
    PKCstar2Rafact = smodel.Spec('PKCstar2Rafact', mdl)    
    PKCstar3Rafact = smodel.Spec('PKCstar3Rafact', mdl)
    PKCstar4Rafact = smodel.Spec('PKCstar4Rafact', mdl)   
    Rafactstar = smodel.Spec('Rafactstar', mdl)    
    Raf = smodel.Spec('Raf', mdl)  
    RafactstarRaf = smodel.Spec('RafactstarRaf', 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)         
    AMPAR_P = smodel.Spec('AMPAR_P', mdl)                  
    PP2AAMPAR_P = smodel.Spec('PP2AAMPAR_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)
    
    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 = 25000e6
    Reac2 = smodel.SReac('Reac2', s, slhs = [Ca1PMCA], irhs = [Ca], srhs = [PMCA])
    Reac2.kcst = 2000 
    Reac3 = smodel.SReac('Reac3', s, slhs = [Ca1PMCA], srhs = [PMCA])	
    Reac3.kcst = 500
    
    # Ca + NCX <->  Ca1NCX + Ca <->  Ca2NCX -> NCX 
    Reac4 = smodel.SReac('Reac4', s, ilhs = [Ca], slhs = [NCX], srhs = [Ca1NCX])
    Reac4.kcst = 93.827e6  
    Reac5 = smodel.SReac('Reac5', s, slhs = [Ca1NCX], irhs = [Ca], srhs = [NCX])
    Reac5.kcst = 612.6
    Reac6 = smodel.SReac('Reac6', s, ilhs = [Ca], slhs = [Ca1NCX], srhs = [Ca2NCX])
    Reac6.kcst = 93.827e6 
    Reac7 = smodel.SReac('Reac7', s, slhs = [Ca2NCX], irhs = [Ca], srhs = [Ca1NCX])
    Reac7.kcst = 612.6
    Reac8 = smodel.SReac('Reac8', s, slhs = [Ca2NCX], srhs = [NCX])
    Reac8.kcst = 1000
    
    # Ca + SERCA <->  Ca1SERCA +Ca <->  Ca2SERCA  ->  SERCA
    Reac9 = smodel.SReac('Reac9', ERs, olhs = [Ca], slhs = [SERCA], srhs = [Ca1SERCA])
    Reac9.kcst = 17147e6 
    Reac10 = smodel.SReac('Reac10', ERs, slhs = [Ca1SERCA], orhs = [Ca], srhs = [SERCA])
    Reac10.kcst = 8426.3
    Reac11 = smodel.SReac('Reac11', ERs, olhs = [Ca], slhs = [Ca1SERCA], srhs = [Ca2SERCA])
    Reac11.kcst = 17147e6 
    Reac12 = smodel.SReac('Reac12', ERs, slhs = [Ca2SERCA], orhs = [Ca], srhs = [Ca1SERCA])
    Reac12.kcst = 8426.3
    Reac13 = smodel.SReac('Reac13', ERs, slhs = [Ca2SERCA], srhs = [SERCA], irhs = [Ca,Ca])
    Reac13.kcst = 250
    
    # Leak
    Reac14 = smodel.Reac('Reac14', vsys, rhs = [Ca])	
    # reaction constant for zero-order reaction units M/s in STEPS 1.2 and above
    Reac14.kcst = 1900/(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* + Rafact <-> Ca3PKC*.Rafact -> Ca3PKC* + Rafact* (Rafactstar)
    Reac59 = smodel.SReac('Reac59', s, ilhs = [Rafact], slhs = [PKCstar], srhs = [PKCstarRafact])
    Reac59.kcst = 5.80e6
    Reac60 = smodel.SReac('Reac60', s, slhs = [PKCstarRafact], irhs = [Rafact], srhs = [PKCstar])
    Reac60.kcst = 3.608
    Reac61 = smodel.SReac('Reac61', s, slhs = [PKCstarRafact], srhs = [PKCstar], irhs = [Rafactstar])
    Reac61.kcst = 4.7
    
    
    
    # AAPKC* + Rafact <-> AAPKC*.Rafact -> AAPKC* + Rafact* (Rafactstar)
    Reac62 = smodel.SReac('Reac62', s, ilhs = [Rafact], slhs = [PKCstar2], srhs = [PKCstar2Rafact])
    Reac62.kcst =  5.80e6  
    Reac63 = smodel.SReac('Reac63', s, slhs = [PKCstar2Rafact], irhs = [Rafact], srhs = [PKCstar2])
    Reac63.kcst = 3.608
    Reac64 = smodel.SReac('Reac64', s, slhs = [PKCstar2Rafact], srhs = [PKCstar2], irhs = [Rafactstar])
    Reac64.kcst =  4.7
    
    
    
    # AACa1PKC* + Rafact <-> AACa1PKC*.Rafact -> AACa1PKC* + Rafact* (Rafactstar)
    Reac65 = smodel.SReac('Reac65', s, ilhs = [Rafact], slhs = [PKCstar3], srhs = [PKCstar3Rafact])
    Reac65.kcst =  5.80e6  
    Reac66 = smodel.SReac('Reac66', s, slhs = [PKCstar3Rafact], irhs = [Rafact], srhs = [PKCstar3])
    Reac66.kcst = 3.608
    Reac67 = smodel.SReac('Reac67', s, slhs = [PKCstar3Rafact], srhs = [PKCstar3], irhs = [Rafactstar])
    Reac67.kcst =  4.7
    
    
    
    # AACa3PKC* + Rafact <-> AACa3PKC*.Rafact -> AACa3PKC* + Rafact* (Rafactstar)
    Reac68 = smodel.SReac('Reac68', s, ilhs = [Rafact], slhs = [PKCstar4], srhs = [PKCstar4Rafact])
    Reac68.kcst =  5.80e6  
    Reac69 = smodel.SReac('Reac69', s, slhs = [PKCstar4Rafact], irhs = [Rafact], srhs = [PKCstar4])
    Reac69.kcst = 3.608
    Reac70 = smodel.SReac('Reac70', s, slhs = [PKCstar4Rafact], srhs = [PKCstar4], irhs =  [Rafactstar])
    Reac70.kcst = 4.7
    
    
    
    # Rafact* ->  Rafact, deactivation 
    Reac71 = smodel.Reac('Reac71', vsys, lhs = [Rafactstar], rhs = [Rafact])
    Reac71.kcst = 1.
    
    
    
    # Rafact* + Raf <-> Rafact*. Raf -> Rafact* + Raf* (Rafstar)
    Reac72 = smodel.Reac('Reac72', vsys, lhs = [Rafactstar, Raf], rhs = [RafactstarRaf])
    Reac72.kcst = 1e6 
    Reac73 = smodel.Reac('Reac73', vsys, lhs = [RafactstarRaf], rhs = [Rafactstar, Raf])
    Reac73.kcst = 2
    Reac74 = smodel.Reac('Reac74', vsys, lhs = [RafactstarRaf], rhs = [Rafactstar, Rafstar])
    Reac74.kcst = 1.5
    
    
    
    # PP5 (PP) + Raf* <-> PP5.Raf* -> PP5 + Raf
    Reac75 = smodel.Reac('Reac75', vsys, lhs = [PP5, Rafstar], rhs = [PP5Rafstar])
    Reac75.kcst = .55e6 
    Reac76 = smodel.Reac('Reac76', vsys, lhs = [PP5Rafstar], rhs = [PP5, Rafstar])
    Reac76.kcst = 2
    Reac77 = smodel.Reac('Reac77', vsys, lhs = [PP5Rafstar], rhs = [PP5, Raf])
    Reac77.kcst = 0.5
    
    
    
    # Raf* + MEK <-> Raf*.MEK -> Raf* + MEKp <->Raf*.MEKp -> Raf* + MEK* (MEKstar)
    Reac78 = smodel.Reac('Reac78', vsys, lhs = [Rafstar, MEK], rhs = [RafstarMEK])
    Reac78.kcst = 0.65e6 
    Reac79 = smodel.Reac('Reac79', vsys, lhs = [RafstarMEK], rhs = [Rafstar, MEK])
    Reac79.kcst = 0.065
    Reac80 = smodel.Reac('Reac80', vsys, lhs = [RafstarMEK], rhs = [Rafstar, MEKp])
    Reac80.kcst = 1.0
    Reac81 = smodel.Reac('Reac81', vsys, lhs = [Rafstar, MEKp], rhs = [RafstarMEKp])
    Reac81.kcst = 0.65e6 
    Reac82 = smodel.Reac('Reac82', vsys, lhs = [RafstarMEKp], rhs = [Rafstar, MEKp])
    Reac82.kcst = 0.065
    Reac83 = smodel.Reac('Reac83', vsys, lhs = [RafstarMEKp], rhs = [Rafstar, MEKstar])
    Reac83.kcst = 1.0
    
    
    
    # PP2A + MEK* <-> PP2A.MEK* -> PP2A + MEKp <-> PP2A.MEKp -> PP2A + MEK
    Reac84 = smodel.Reac('Reac84', vsys, lhs = [PP2A, MEKstar], rhs = [PP2AMEKstar])
    Reac84.kcst = 0.75e6 
    Reac85 = smodel.Reac('Reac85', vsys, lhs = [PP2AMEKstar], rhs = [PP2A, MEKstar])
    Reac85.kcst = 2
    Reac86 = smodel.Reac('Reac86', vsys, lhs = [PP2AMEKstar], rhs = [PP2A, MEKp])
    Reac86.kcst = 0.5
    Reac87 = smodel.Reac('Reac87', vsys, lhs = [PP2A, MEKp], rhs = [PP2AMEKp])
    Reac87.kcst = 0.75e6 
    Reac88 = smodel.Reac('Reac88', vsys, lhs = [PP2AMEKp], rhs = [PP2A, MEKp])
    Reac88.kcst = 2
    Reac89 = smodel.Reac('Reac89', vsys, lhs = [PP2AMEKp], rhs = [PP2A, MEK])
    Reac89.kcst = 0.5
    
    
    
    # MEK* + ERK <-> MEK*.ERK -> MEK* + ERKp <-> MEK*.ERKp -> MEK* + ERK* (ERKstar)
    Reac90 = smodel.Reac('Reac90', vsys, lhs = [MEKstar, ERK], rhs = [MEKstarERK])
    Reac90.kcst = 16.2e6 
    Reac91 = smodel.Reac('Reac91', vsys, lhs = [MEKstarERK], rhs = [MEKstar, ERK])
    Reac91.kcst = 0.6
    Reac92 = smodel.Reac('Reac92', vsys, lhs = [MEKstarERK], rhs = [MEKstar, ERKp])
    Reac92.kcst = 0.15
    Reac93 = smodel.Reac('Reac93', vsys, lhs = [MEKstar, ERKp], rhs = [MEKstarERKp])
    Reac93.kcst = 16.2e6 
    Reac94 = smodel.Reac('Reac94', vsys, lhs = [MEKstarERKp], rhs = [MEKstar, ERKp])
    Reac94.kcst = 0.6
    Reac95 = smodel.Reac('Reac95', vsys, lhs = [MEKstarERKp], rhs = [MEKstar, ERKstar])
    Reac95.kcst = 0.3 
    
    
    
    # MKP + ERK* <-> MKP.ERK* -> MKP + ERKp <-> MKP.ERKp -> MKP + ERK
    Reac96 = smodel.Reac('Reac96', vsys, lhs = [MKP, ERKstar], rhs = [MKPERKstar])
    Reac96.kcst = 13e6 
    Reac97 = smodel.Reac('Reac97', vsys, lhs = [MKPERKstar], rhs = [MKP, ERKstar])
    Reac97.kcst = 0.396
    Reac98 = smodel.Reac('Reac98', vsys, lhs = [MKPERKstar], rhs = [MKP, ERKp])
    Reac98.kcst = 0.099 
    Reac99 = smodel.Reac('Reac99', vsys, lhs = [MKP, ERKp], rhs = [MKPERKp])
    Reac99.kcst = 28e6 
    Reac100 = smodel.Reac('Reac100', vsys, lhs = [MKPERKp], rhs = [MKP, ERKp])
    Reac100.kcst = 0.56
    Reac101 = smodel.Reac('Reac101', vsys, lhs = [MKPERKp], rhs = [MKP, ERK])
    Reac101.kcst = 0.14 
    
    
    
    # Ca + PLA2 <-> Ca1PLA2 + Ca <-> Ca2PLA2 <-> Ca2PLA2* (PLA2star1)
    Reac102 = smodel.Reac('Reac102', vsys, lhs = [PLA2, Ca], rhs = [Ca1PLA2])
    Reac102.kcst = 1.93e6
    Reac103 = smodel.Reac('Reac103', vsys, lhs = [Ca1PLA2], rhs = [PLA2, Ca])
    Reac103.kcst = 108
    Reac104 = smodel.Reac('Reac104', vsys, lhs = [Ca, Ca1PLA2], rhs = [Ca2PLA2])
    Reac104.kcst = 10.8e6 
    Reac105 = smodel.Reac('Reac105', vsys, lhs = [Ca2PLA2], rhs = [Ca1PLA2, Ca])
    Reac105.kcst = 108
    Reac106 = smodel.SReac('Reac106', s, ilhs = [Ca2PLA2], srhs = [PLA2star1])
    Reac106.kcst = 300 
    Reac107 = smodel.SReac('Reac107', s, slhs = [PLA2star1], irhs = [Ca2PLA2])
    Reac107.kcst = 15
    
    
    
    # Ca1PLA2 <-> Ca1PLA2memb
    Reac108 = smodel.SReac('Reac108', s, ilhs = [Ca1PLA2],  srhs = [Ca1PLA2memb])
    Reac108.kcst = 30
    Reac109 = smodel.SReac('Reac109', s, slhs = [Ca1PLA2memb], irhs = [Ca1PLA2])
    Reac109.kcst = 15
    
    
    
    # PLA2 <-> PLA2memb 
    Reac110 = smodel.SReac('Reac110', s, ilhs = [PLA2], srhs = [PLA2memb])
    Reac110.kcst = 3
    Reac111 = smodel.SReac('Reac111', s, slhs = [PLA2memb], irhs = [PLA2])
    Reac111.kcst = 15
    
    
    
    # PLA2memb + Ca <-> Ca1PLA2memb + Ca <-> Ca2PLA2* (PLA2star1)
    Reac112 = smodel.SReac('Reac112', s, slhs = [PLA2memb], ilhs = [Ca], srhs = [Ca1PLA2memb])
    Reac112.kcst = 1.93e6
    Reac113 = smodel.SReac('Reac113', s, ilhs = [Ca1PLA2memb], srhs = [PLA2memb], irhs = [Ca])
    Reac113.kcst = 0.41
    Reac114 = smodel.SReac('Reac114', s, slhs = [Ca1PLA2memb], ilhs = [Ca], srhs = [PLA2star1])
    Reac114.kcst = 10.8e6 
    Reac115 = smodel.SReac('Reac115', s, ilhs = [PLA2star1], srhs = [Ca1PLA2memb], irhs = [Ca])
    Reac115.kcst = 2.5
    
    
    
    # PLA2star1 <-> PLA2star1APC -> PLA2star1 + AA
    Reac116 = smodel.SReac('Reac116', s, slhs = [PLA2star1],  srhs = [PLA2star1APC])
    Reac116.kcst = 43  
    Reac117 = smodel.SReac('Reac117', s, slhs = [PLA2star1APC],  srhs = [PLA2star1])
    Reac117.kcst = 600
    Reac118 = smodel.SReac('Reac118', s, slhs = [PLA2star1APC], irhs = [AA], srhs = [PLA2star1])
    Reac118.kcst = 450
    
    
    
    # ERK* + PLA2 <-> ERK*.PLA2 -> ERK* + PLA2p (PLA2star2)
    Reac119 = smodel.Reac('Reac119', vsys, lhs = [ERKstar, PLA2], rhs = [ERKstarPLA2])
    Reac119.kcst = 4e6 
    Reac120 = smodel.Reac('Reac120', vsys, lhs = [ERKstarPLA2], rhs = [ERKstar, PLA2])
    Reac120.kcst = 1
    Reac121 = smodel.Reac('Reac121', vsys, lhs = [ERKstarPLA2], rhs = [ERKstar, PLA2star2])
    Reac121.kcst = 14
    
    
    
    # ERK* + Ca1PLA2 <-> ERK*.Ca1PLA2 -> ERK* + Ca1PLA2p (Ca1PLA2star2)
    Reac122 = smodel.Reac('Reac122', vsys, lhs = [ERKstar, Ca1PLA2], rhs = [ERKstarCa1PLA2])
    Reac122.kcst = 4e6 
    Reac123 = smodel.Reac('Reac123', vsys, lhs = [ERKstarCa1PLA2], rhs = [ERKstar, Ca1PLA2])
    Reac123.kcst = 1
    Reac124 = smodel.Reac('Reac124', vsys, lhs = [ERKstarCa1PLA2], rhs = [ERKstar, Ca1PLA2star2])
    Reac124.kcst = 14
    
    
    
    # ERK* + Ca2PLA2 <-> ERK*.Ca2PLA2 -> ERK* + Ca2PLA2p (Ca2PLA2star2)
    Reac125 = smodel.Reac('Reac125', vsys, lhs = [ERKstar, Ca2PLA2], rhs = [ERKstarCa2PLA2])
    Reac125.kcst = 4e6 
    Reac126 = smodel.Reac('Reac126', vsys, lhs = [ERKstarCa2PLA2], rhs = [ERKstar, Ca2PLA2])
    Reac126.kcst = 1
    Reac127 = smodel.Reac('Reac127', vsys, lhs = [ERKstarCa2PLA2], rhs = [ERKstar, Ca2PLA2star2])
    Reac127.kcst = 14
    
    
    # PP2A + PLA2p <-> PP2A.PLA2p -> PP2A + PLA2
    Reac128 = smodel.Reac('Reac128', vsys, lhs = [PP2A, PLA2star2], rhs = [PP2APLA2star2])
    Reac128.kcst = 1.4e6
    Reac129 = smodel.Reac('Reac129', vsys, lhs = [PP2APLA2star2], rhs = [PP2A, PLA2star2])
    Reac129.kcst = 1.5
    Reac130 = smodel.Reac('Reac130', vsys, lhs = [PP2APLA2star2], rhs = [PLA2, PP2A])
    Reac130.kcst = 2.5
    
    
    
    # PP2A + Ca1PLA2p <-> PP2A.Ca1PLA2p -> PP2A + Ca1PLA2
    Reac131 = smodel.Reac('Reac131', vsys, lhs = [PP2A, Ca1PLA2star2], rhs = [PP2ACa1PLA2star2])
    Reac131.kcst = 1.4e6
    Reac132 = smodel.Reac('Reac132', vsys, lhs = [PP2ACa1PLA2star2], rhs = [PP2A, Ca1PLA2star2])
    Reac132.kcst = 1.5 
    Reac133 = smodel.Reac('Reac133', vsys, lhs = [PP2ACa1PLA2star2], rhs = [PP2A, Ca1PLA2])
    Reac133.kcst = 2.5
    
    
    
    # PP2A + Ca2PLA2p <-> PP2A.Ca2PLA2p -> PP2A + Ca2PLA2
    Reac134 = smodel.Reac('Reac134', vsys, lhs = [PP2A, Ca2PLA2star2], rhs = [PP2ACa2PLA2star2])
    Reac134.kcst = 1.4e6
    Reac135 = smodel.Reac('Reac135', vsys, lhs = [PP2ACa2PLA2star2], rhs = [PP2A, Ca2PLA2star2])
    Reac135.kcst = 1.5 
    Reac136 = smodel.Reac('Reac136', vsys, lhs = [PP2ACa2PLA2star2], rhs = [PP2A, Ca2PLA2])
    Reac136.kcst = 2.5    
    
    
    
    # PP1 + PLA2p <-> PP1.PLA2p -> PP1 + PLA2
    Reac137 = smodel.Reac('Reac137', vsys, lhs = [PP1, PLA2star2], rhs = [PP1PLA2star2])
    Reac137.kcst = 1.4e6
    Reac138 = smodel.Reac('Reac138', vsys, lhs = [PP1PLA2star2], rhs = [PP1, PLA2star2])
    Reac138.kcst = 1.5
    Reac139 = smodel.Reac('Reac139', vsys, lhs = [PP1PLA2star2], rhs = [PLA2, PP1])
    Reac139.kcst = 2.5
    
    
    
    # PP1 + Ca1PLA2p <-> PP1.Ca1PLA2p -> PP1 + Ca1PLA2
    Reac140 = smodel.Reac('Reac140', vsys, lhs = [PP1, Ca1PLA2star2], rhs = [PP1Ca1PLA2star2])
    Reac140.kcst = 1.4e6
    Reac141 = smodel.Reac('Reac141', vsys, lhs = [PP1Ca1PLA2star2], rhs = [PP1, Ca1PLA2star2])
    Reac141.kcst = 1.5
    Reac142 = smodel.Reac('Reac142', vsys, lhs = [PP1Ca1PLA2star2], rhs = [PP1, Ca1PLA2])
    Reac142.kcst = 2.5
    
    
    
    # PP1 + Ca2PLA2p <-> PP1.Ca2PLA2p -> PP1 + Ca2PLA2
    Reac143 = smodel.Reac('Reac143', vsys, lhs = [PP1, Ca2PLA2star2], rhs = [PP1Ca2PLA2star2])
    Reac143.kcst = 1.4e6
    Reac144 = smodel.Reac('Reac144', vsys, lhs = [PP1Ca2PLA2star2], rhs = [PP1, Ca2PLA2star2])
    Reac144.kcst = 1.5
    Reac145 = smodel.Reac('Reac145', vsys, lhs = [PP1Ca2PLA2star2], rhs = [PP1, Ca2PLA2])
    Reac145.kcst = 2.5
    
    
    
    # PLA2p <-> PLA2** (PLA2star2memb) <-> PLA2**APC -> PLA2** + AA 
    Reac146 = smodel.SReac('Reac146', s, ilhs = [PLA2star2],  srhs = [PLA2star2memb])
    Reac146.kcst = 50
    Reac147 = smodel.SReac('Reac147', s, slhs = [PLA2star2memb], irhs = [PLA2star2])
    Reac147.kcst = 15
    Reac148 = smodel.SReac('Reac148', s, slhs = [PLA2star2memb],  srhs = [PLA2star2membAPC])
    Reac148.kcst = 43 
    Reac149 = smodel.SReac('Reac149', s, slhs = [PLA2star2membAPC],  srhs = [PLA2star2memb])
    Reac149.kcst = 600
    Reac150 = smodel.SReac('Reac150', s, slhs = [PLA2star2membAPC], srhs = [PLA2star2memb], irhs = [AA])
    Reac150.kcst = 3600 
    
    
    
    # Ca1PLA2p <-> Ca1PLA2** (Ca1PLA2star2memb) <-> Ca1PLA2**APC -> Ca1PLA2** + AA
    Reac151 = smodel.SReac('Reac151', s, ilhs = [Ca1PLA2star2], srhs = [Ca1PLA2star2memb])
    Reac151.kcst = 30
    Reac152 = smodel.SReac('Reac152', s, slhs = [Ca1PLA2star2memb], irhs = [Ca1PLA2star2])
    Reac152.kcst = 15
    Reac153 = smodel.SReac('Reac153', s, slhs = [Ca1PLA2star2memb],  srhs = [Ca1PLA2star2membAPC])
    Reac153.kcst = 43 
    Reac154 = smodel.SReac('Reac154', s, slhs = [Ca1PLA2star2membAPC],  srhs = [Ca1PLA2star2memb])
    Reac154.kcst = 600
    Reac155 = smodel.SReac('Reac155', s, slhs = [Ca1PLA2star2membAPC], srhs = [Ca1PLA2star2memb], irhs = [AA])
    Reac155.kcst = 3600
    
    
    
    # Ca2PLA2p <-> Ca2PLA2** (Ca2PLA2star2memb) <-> Ca2PLA2**APC -> Ca2PLA2** + AA
    Reac156 = smodel.SReac('Reac156', s, ilhs = [Ca2PLA2star2], srhs = [Ca2PLA2star2memb])
    Reac156.kcst = 300
    Reac157 = smodel.SReac('Reac157', s, slhs = [Ca2PLA2star2memb], irhs = [Ca2PLA2star2])
    Reac157.kcst = 15   
    Reac158 = smodel.SReac('Reac158', s, slhs = [Ca2PLA2star2memb],  srhs = [Ca2PLA2star2membAPC])
    Reac158.kcst = 43 
    Reac159 = smodel.SReac('Reac159', s, slhs = [Ca2PLA2star2membAPC],  srhs = [Ca2PLA2star2memb])
    Reac159.kcst = 600
    Reac160 = smodel.SReac('Reac160', s, slhs = [Ca2PLA2star2membAPC], srhs = [Ca2PLA2star2memb], irhs = [AA])
    Reac160.kcst = 3600
    
    
    
    # Ca + PLA2p <-> Ca1PLA2p
    Reac161 = smodel.Reac('Reac161', vsys, lhs = [PLA2star2, Ca], rhs = [Ca1PLA2star2])
    Reac161.kcst = 1.93e6 
    Reac162 = smodel.Reac('Reac162', vsys, lhs = [Ca1PLA2star2], rhs = [PLA2star2, Ca])
    Reac162.kcst = 108
    
    
    
    # Ca + PLA2** <-> Ca1PLA2**
    Reac163 = smodel.SReac('Reac163', s, slhs = [PLA2star2memb], ilhs = [Ca], srhs = [Ca1PLA2star2memb])
    Reac163.kcst = 1.93e6 
    Reac164 = smodel.SReac('Reac164', s, ilhs = [Ca1PLA2star2memb], srhs = [PLA2star2memb], irhs = [Ca])
    Reac164.kcst = 0.41
    
    
    
    # Ca + Ca1PLA2p <-> Ca2PLA2p
    Reac165 = smodel.Reac('Reac165', vsys, lhs = [Ca, Ca1PLA2star2], rhs = [Ca2PLA2star2])
    Reac165.kcst = 10.8e6 
    Reac166 = smodel.Reac('Reac166', vsys, lhs = [Ca2PLA2star2], rhs = [Ca1PLA2star2, Ca])
    Reac166.kcst = 108
    
    
    
    # Ca + Ca1PLA2** <-> Ca2PLA2**
    Reac167 = smodel.SReac('Reac167', s, slhs = [Ca1PLA2star2memb], ilhs = [Ca], srhs = [Ca2PLA2star2memb])
    Reac167.kcst = 10.8e6
    Reac168 = smodel.SReac('Reac168', s, ilhs = [Ca2PLA2star2memb], srhs = [Ca1PLA2star2memb], irhs = [Ca])
    Reac168.kcst = 2.5
    
    
    
    # AA -> 0 ,degradation
    Reac169 = smodel.SReac('Reac169', s, ilhs = [AA])
    Reac169.kcst = .4
    
    
    
    # Ca3PKC* + AMPAR <-> Ca3PKC*.AMPAR -> Ca3PKC* + AMPARp 
    Reac170 = smodel.SReac('Reac170', 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
    Reac170.kcst = 0.4e6*(1.02e-12/1.87e-19)
    Reac171 = smodel.SReac('Reac171', s, slhs = [PKCstarAMPAR], srhs = [PKCstar, AMPAR])
    Reac171.kcst = 0.8
    Reac172 = smodel.SReac('Reac172', s, slhs = [PKCstarAMPAR], srhs = [PKCstar, AMPAR_P])
    Reac172.kcst = 0.3
    
    
    
    # AAPKC* + AMPAR <-> AAPKC*.AMPAR -> AAPKC* + AMPARp
    Reac173 = smodel.SReac('Reac173', 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
    Reac173.kcst = 0.4e6*(1.02e-12/1.87e-19)
    Reac174 = smodel.SReac('Reac174', s, slhs = [PKCstar2AMPAR], srhs = [PKCstar2, AMPAR])
    Reac174.kcst = 0.8
    Reac175 = smodel.SReac('Reac175', s, slhs = [PKCstar2AMPAR], srhs = [PKCstar2, AMPAR_P])
    Reac175.kcst = 0.3
    
    
    
    # AACa1PKC* + AMPAR <-> AACa1PKC*.AMPAR -> AACa1PKC* + AMPARp
    Reac176 = smodel.SReac('Reac176', 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
    Reac176.kcst = 0.4e6*(1.02e-12/1.87e-19)   
    Reac177 = smodel.SReac('Reac177', s, slhs = [PKCstar3AMPAR], srhs = [PKCstar3, AMPAR])
    Reac177.kcst = 0.8
    Reac178 = smodel.SReac('Reac178', s, slhs = [PKCstar3AMPAR], srhs = [PKCstar3, AMPAR_P])
    Reac178.kcst =  0.3
    
    
    
    # AACa3PKC* + AMPAR <-> AACa3PKC*.AMPAR -> AACa3PKC* + AMPARp 
    Reac179 = smodel.SReac('Reac179', 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
    Reac179.kcst = 0.4e6*(1.02e-12/1.87e-19)  
    Reac180 = smodel.SReac('Reac180', s, slhs = [PKCstar4AMPAR], srhs = [PKCstar4, AMPAR])
    Reac180.kcst = 0.8 
    Reac181 = smodel.SReac('Reac181', s, slhs = [PKCstar4AMPAR], srhs = [PKCstar4, AMPAR_P])
    Reac181.kcst = 0.3
    
    
    
    # PP2A + AMPARp <-> PP2A.AMPARp -> PP2A + AMPAR
    Reac182 = smodel.SReac('Reac182', s, ilhs = [PP2A], slhs = [AMPAR_P], srhs = [PP2AAMPAR_P])
    Reac182.kcst = 0.6e6
    Reac183 = smodel.SReac('Reac183', s, slhs = [PP2AAMPAR_P], irhs = [PP2A], srhs = [AMPAR_P])
    Reac183.kcst = 0.17
    Reac184 = smodel.SReac('Reac184', s, slhs = [PP2AAMPAR_P], srhs = [AMPAR], irhs = [PP2A])
    Reac184.kcst = 0.25 
    
    
    
    # AMPAR + GRIP <-> GRIP.AMPA
    Reac185 = smodel.SReac('Reac185', 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
    Reac185.kcst = 1e6*(1.02e-12/1.87e-19)
    Reac186 = smodel.SReac('Reac186', s, slhs = [GRIPAMPAR], srhs = [GRIP, AMPAR])
    Reac186.kcst = 7
    
    
    
    # AMPARp + GRIP <-> GRIP.AMPAp (AMPAR = synaptic AMPAR)
    Reac187 = smodel.SReac('Reac187', 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
    Reac187.kcst = 1e6*(1.02e-12/1.87e-19)
    Reac188 = smodel.SReac('Reac188', s, slhs = [GRIPAMPAR_P], srhs = [GRIP, AMPAR_P])
    Reac188.kcst = 70
    
    
    
    # AMPAR <-> AMPARextra (AMPARextra = extra-synaptic AMPAR)  
    Reac189 = smodel.SReac('Reac189', s, slhs = [AMPAR], srhs = [AMPARextra])
    Reac189.kcst = 0.1
    Reac190 = smodel.SReac('Reac190', s, slhs = [AMPARextra], srhs = [AMPAR])
    Reac190.kcst = 0.02
    
    
    
    # AMPARp <-> AMPARextra_P
    Reac191 = smodel.SReac('Reac191', s, slhs = [AMPAR_P], srhs = [AMPARextra_P])
    Reac191.kcst =  0.1
    Reac192 = smodel.SReac('Reac192', s, slhs = [AMPARextra_P], srhs = [AMPAR_P])
    Reac192.kcst = 0.02
    
    
    
    # PP2A + AMPARextra_P  <-> PP2A.AMPARextra_P  -> PP2A + AMPARextra
    Reac193 = smodel.SReac('Reac193', s, ilhs = [PP2A], slhs = [AMPARextra_P], srhs = [PP2AAMPARextra_P])
    Reac193.kcst = 0.6e6
    Reac194 = smodel.SReac('Reac194', s, slhs = [PP2AAMPARextra_P], irhs = [PP2A], srhs = [AMPARextra_P])
    Reac194.kcst = 0.17
    Reac195 = smodel.SReac('Reac195', s, slhs = [PP2AAMPARextra_P], srhs = [AMPARextra], irhs = [PP2A])
    Reac195.kcst = 0.25
    
    
    
    # AMPARextra <-> AMPARdend (AMPARdend = dendritic AMPAR)  
    Reac196 = smodel.SReac('Reac196', s, slhs = [AMPARextra], srhs = [AMPARdend])
    Reac196.kcst = 0.02
    Reac197 = smodel.SReac('Reac197', s, slhs = [AMPARdend], srhs = [AMPARextra])
    Reac197.kcst = 0.00025
    
    
    
    # AMPARextra_P <-> AMPARdend_P 
    Reac198 = smodel.SReac('Reac198', s, slhs = [AMPARextra_P], srhs = [AMPARdend_P])
    Reac198.kcst =  0.02
    Reac199 = smodel.SReac('Reac199', s, slhs = [AMPARdend_P], srhs = [AMPARextra_P])
    Reac199.kcst = 0.00025
    
    
    
    # PP2A + AMPARp_dend  <-> PP2A.AMPARp_dend  -> PP2A + AMPAR_dend
    Reac200 = smodel.SReac('Reac200', s, ilhs = [PP2A], slhs = [AMPARdend_P], srhs = [PP2AAMPARdend_P])
    Reac200.kcst = 0.6e6
    Reac201 = smodel.SReac('Reac201', s, slhs = [PP2AAMPARdend_P], irhs = [PP2A], srhs = [AMPARdend_P])
    Reac201.kcst = 0.17
    Reac202 = smodel.SReac('Reac202', s, slhs = [PP2AAMPARdend_P], srhs = [AMPARdend], irhs = [PP2A])
    Reac202.kcst = 0.25
    
    
    
    # AMPARdend_P <-> AMPARcyt_P (AMPARcyt = cytosolic AMPAR) 
    Reac203 = smodel.SReac('Reac203', s, slhs = [AMPARdend_P], irhs = [AMPARcyt_P])
    Reac203.kcst =  0.003
    Reac204 = smodel.SReac('Reac204', s, ilhs = [AMPARcyt_P], srhs = [AMPARdend_P])
    Reac204.kcst =  0.002
    
    
    
    # PP2A + AMPARcyt_P  <-> PP2A.AMPARcyt_P  -> PP2A + AMPARcyt 
    Reac205 = smodel.Reac('Reac205', vsys, lhs = [PP2A, AMPARcyt_P], rhs = [PP2AAMPARcyt_P])
    Reac205.kcst = 0.6e6
    Reac206 = smodel.Reac('Reac206', vsys, lhs = [PP2AAMPARcyt_P], rhs = [PP2A, AMPARcyt_P])
    Reac206.kcst = 0.17
    Reac207 = smodel.Reac('Reac207', vsys, lhs = [PP2AAMPARcyt_P], rhs = [AMPARcyt, PP2A])
    Reac207.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_gauss_mean():
    rng=srng.create('mt19937',512)
    m=gen_model() 
    g=gen_geom() 
    tpnts=np.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))
    Rafact = numpy.zeros((NITER, ntpnts))
    PKCstarRafact = numpy.zeros((NITER, ntpnts))
    PKCstar2Rafact = numpy.zeros((NITER, ntpnts))
    PKCstar3Rafact = numpy.zeros((NITER, ntpnts))
    PKCstar4Rafact = numpy.zeros((NITER, ntpnts))
    Rafactstar = numpy.zeros((NITER, ntpnts))
    Raf = numpy.zeros((NITER, ntpnts))
    RafactstarRaf = 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))
    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))

    ica=gaussian_ica(tpnts)
    
    sim = swmdirect.Wmdirect(m, g, rng)
    
    
    for j in arange(NITER): 
        print "Run number {0}".format(j)
        rng.initialize(datetime.datetime.now().microsecond) 
        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', 48) 
        sim.setCompConc('vsys', 'Raf',0.1e-6) 
        sim.setCompConc('vsys', 'Rafact', 0.5e-6)
        sim.setCompConc('vsys', 'MEK',1.5e-6) 
        sim.setCompConc('vsys', 'ERK',1.0e-6)    
        sim.setCompConc('vsys', 'MKP',0.26e-6) 
        sim.setCompConc('vsys', 'PP5', 1.0e-6) 
        sim.setCompConc('vsys', 'PP2A', 1.5e-6) 
        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', 100)
        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', 119) 
        sim.setPatchCount('memb', 'AMPARdend', 1600)
        
        for i in range(ntpnts):
            sim.run(tpnts[i]) 
            print "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')
            Rafact[j,i] = sim.getCompConc('vsys', 'Rafact')
            PKCstarRafact[j,i] = sim.getPatchCount('memb', 'PKCstarRafact')
            PKCstar2Rafact[j,i] = sim.getPatchCount('memb', 'PKCstar2Rafact')
            PKCstar3Rafact[j,i] = sim.getPatchCount('memb', 'PKCstar3Rafact')
            PKCstar4Rafact[j,i] = sim.getPatchCount('memb', 'PKCstar4Rafact')
            Rafactstar[j,i] = sim.getCompCount('vsys', 'Rafactstar')
            RafactstarRaf[j,i] = sim.getCompCount('vsys', 'RafactstarRaf')
            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')
            RafstarMEKp[j,i] = sim.getCompCount('vsys', 'RafstarMEKp')
            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.getCompConc('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')
            Ca1NCX[j,i] = sim.getPatchCount('memb', 'Ca1NCX')
            Ca2NCX[j,i] = sim.getPatchCount('memb', 'Ca2NCX')
            SERCA[j,i] = sim.getPatchCount('ERmemb', 'SERCA')
            Ca1SERCA[j,i] = sim.getPatchCount('ERmemb', 'Ca1SERCA')
            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')   
            
            sim.setCompReacK('vsys', 'Cainflux', ica[i]) # changes reaction constant to be ...
    
####################################################################################

    # Plot Calcium: these lines must be commented if pylab is not installed

    Ca_mean = mean(Ca, axis=0)
                
    plot(linspace(0,INT-tinit,(INT-tinit)/DT),Ca_mean[tinit/DT-1:INT/DT-1]*1e6,'k')
    xlabel('Time (s)')
    ylabel('Calcium (umol.L^-^1)')
    show()
    
    # Plot PKC: these lines must be commented if pylab is not installed
                
    PKCstar_mean = mean(PKCstar, axis=0)
    PKCstar2_mean = mean(PKCstar2, axis=0)
    PKCstar3_mean = mean(PKCstar3, axis=0)
    PKCstar4_mean = mean(PKCstar4, axis=0)
    PKCstarRafact_mean = mean(PKCstarRafact, axis=0)
    PKCstar2Rafact_mean = mean(PKCstar2Rafact, axis=0)
    PKCstar3Rafact_mean = mean(PKCstar3Rafact, axis=0)
    PKCstar4Rafact_mean = mean(PKCstar4Rafact, 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)
            
    PKC_mean=PKCstar_mean+PKCstar2_mean+PKCstar3_mean+PKCstar4_mean+PKCstarRafact_mean+PKCstar2Rafact_mean+PKCstar3Rafact_mean+PKCstar4Rafact_mean+PKCstarAMPAR_mean+PKCstar2AMPAR_mean+PKCstar3AMPAR_mean+PKCstar4AMPAR_mean
    
    plot(linspace(0,INT-tinit,(INT-tinit)/DT), PKC_mean[tinit/DT-1:INT/DT-1], 'r')
    xlabel('Time (s)')
    ylabel('PKC (#)')
    show()
                
                
    # Plot ERK: these lines must be commented if pylab is not installed
                
    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)
                
    ERK_mean=ERKstar_mean+MKPERKstar_mean+ERKstarPLA2_mean+ERKstarCa1PLA2_mean+ERKstarCa2PLA2_mean;
    plot(linspace(0,INT-tinit,(INT-tinit)/DT), ERK_mean[tinit/DT-1:INT/DT-1],'c')
    xlabel('Time (s)')
    ylabel('ERK (#)')
    show()
                
                
    # Plot cPLA2: these lines must be commented if pylab is not installed
        
    PLA2star1_mean = mean(PLA2star1, axis=0)
    PLA2star1APC_mean = mean(PLA2star1APC, axis=0)            
    PLA2star2_mean = mean(PLA2star2, 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+Ca1PLA2star2_mean+Ca2PLA2star2_mean+Ca1PLA2star2memb_mean+Ca2PLA2star2memb_mean+Ca1PLA2star2membAPC_mean+Ca2PLA2star2membAPC_mean+PP1PLA2star2_mean+PP1Ca1PLA2star2_mean+PP1Ca2PLA2star2_mean+PP2APLA2star2_mean+PP2ACa1PLA2star2_mean+PP2ACa2PLA2star2_mean;
    
    plot(linspace(0,INT-tinit,(INT-tinit)/DT), PLA2_mean[tinit/DT-1:INT/DT-1], 'b')
    xlabel('Time (s)')
    ylabel('cPLA2 (#)')
    show()
    
        
    # Plot Synaptic AMPAR: these lines must be commented if pylab is not installed
            
    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)
                
    AMPAR_syn_mean=AMPAR_mean+AMPAR_P_mean+PP2AAMPAR_P_mean+GRIPAMPAR_mean+GRIPAMPAR_P_mean+PKCstarAMPAR_mean+PKCstar2AMPAR_mean+PKCstar3AMPAR_mean+PKCstar4AMPAR_mean;
    plot(linspace(0,INT-tinit,(INT-tinit)/DT), AMPAR_syn_mean[tinit/DT-1:INT/DT-1], 'g')
    xlabel('Time (s)')
    ylabel('Synaptic AMPAR (#)')
    show()
    
    return AMPAR_syn_mean

ampar_syn_mean = run_gauss_mean()