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


# IMPORTANT: this script has not been released. 


# WARNING: THIS SCRIPT RUNS ON STEPS 1.0  and 1.1.
# DO NOT ATTEMPT TO RUN THIS SCRIPT ON ANY OTHER VERSION OF STEPS.
# The model was developed in a pre-released version of STEPS. It has been tested in STEPS 1.0 and 1.1
# and no incompatibilities were found. However, there are incompatibilities with other versions.


# The reactions described in this script have the same identities and are presented in the same order
# of 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
from pylab import* # imports pylab to plot the results, if pylab is not installed this line must be commented


DT = 0.05 # time step for the plots
INT = 3600  # simulated time                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       
NITER = 1 # number of runs that will be 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. 

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

# 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 par1 and par2 alter the duration and magnitude of the pulse respectively, crt defines the moment when the pulse will occur

def gaussian_ica(t, par1 = 110000, ctr = 1200000, par2 = 300):
    t2 = ((t*1.0e3)-ctr)/par2
    return (par1/par2) * np.exp(-t2 * t2)


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

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.SReac('Cainflux', s, irhs=[Ca])
    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])	
    Reac14.kcst = 1900


  
    # 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])
    Reac170.kcst = 0.4e6  
    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])
    Reac173.kcst = 0.4e6  
    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])
    Reac176.kcst = 0.4e6   
    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])
    Reac179.kcst = 0.4e6  
    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])
    Reac185.kcst = 1e6
    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])
    Reac187.kcst = 1e6
    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. The leak current (Reac14) must be scaled by the same ratio.

import steps.geom as swm

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

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


import steps.rng as srng

r = srng.create('mt19937', 512)
r.initialize(datetime.datetime.now().microsecond)

def run_gauss_mean():
    rng = srng.create('mt19937', 512)    

    m = gen_model()
    g = gen_geom()

    sim = swmdirect.Wmdirect(m, g, r)

    tpnts = np.arange(0.0, INT, DT)
    ntpnts = tpnts.shape[0]

    ica = 1.02e-12 * 1.5e14 * gaussian_ica(tpnts)     


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



    for j in xrange(NITER):
        rng.initialize(datetime.datetime.now().microsecond)
        sim.reset()

        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)
        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 xrange(ntpnts):

            sim.run(tpnts[i])

            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.setPatchSReacK('memb', 'Cainflux', ica[i])


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


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

run_gauss_mean()

Loading data, please wait...