Cerebellar long-term depression (LTD) (Antunes and De Schutter 2012)

 Download zip file 
Help downloading and running models
Accession:141270
Many cellular processes involve small number of molecules and undergo stochastic fluctuations in their levels of activity. Among these processes is cerebellar long-term depression (LTD), a form of synaptic plasticity expressed as a reduction in the number of synaptic AMPA receptors (AMPARs) in Purkinje cells. Using a stochastic model of the signaling network and mechanisms of AMPAR trafficking involved in LTD, we show that the network activity in single synapses switches between two discrete stable states (LTD and non-LTD). Stochastic fluctuations affecting more intensely the level of activity of a few components of the network lead to the probabilistic induction of LTD and threshold dithering. The non-uniformly distributed stochasticity of the network allows the stable occurrence of several different macroscopic levels of depression, determining the experimentally observed sigmoidal relationship between the magnitude of depression and the concentration of the triggering signal.
Reference:
1 . Antunes G, De Schutter E (2012) A stochastic signaling network mediates the probabilistic induction of cerebellar long-term depression. J Neurosci 32:9288-300 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Synapse;
Brain Region(s)/Organism:
Cell Type(s): Cerebellum Purkinje GABA cell;
Channel(s):
Gap Junctions:
Receptor(s): AMPA;
Gene(s):
Transmitter(s):
Simulation Environment: STEPS;
Model Concept(s): Synaptic Plasticity;
Implementer(s): De Schutter, Erik [erik at oist.jp]; Antunes, Gabriela [gabri_antunes at hotmail.com];
Search NeuronDB for information about:  Cerebellum Purkinje GABA cell; AMPA;
/
LTDmodel
README.txt
Info_on_the_model.pdf
ltd.py
ltd_updated.py
                            
	
#############################################################################
#                                                                           #
#                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()


Loading data, please wait...