####################################################################################################
# Vesicle cycle model STEPS script                                                                 #
#                                                                                                  #
# Gallimore, A. R., Hepburn, I., Georgiev, S. V., Rizzoli, S. O., & De Schutter, E. (2025)         #
# Dynamic Regulation of Vesicle Pools in a Detailed Spatial Model of the Complete Synaptic Vesicle #
# Cycle. Science Advances                                                                          #
#                                                                                                  #
# Script versions:                                                                                 #
#     Original author: Andrew Gallimore                                                            #
#     Conversion to new STEPS API: Jules Lallouette                                                #
####################################################################################################

import steps.interface

from steps.model import *
from steps.geom import *
from steps.rng import *
from steps.sim import *
from steps.saving import *

import argparse
import collections
import copy
import math
import numpy as np
import os
import pickle
import random
import sys
import time

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

# Example of command line arguments:
# python3 {script.py} . ./meshes/pyr_axon_2021 ./logs ./data SimName 1234 1 1 5 10
# With Blender data being recorded:
# python3 {script.py} . ./meshes/pyr_axon_2021 ./logs ./data SimName 1234 1 1 5 10 true


MODEL_PATH = sys.argv[1] # Path to the folder containing model files (like actin_tets_215nm.txt)
MESH_PATH = sys.argv[2] # Path to the xml mesh file that should be used for the simulation
LOG_PATH = sys.argv[3] # Path to the folder to which logs should be written
DATA_PATH = sys.argv[4] # Path to the folder to which data should be written
DATA_ID = sys.argv[5] # Identifier for the simulation group
JOB_ID = int(sys.argv[6]) # ID of the SLURM job
JOB_INDEX = int(sys.argv[7]) # Index of the SLURM job (start numbering at 1)
kcat_can_tomo = float(sys.argv[8])
RATE = float(sys.argv[9]) # Rate of stimulation in Hz
INT = float(sys.argv[10]) # The simulation endtime (s)
save_blender_data = sys.argv[11].lower() == 'true' if len(sys.argv) > 11 else False
checkpoint_mode = sys.argv[12].lower() if len(sys.argv) > 12 else None # 'checkpoint' to checkpoint, 'restore' to restore, 'both' to do both...
checkpoint_path = sys.argv[13] if len(sys.argv) > 13 else None # Path to the checkpoint file
if checkpoint_mode == 'both':
    checkpoint_path_2 = sys.argv[14] if len(sys.argv) > 14 else None # Path to the second checkpoint file

rng_seed = JOB_ID * JOB_INDEX

print('JOB_ID:', JOB_ID)
print('seed:', rng_seed)

random.seed(1)

exportParam = False

# Only output on rank 0
if MPI.rank == 0:
    if LOG_PATH.lower() != 'none':
        os.makedirs(LOG_PATH, exist_ok=True)
        LOG_FILE = open(os.path.join(LOG_PATH, f'{DATA_ID}_{JOB_INDEX}_{INT}_{rng_seed}_{JOB_ID}_stdout.txt'), 'w')
        sys.stdout = LOG_FILE
    else:
        LOG_FILE = None
else:
    LOG_FILE = None
    sys.stdout = None

print('FINAL VESICLE MODEL VERSION OCT 28 2021, API_2 version NOV 2023')

#Boolean for either static or mobile Ca channels
cachan_static = True
print('CDK5 kcat = 5')
if cachan_static:
    print('CA CHANNELS ARE STATIC')
else:
    print('CA CHANNELS ARE MOBILE')

print('TETHERS ARE ACTIVE')
print('TESTING DOCK SITES 3')

#number of Ca channels
cachan_n = 15

#Choose dock triangle set
dock_tris = [42155, 76431, 84872, 61803, 21332, 71084, 79446, 89381]

#vesicle membrane protein binding
ves_protein_bind_on = True

tomo_init = 1500

can_init = 1500
print('TOTAL CALCINEURIN: ', can_init)

print('kcat_can_tomo: ', kcat_can_tomo)

kcat_can_init = 0.5
print('kcat_can:', kcat_can_init)

#################################### CONSTANTS ####################################

# Number of iterations
NITER = 1

# The data collection time increment (s)
DT = 0.001

# The EField dt
EF_DT = 1.0e-3

########## Ca concentrations (M) #########

Ca_oconc = 2e-3

########## CaP channels density & permiability per channel ##########

# CaP_P is permeability per channel (m3/s)
# CaP_ro is channel/surface area (/m2)
# P in Ca Dynamics model is 0.95e-4 cm/s --> 0.95e-6 m/s

CaP_P = 2.5e-20
print('Ca Chan Perm:', CaP_P)
CaP_ro = 3.8e13

#membrane SNARE proteins diffusiion scaling factor
#######SCALE SNARE MEM DIFF CONSTANTS#######
snare_diff_scale = 1

#Boolean for either presence or absence of Calmodulin
cam_present = True

########## Initial Membrane Potential #########
init_pot = -60e-3

########## BULK RESISTIVITY ##########

Ra = 235.7*1.0e-2

########## MEMBRANE CAPACITANCE ##########

memb_capac = 1.5e-2

##########CaP channel parameters ####################

def minf_cap(mV):
    vhalfm = -29.458
    cvm = 8.429
    return 1 / (1 + math.exp(-(mV - vhalfm) / cvm))

def tau_cap(mV):
    if mV >= -40:
        return 0.2702 + 1.1622 * math.exp(-(mV + 26.798) ** 2 / 164.19)
    else:
        return 0.6923 * math.exp(mV / 1089.372)

alpha_cap = VDepRate.Create(
    lambda V: (minf_cap(V * 1e3) / tau_cap(V * 1e3)) * 1e3
)
beta_cap = VDepRate.Create(
    lambda V: (1 - minf_cap(V * 1e3)) / tau_cap(V * 1e3) * 1e3
)

# Initial conditions
CaP_pinit = [0.92402, 0.073988, 0.0019748, 1.7569e-05]

###Model Definition###

syt_model = Model()
r = ReactionManager()
with syt_model:
    
    vsys, syn_sys, cytERsys = VolumeSystem.Create()
    memsys, ERsys = SurfaceSystem.Create()
    
    #General cytosolic DCST
    DCST_CYT = 2.0e-12
    DCST_CYT_low = 0.028e-12
    # Diffusion constant of Calcium
    DCST_CA = 0.223e-9*1 #CHANGED
    #General membrane DCST
    DCST_MEM = 0.25e-12 #0.05 um^2/s
    DCST_RAFT = 0#0.1e-12#0.01e-12 #0.05 um^2/s
    DCST_MEM_2 = 0.05e-12 #0.05 um^2/s
    DCST_MEM_VES = 0.05e-12 #0.05 um^2/s
    #Ca model
    k1_pmca = 1.5e8
    k2_pmca = 15
    k3_pmca = 12
    kl_pmca = 4.3
    kon_serca_ca = 17147e6
    koff_serca_ca = 8426.3
    kflux_serca = 250
    kleak = 1.8e-06
    kon_pv_ca = 107e06
    koff_pv_ca = 0.95
    kon_pv_mg = 472e06
    koff_pv_mg = 25.0
    kon_cbs_ca = 5.5e06
    koff_cbs_ca = 2.6
    kon_cbf_ca = 43.5e06
    koff_cbf_ca = 35.8
    kinflux = 0.0
    
    ##### Species #####
    Ca = Species.Create(valence=2)
    PMCA_P0, PMCA_P1, SERCA, CaSERCA, Ca2SERCA, PV, PV_Ca, PV_2Ca, MgPV, Mg2PV, Mg, CBs, CaCBs, Ca2CBs, CBf, CaCBf, Ca2CBf = Species.Create()

    # CALBINDIN 1 (D-28K) (we consider the 2 (hi-aff) : 2 (lo-aff) scenario)
    # Binding Kinetics of Calbindin-D28k Determined by Flash Photolysis of Caged Ca2 (Nagerl 2000)
    CBhi, CBhi_Ca, CBhi_2Ca, CBlo, CBlo_Ca, CBlo_2Ca = Species.Create()

    CRTT, CRTR_Ca, CRRR_2Ca, CRind, CRind_Ca = Species.Create()
    syt, glu, SNARE_syt, SNARE_syt_CXN, SNARE_syt_CXN_Ca, SNARE_syt_CXN_Ca2, SNARE_syt_CXN_Ca3 = Species.Create()
    SNARE_syt_CXN_bCa, SNARE_syt_CXN_bCa2, SNARE_syt_CXN_Ca_bCa, SNARE_syt_CXN_Ca_bCa2, SNARE_syt_CXN_Ca2_bCa = Species.Create()
    SNARE_syt_CXN_Ca2_bCa2, SNARE_syt_CXN_Ca3_bCa, SNARE_syt_CXN_Ca3_bCa2 = Species.Create()
    # Vesicle docking
    RIM, M13, RIM_M13, Rab3, RIM_M13_Rab3, SYX, M18, SYX_M18, RIM_M13_Rab3_SYX_M18, RIM_M13_Rab3_SYX_M18_SNP25, SNP25, SYB, SNARE, CXN = Species.Create()
    aSNAP, SNARE_aSNAP, NSF, cisSNARE, SNARE_aSNAP_NSF, SNARE_aSNAP_NSF_1, SNARE_aSNAP_NSF_2, SNARE_aSNAP_NSF_3, SNARE_aSNAP_NSF_4 = Species.Create()
    SNARE_aSNAP_NSF_5, SNARE_aSNAP_NSF_6, SNARE_aSNAP_NSF_7, SNARE_aSNAP_NSF_8, SNARE_aSNAP_NSF_9, SNARE_DISS = Species.Create()

    N0, N1, N2, C0, C1, C2 = SubUnitState.Create()
    CaM_N_SU = SubUnit.Create([N0, N1, N2])
    CaM_C_SU = SubUnit.Create([C0, C1, C2])
    CaM = Complex.Create([CaM_N_SU, CaM_C_SU], statesAsSpecies=True)

    CaN, CaN_CaM, DYNpp, DYNp, DYN, CaN_DYNpp, CaN_DYNp = Species.Create()
    SYN1, DYN_SYN1, CDK5, DYN_CDK5, DYNp_CDK5, DYN_SYN1_CDK5, AP180, AP2, RabAd, DYN_AD, AC18, AC18_CaM, cAMP, PKA = Species.Create()
    R2C2, R2C2_cAMP, R2C2_2cAMPbb, R2C2_2cAMPab, R2C2_3cAMP, R2C2_4cAMP, R2C_2cAMPab, R2C_3cAMP, R2C_4cAMP, R2_4cAMP = Species.Create()
    synapsin, synapsin_p, PKA_synapsin, PP2A, PP2A_synapsin_p, actin, synapsin_actin, synapsin_actin_p, Rab3_TOMO_synapsin_actin = Species.Create()
    Rab3_TOMO_p_synapsin_actin, Rab3_TOMO_synapsin_actin_p, Rab3_TOMO_p_synapsin_actin_p, ves_syn1_site, Rab3_TOMO_synapsin = Species.Create()
    Rab3_TOMO_synapsin_p, Rab3_TOMO_p_synapsin, Rab3_TOMO_p_synapsin_p, TOMO_p, Rab3_TOMO_p, TOMO, Rab3_TOMO, TOMO_CDK5, Rab3_TOMO_CDK5 = Species.Create()
    TOMO_p_CaN_CaM, Rab3_TOMO_p_CaN_CaM, Rab3_TOMO_synapsin_PKA, Rab3_TOMO_p_synapsin_PKA, Rab3_TOMO_p_synapsin_p_PP2A, Rab3_TOMO_synapsin_p_PP2A = Species.Create()
    Rab3_TOMO_synapsin_CDK5, Rab3_TOMO_synapsin_p_CDK5, Rab3_TOMO_p_synapsin_CaN_CaM, Rab3_TOMO_p_synapsin_p_CaN_CaM = Species.Create()
    RIM_M13_Rab3_SYX_M18_SNP25_TOMO, RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p, RIM_M13_Rab3_SYX_M18_SNP25_TOMOx, RIM_M13_Rab3_SYX_M18_SNP25_TOMOx_p = Species.Create()
    RIM_M13_Rab3_SYX_M18_SNP25_TOMO_P, RIM_M13_Rab3_SYX_M18_SNP25_TOMOx_P, RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p_P, RIM_M13_Rab3_SYX_M18_SNP25_TOMOx_p_P = Species.Create()
    RIM_M13_Rab3_SYX_M18_SNP25_TOMO_PKA, RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p_PKA, RIM_M13_Rab3_SYX_M18_SNP25_TOMO_P_PP2A, RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p_P_PP2A = Species.Create()
    ves_CXN_site, ves_rab3_site, ves_m13_site, ves_m18_site, ves_syndap_site, ves_asnap_site, ves_nsf_site, REL_IND = Species.Create()

    #Basic states of dimerised synapsins (each "half" of the link between vesicles)
    #None of the dimers need explicitly defining.
    synapsin_dimer = LinkSpecies.Create(1e-13)
    synapsin_dimer_p = LinkSpecies.Create(1e-13)
    Rab3_TOMO_synapsin_dimer = LinkSpecies.Create(1e-13) #(2 formation routes)
    Rab3_TOMO_synapsin_dimer_p = LinkSpecies.Create(1e-13) #(2 formation routes)
    Rab3_TOMO_p_synapsin_dimer = LinkSpecies.Create(1e-13) #(2 formation routes)
    Rab3_TOMO_p_synapsin_dimer_p = LinkSpecies.Create(1e-13) #(2 formation routes)
    #Enzyme-complex states of ves-bound dimerised synapsins (assume the complexed enzyme blocks the dimerisation/unbinding reaction)
    PKA_synapsin_dimer = LinkSpecies.Create(1e-13)
    Rab3_TOMO_synapsin_dimer_PKA = LinkSpecies.Create(1e-13)
    Rab3_TOMO_p_synapsin_dimer_PKA = LinkSpecies.Create(1e-13)
    #PP2A dephos.
    PP2A_synapsin_dimer_p = LinkSpecies.Create(1e-13)
    Rab3_TOMO_synapsin_dimer_p_PP2A = LinkSpecies.Create(1e-13)
    Rab3_TOMO_p_synapsin_dimer_p_PP2A = LinkSpecies.Create(1e-13)
    #CDK5 phos.
    Rab3_TOMO_synapsin_dimer_CDK5 = LinkSpecies.Create(1e-13)
    Rab3_TOMO_synapsin_dimer_p_CDK5 = LinkSpecies.Create(1e-13)
    #CaN dephos.
    Rab3_TOMO_p_synapsin_dimer_CaN_CaM = LinkSpecies.Create(1e-13)
    Rab3_TOMO_p_synapsin_dimer_p_CaN_CaM = LinkSpecies.Create(1e-13)

    # Channels
    CaPc, CaPo = SubUnitState.Create()
    CaP_SU = SubUnit.Create([CaPc, CaPo])
    CaPchan = Channel.Create([CaP_SU, CaP_SU, CaP_SU], statesAsSpecies=True)

    #######VESICLES#########
    ves_diameter = 40e-9
    ves_radius = ves_diameter/2
    ves_diff_k = 0.11e-12 #Gaffield 2006 Neuron
    ves_diff_k_slow = 0.0011e-12 #Rizzoli
        
    ves_ssys = VesicleSurfaceSystem.Create()
    #free vesicles
    ves = Vesicle.Create(ves_diameter, ves_diff_k, ves_ssys)#0.003e-12

    raft_sys = RaftSurfaceSystem.Create()
    raft = Raft.Create(3e-08, DCST_RAFT, raft_sys)

    kon_dyn_syn = 1000e06
    koff_dyn_syn = 1
    kon_cdk_dyn = 3e06
    koff_cdk_dyn = 10
    kcat_cdk = 0.5

    kon_syb_ap180 = 2000e06
    kon_syt_ap2 = 2000e06
    kon_rab3_rabad = 2000e06
    kon_dyn_raft = 2000e06
    koff_dyn_raft = 0
    #Parameters
    kon_ac_cam = 500e06
    koff_ac_cam = 0.1
    kcat_ac18 = 18
    kdeg_camp = 1
    #Parameters
    kon_r2c2_camp1 = 2e06
    koff_r2c2_camp1 = 0.75
    
    kon_r2c2_camp2bb = 1e06
    koff_r2c2_camp2bb = 1.5
    
    kon_r2c2_camp2ab = 10e06
    koff_r2c2_camp2ab = 7.5
    
    kon_r2c2_camp3bb = 20e06
    koff_r2c2_camp3bb = 7.5
    
    kon_r2c2_camp3ab = 1e06
    koff_r2c2_camp3ab = 0.75
    
    kon_r2c2_camp4 = 10e06
    koff_r2c2_camp4 = 15
    
    kact_pka_23 = 0.005
    kdeact_pka_23 = 5e06
    
    kact_pka_4 = 6
    kdeact_pka_4 = 5e06
    
    kon_r2c_camp3 = 1e06
    koff_r2c_camp3 = 0.75
    
    kon_r2c_camp4 = 10e06
    koff_r2c_camp4 = 7.5
    
    kact_pka_r2c = 3
    kdeact_pka_r2c = 10e06

    kon_PKA_syn = 10e06
    koff_PKA_syn = 1
    kcat_PKA_syn = 10

    kon_pp2a_synapsin = 3e06 #Gallimore, Cell Reports, 2018 (from PP2A CaMKII dephos parameters)
    koff_pp2a_synapsin = 10
    kcat_pp2a_syn = 0.5

    #synapsin bound to vesicle binds to actin
    kon_synapsin_actin = 1000e06
    koff_synapsin_actin = 1
    koff_synapsin_actin_p = 1000

    kon_cdk_tomo = 3e06
    koff_cdk_tomo = 10
    kcat_cdk_tomo = 0.5

    kon_can_tomo = 1e06
    koff_can_tomo = 10
    #TOMOSYN binds Rab3
    kon_rab3_TOMO_p = 1e06
    koff_rab3_TOMO_p = 10
    #Rab3-TOMO binds synapsin (non-dimerised)
    kon_TOMO_p_synapsin = 100e06
    koff_TOMO_p_synapsin = 1
    kon_TOMO_synapsin = 100e06
    koff_TOMO_synapsin = 100 #increase 100x with TOMO_p dephosphorylation

    kon_synapsin_dimer = 10e06
    koff_synapsin_dimer = 10
    koff_synapsin_dimer_p = 10000

    # #2. Inhibition of vesicle priming (competition with SYB, same rate)
    kon_tomo_syx = 2.35e06
    koff_tomo_syx = 0.0047
    kdisp_tomo_syx = 10e06
    kflip_tomo_f = 0.01
    kflip_tomo_P_f = 10
    kflip_tomo_r = 1
    kon_tomo_pka = 10e06
    koff_tomo_pka = 1
    kcat_pka_tomo = 50

    kon_tomo_PP2A = 2e06 #Gallimore, Cell Reports, 2018 (from PP2A CaMKII dephos parameters)
    koff_tomo_PP2A = 8
    kcat_PP2A_tomo = 2

    if ves_protein_bind_on:
        kon_cxn_ves = 1.2e06
    else:
        kon_cxn_ves = 0
    koff_cxn_ves = 1000
    if ves_protein_bind_on:
        kon_rab3_ves = 10e06
    else:
        kon_rab3_ves = 0
    koff_rab3_ves = 1
    if ves_protein_bind_on:
        kon_m13_ves = 1.6e06
    else:
        kon_m13_ves = 0
    koff_m13_ves = 1000
    if ves_protein_bind_on:
        kon_m18_ves = 1.6e06
    else:
        kon_m18_ves = 0
    koff_m18_ves = 1000
    if ves_protein_bind_on:
        kon_syn1_ves = 2.8e06
    else:
        kon_syn1_ves = 0
    koff_syn1_ves = 1000
    if ves_protein_bind_on:
        kon_asnap_ves = 2.8e06
    else:
        kon_asnap_ves = 0
    koff_asnap_ves = 1000
    if ves_protein_bind_on:
        kon_nsf_ves = 0.9e06
    else:
        kon_nsf_ves = 0
    koff_nsf_ves = 1000


    with vsys:
        
        #Diffusion rules
        diff_Ca = Diffusion.Create(Ca, DCST_CA)
        for spec in [PV, PV_Ca, PV_2Ca, MgPV, Mg2PV, Mg, CBs, CaCBs, Ca2CBs, CBf, CaCBf, Ca2CBf]:
            Diffusion(spec, DCST_CYT)

        #Reactions
        #Ca uncage
        None >r['Ca_uncage']> Ca
        r['Ca_uncage'].K = 0.0

        (Ca + PV <r[1]> PV_Ca) + Ca <r[2]> PV_2Ca
        r[1].K = 107000000.0, 0.95
        r[2].K = 107000000.0, 0.95
        
        (CBhi + Ca <r[1]> CBhi_Ca) + Ca <r[2]> CBhi_2Ca
        r[1].K = 11000000.0, 2.607
        r[2].K = 11000000.0, 2.607
        
        (CBlo + Ca <r[1]> CBlo_Ca) + Ca <r[2]> CBlo_2Ca
        r[1].K = 87000000.0, 35.76
        r[2].K = 87000000.0, 35.76

        kon_T  = 1.8e6
        koff_T = 53
        kon_R  = 3.1e8
        koff_R = 20
        kon_ind  = 7.3e6
        koff_ind = 252
        
        # pair 1 (you MUST multiply its concentration by 2) because WE'VE GOT TWO PAIRS OF COOPERATIVE CA2+ BINDING SITES
        (CRTT + Ca <r[1]> CRTR_Ca) + Ca <r[2]> CRRR_2Ca
        r[1].K = 2*kon_T, koff_T
        r[2].K = kon_T, 2*koff_R
        
        # independent Ca2+ binding site:
        CRind + Ca <r[1]> CRind_Ca
        r[1].K = kon_ind, koff_ind

        diff_M13 = Diffusion.Create(M13, 2.2796e-12)
        diff_Rab3 = Diffusion.Create(Rab3, 1.7793e-12)
        diff_M18 = Diffusion.Create(M18, 2.2357e-12)
        diff_CXN = Diffusion.Create(CXN, 1.7793e-12)

        diff_aSNAP = Diffusion.Create(aSNAP, 2.1959e-12)
        diff_NSF = Diffusion.Create(NSF, 2.2788e-12)


        kon_cam_NT = 7.7e08
        kon_cam_NR = 3.2e10
        kon_cam_CT = 8.4e07
        kon_cam_CR = 2.5e07
        
        koff_cam_NT = 1.6e05
        koff_cam_CT = 2.2e04
        koff_cam_NR = 2.6e03
        koff_cam_CR = 6.5
        
        Diffusion(CaM[...], 1.6496e-12)
        with CaM[...]:
            N0 + Ca <r[1]> N1
            N1 + Ca <r[2]> N2
            r[1].K = 2 * kon_cam_NT, koff_cam_NT
            r[2].K = kon_cam_NR, 2 * koff_cam_NR

            C0 + Ca <r[3]> C1
            C1 + Ca <r[4]> C2
            r[3].K = 2 * kon_cam_CT, koff_cam_CT
            r[4].K = kon_cam_CR, 2 * koff_cam_CR
        # NOTE: The following two reactions are added to be identical to the original API_1 script, but they might not be required
        Ca + CaM[N1, C1] >r[1]> CaM[N2, C1] # Already declared earlier
        r[1].K = kon_cam_NR
        Ca + CaM[N1, C1] >r[1]> CaM[N1, C2] # Already declared earlier
        r[1].K = kon_cam_CR

        diff_CaN = Diffusion.Create(CaN, DCST_CYT)
        diff_CaN_CaM = Diffusion.Create(CaN_CaM, DCST_CYT)

        
        kon_can_cam = 46e06
        koff_can_cam = 1.2
        
        kon_can_dyn = 10e06
        koff_can_dyn = 1
        kcat_can = 0.5
        
        #Calcineurin binding to CaM
        CaM[N2, C2] + CaN <r[1]> CaN_CaM
        r[1].K = kon_can_cam, koff_can_cam

        diff_DYNpp = Diffusion.Create(DYNpp, 3.2464e-12)
        diff_DYNp = Diffusion.Create(DYNp, 3.2464e-12)
        diff_DYN = Diffusion.Create(DYN, 3.2464e-12)

        #cytosolic dephosphorylation
        CaN_CaM + DYNpp <r[1]> CaN_DYNpp
        r[1].K = kon_can_dyn, koff_can_dyn
        CaN_DYNpp >r[1]> CaN_CaM + DYNp
        r[1].K = kcat_can
        
        CaN_CaM + DYNp <r[1]> CaN_DYNp
        r[1].K = kon_can_dyn, koff_can_dyn
        CaN_DYNp >r[1]> CaN_CaM + DYN
        r[1].K = kcat_can

        diff_SYN1 = Diffusion.Create(SYN1, 2.1959e-12)
        diff_DYN_SYN1 = Diffusion.Create(DYN_SYN1, 3.2464e-12)

        DYN + SYN1 <r[1]> DYN_SYN1
        r[1].K = kon_dyn_syn, koff_dyn_syn
        
        diff_CDK5 = Diffusion.Create(CDK5, DCST_CYT)
        diff_DYN_CDK5 = Diffusion.Create(DYN_CDK5, DCST_CYT)
        diff_DYNp_CDK5 = Diffusion.Create(DYNp_CDK5, DCST_CYT)
        diff_DYN_SYN1_CDK5 = Diffusion.Create(DYN_SYN1_CDK5, DCST_CYT)
        
        #phosphorylation of dynamin in the cytosol
        CDK5 + DYN <r[1]> DYN_CDK5
        r[1].K = kon_cdk_dyn, koff_cdk_dyn
        DYN_CDK5 >r[1]> CDK5 + DYNp
        r[1].K = kcat_cdk

        CDK5 + DYN_SYN1 <r[1]> DYN_SYN1_CDK5
        r[1].K = kon_cdk_dyn, koff_cdk_dyn
        DYN_SYN1_CDK5 >r[1]> CDK5 + DYNp + SYN1
        r[1].K = kcat_cdk

        diff_AC18 = Diffusion.Create(AC18, DCST_CYT)
        diff_AC18_CaM = Diffusion.Create(AC18_CaM, DCST_CYT)
        diff_cAMP = Diffusion.Create(cAMP, DCST_CYT)

        AC18 + CaM[N2, C2] <r[1]> AC18_CaM
        r[1].K = kon_ac_cam, koff_ac_cam
        AC18_CaM >r[1]> AC18_CaM + cAMP
        r[1].K = kcat_ac18
        cAMP >r[1]> None
        r[1].K = kdeg_camp

        diff_PKA = Diffusion.Create(PKA, DCST_CYT)
        diff_R2C2 = Diffusion.Create(R2C2, DCST_CYT)
        diff_R2C2_cAMP = Diffusion.Create(R2C2_cAMP, DCST_CYT)
        diff_R2C2_2cAMPbb = Diffusion.Create(R2C2_2cAMPbb, DCST_CYT)
        diff_R2C2_2cAMPab = Diffusion.Create(R2C2_2cAMPab, DCST_CYT)
        diff_R2C2_3cAMP = Diffusion.Create(R2C2_3cAMP, DCST_CYT)
        diff_R2C2_4cAMP = Diffusion.Create(R2C2_4cAMP, DCST_CYT)
        diff_R2C_2cAMPab = Diffusion.Create(R2C_2cAMPab, DCST_CYT)
        diff_R2C_3cAMP = Diffusion.Create(R2C_3cAMP, DCST_CYT)
        diff_R2C_4cAMP = Diffusion.Create(R2C_4cAMP, DCST_CYT)
        diff_R2_4cAMP = Diffusion.Create(R2_4cAMP, DCST_CYT)

        #Reactions
        R2C2 + cAMP <r[1]> R2C2_cAMP
        r[1].K = kon_r2c2_camp1, koff_r2c2_camp1
        R2C2_cAMP + cAMP <r[1]> R2C2_2cAMPbb
        r[1].K = kon_r2c2_camp2bb, koff_r2c2_camp2bb
        R2C2_cAMP + cAMP <r[1]> R2C2_2cAMPab
        r[1].K = kon_r2c2_camp2ab, koff_r2c2_camp2ab
        R2C2_2cAMPbb + cAMP <r[1]> R2C2_3cAMP
        r[1].K = kon_r2c2_camp3bb, koff_r2c2_camp3bb
        R2C2_2cAMPab + cAMP >r[1]> R2C2_3cAMP
        r[1].K = kon_r2c2_camp3ab
        R2C2_3cAMP >r[1]> R2C2_2cAMPab
        r[1].K = koff_r2c2_camp3ab
        R2C2_3cAMP + cAMP <r[1]> R2C2_4cAMP
        r[1].K = kon_r2c2_camp4, koff_r2c2_camp4
        R2C2_2cAMPab <r[1]> PKA + R2C_2cAMPab
        r[1].K = kact_pka_23, kdeact_pka_23
        R2C2_3cAMP <r[1]> PKA + R2C_3cAMP
        r[1].K = kact_pka_23, kdeact_pka_23
        R2C2_4cAMP <r[1]> PKA + R2C_4cAMP
        r[1].K = kact_pka_4, kdeact_pka_4
        R2C_2cAMPab + cAMP <r[1]> R2C_3cAMP
        r[1].K = kon_r2c_camp3, koff_r2c_camp3
        R2C_3cAMP + cAMP <r[1]> R2C_4cAMP
        r[1].K = kon_r2c_camp4, koff_r2c_camp4
        R2C_4cAMP <r[1]> PKA + R2_4cAMP
        r[1].K = kact_pka_r2c, kdeact_pka_r2c

        diff_synapsin = Diffusion.Create(synapsin, 2.65e-12)
        diff_synapsin_p = Diffusion.Create(synapsin_p, 2.65e-12)
        diff_synapsin_pka = Diffusion.Create(PKA_synapsin, 2.65e-12)
        diff_synapsin_pp2a = Diffusion.Create(PP2A_synapsin_p, 2.65e-12)

        #dephos of cytosolic synapsin_p
        PP2A + synapsin_p <r[1]> PP2A_synapsin_p
        r[1].K = kon_pp2a_synapsin, koff_pp2a_synapsin
        PP2A_synapsin_p >r[1]> PP2A + synapsin
        r[1].K = kcat_pp2a_syn

        diff_actin = Diffusion.Create(actin, 0.0)
        diff_synapsin_actin = Diffusion.Create(synapsin_actin, 0.0)
        diff_synapsin_p_actin = Diffusion.Create(synapsin_actin_p, 0.0)
        diff_Rab3_TOMO_synapsin_actin = Diffusion.Create(Rab3_TOMO_synapsin_actin, 0.0)
        diff_Rab3_TOMO_p_synapsin_actin = Diffusion.Create(Rab3_TOMO_p_synapsin_actin, 0.0)
        diff_TOMO_p = Diffusion.Create(TOMO_p, 1.7793e-12)
        diff_TOMO = Diffusion.Create(TOMO, 1.7793e-12)
        diff_TOMO_CDK5 = Diffusion.Create(TOMO_CDK5, 1.7793e-12)

        CDK5 + TOMO <r[1]> TOMO_CDK5
        r[1].K = kon_cdk_tomo, koff_cdk_tomo
        TOMO_CDK5 >r[1]> CDK5 + TOMO_p
        r[1].K = kcat_cdk_tomo

        diff_TOMO_p_CaN = Diffusion.Create(TOMO_p_CaN_CaM, 1.7793e-12)
        
        CaN_CaM + TOMO_p <r[1]> TOMO_p_CaN_CaM
        r[1].K = kon_can_tomo, koff_can_tomo
        TOMO_p_CaN_CaM >r[1]> CaN_CaM + TOMO
        r[1].K = kcat_can_tomo

        #initial dimerisation depends on synapsin_actin dimer [5 different forms] (ensures actin is the focus of the clustering)
        synapsin_dimerise_1a = VesicleBind.Create((ves, ves), (synapsin, synapsin), (synapsin_dimer, synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_actin.v, None))
        synapsin_dimerise_1b = VesicleBind.Create((ves, ves), (synapsin, synapsin), (synapsin_dimer, synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_actin.v, None))
        synapsin_dimerise_1c = VesicleBind.Create((ves, ves), (synapsin, synapsin), (synapsin_dimer, synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_actin.v, None))
        synapsin_dimerise_1d = VesicleBind.Create((ves, ves), (synapsin, synapsin), (synapsin_dimer, synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_actin_p.v, None))
        synapsin_dimerise_1e = VesicleBind.Create((ves, ves), (synapsin, synapsin), (synapsin_dimer, synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_actin_p.v, None))
        
        #Dimer states (need both VesBind and VesUnbind reactions for each, and ldeps for VesBinds -- 6, one for each state)
        #synapsin - synapsin -> synapsin_dimer - synapsin_dimer
        synapsin_synapsin_dimerise_a = VesicleBind.Create((ves, ves), (synapsin, synapsin), (synapsin_dimer, synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer.v, None))
        synapsin_synapsin_dimerise_b = VesicleBind.Create((ves, ves), (synapsin, synapsin), (synapsin_dimer, synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer_p.v, None))
        synapsin_synapsin_dimerise_c = VesicleBind.Create((ves, ves), (synapsin, synapsin), (synapsin_dimer, synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer.v, None))
        synapsin_synapsin_dimerise_d = VesicleBind.Create((ves, ves), (synapsin, synapsin), (synapsin_dimer, synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer_p.v, None))
        synapsin_synapsin_dimerise_e = VesicleBind.Create((ves, ves), (synapsin, synapsin), (synapsin_dimer, synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer.v, None))
        synapsin_synapsin_dimerise_f = VesicleBind.Create((ves, ves), (synapsin, synapsin), (synapsin_dimer, synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer_p.v, None))
        
        synapsin_synapsin_unbind = VesicleUnbind.Create((ves, ves), (synapsin_dimer, synapsin_dimer), (synapsin, synapsin), kcst=koff_synapsin_dimer, immobilization=NO_EFFECT)
        
        #synapsin - synapsin_p -> synapsin_dimer - synapsin_dimer_p
        synapsin_synapsin_p_dimerise_a = VesicleBind.Create((ves, ves), (synapsin, synapsin_p), (synapsin_dimer, synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer.v, None))
        synapsin_synapsin_p_dimerise_b = VesicleBind.Create((ves, ves), (synapsin, synapsin_p), (synapsin_dimer, synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer_p.v, None))
        synapsin_synapsin_p_dimerise_c = VesicleBind.Create((ves, ves), (synapsin, synapsin_p), (synapsin_dimer, synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer.v, None))
        synapsin_synapsin_p_dimerise_d = VesicleBind.Create((ves, ves), (synapsin, synapsin_p), (synapsin_dimer, synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer_p.v, None))
        synapsin_synapsin_p_dimerise_e = VesicleBind.Create((ves, ves), (synapsin, synapsin_p), (synapsin_dimer, synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer.v, None))
        synapsin_synapsin_p_dimerise_f = VesicleBind.Create((ves, ves), (synapsin, synapsin_p), (synapsin_dimer, synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer_p.v, None))
        
        synapsin_synapsin_p_unbind = VesicleUnbind.Create((ves, ves), (synapsin_dimer, synapsin_dimer_p), (synapsin, synapsin_p), kcst=koff_synapsin_dimer_p, immobilization=NO_EFFECT)
        
        #synapsin - Rab3_TOMO_synapsin -> synapsin_dimer - Rab3_TOMO_synapsin_dimer
        synapsin_Rab3_TOMO_synapsin_dimerise_a = VesicleBind.Create((ves, ves), (synapsin, Rab3_TOMO_synapsin), (synapsin_dimer, Rab3_TOMO_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer.v, None))
        synapsin_Rab3_TOMO_synapsin_dimerise_b = VesicleBind.Create((ves, ves), (synapsin, Rab3_TOMO_synapsin), (synapsin_dimer, Rab3_TOMO_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer_p.v, None))
        synapsin_Rab3_TOMO_synapsin_dimerise_c = VesicleBind.Create((ves, ves), (synapsin, Rab3_TOMO_synapsin), (synapsin_dimer, Rab3_TOMO_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer.v, None))
        synapsin_Rab3_TOMO_synapsin_dimerise_d = VesicleBind.Create((ves, ves), (synapsin, Rab3_TOMO_synapsin), (synapsin_dimer, Rab3_TOMO_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer_p.v, None))
        synapsin_Rab3_TOMO_synapsin_dimerise_e = VesicleBind.Create((ves, ves), (synapsin, Rab3_TOMO_synapsin), (synapsin_dimer, Rab3_TOMO_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer.v, None))
        synapsin_Rab3_TOMO_synapsin_dimerise_f = VesicleBind.Create((ves, ves), (synapsin, Rab3_TOMO_synapsin), (synapsin_dimer, Rab3_TOMO_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer_p.v, None))
        
        synapsin_Rab3_TOMO_synapsin_unbind = VesicleUnbind.Create((ves, ves), (synapsin_dimer, Rab3_TOMO_synapsin_dimer), (synapsin, Rab3_TOMO_synapsin), kcst=koff_synapsin_dimer, immobilization=NO_EFFECT)
        
        #synapsin - Rab3_TOMO_synapsin_p -> synapsin_dimer - Rab3_TOMO_synapsin_dimer_p
        synapsin_Rab3_TOMO_synapsin_p_dimerise_a = VesicleBind.Create((ves, ves), (synapsin, Rab3_TOMO_synapsin_p), (synapsin_dimer, Rab3_TOMO_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer.v, None))
        synapsin_Rab3_TOMO_synapsin_p_dimerise_b = VesicleBind.Create((ves, ves), (synapsin, Rab3_TOMO_synapsin_p), (synapsin_dimer, Rab3_TOMO_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer_p.v, None))
        synapsin_Rab3_TOMO_synapsin_p_dimerise_c = VesicleBind.Create((ves, ves), (synapsin, Rab3_TOMO_synapsin_p), (synapsin_dimer, Rab3_TOMO_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer.v, None))
        synapsin_Rab3_TOMO_synapsin_p_dimerise_d = VesicleBind.Create((ves, ves), (synapsin, Rab3_TOMO_synapsin_p), (synapsin_dimer, Rab3_TOMO_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer_p.v, None))
        synapsin_Rab3_TOMO_synapsin_p_dimerise_e = VesicleBind.Create((ves, ves), (synapsin, Rab3_TOMO_synapsin_p), (synapsin_dimer, Rab3_TOMO_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer.v, None))
        synapsin_Rab3_TOMO_synapsin_p_dimerise_f = VesicleBind.Create((ves, ves), (synapsin, Rab3_TOMO_synapsin_p), (synapsin_dimer, Rab3_TOMO_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer_p.v, None))
        
        synapsin_Rab3_TOMO_synapsin_p_unbind = VesicleUnbind.Create((ves, ves), (synapsin_dimer, Rab3_TOMO_synapsin_dimer_p), (synapsin, Rab3_TOMO_synapsin_p), kcst=koff_synapsin_dimer_p, immobilization=NO_EFFECT)
        
        #synapsin - Rab3_TOMO_p_synapsin -> synapsin_dimer - Rab3_TOMO_p_synapsin_dimer
        synapsin_Rab3_TOMO_p_synapsin_dimerise_a = VesicleBind.Create((ves, ves), (synapsin, Rab3_TOMO_p_synapsin), (synapsin_dimer, Rab3_TOMO_p_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer.v, None))
        synapsin_Rab3_TOMO_p_synapsin_dimerise_b = VesicleBind.Create((ves, ves), (synapsin, Rab3_TOMO_p_synapsin), (synapsin_dimer, Rab3_TOMO_p_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer_p.v, None))
        synapsin_Rab3_TOMO_p_synapsin_dimerise_c = VesicleBind.Create((ves, ves), (synapsin, Rab3_TOMO_p_synapsin), (synapsin_dimer, Rab3_TOMO_p_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer.v, None))
        synapsin_Rab3_TOMO_p_synapsin_dimerise_d = VesicleBind.Create((ves, ves), (synapsin, Rab3_TOMO_p_synapsin), (synapsin_dimer, Rab3_TOMO_p_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer_p.v, None))
        synapsin_Rab3_TOMO_p_synapsin_dimerise_e = VesicleBind.Create((ves, ves), (synapsin, Rab3_TOMO_p_synapsin), (synapsin_dimer, Rab3_TOMO_p_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer.v, None))
        synapsin_Rab3_TOMO_p_synapsin_dimerise_f = VesicleBind.Create((ves, ves), (synapsin, Rab3_TOMO_p_synapsin), (synapsin_dimer, Rab3_TOMO_p_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer_p.v, None))
        
        synapsin_Rab3_TOMO_p_synapsin_unbind = VesicleUnbind.Create((ves, ves), (synapsin_dimer, Rab3_TOMO_p_synapsin_dimer), (synapsin, Rab3_TOMO_p_synapsin), kcst=koff_synapsin_dimer, immobilization=NO_EFFECT)
        
        #synapsin - Rab3_TOMO_p_synapsin_p -> synapsin_dimer - Rab3_TOMO_p_synapsin_dimer_p
        synapsin_Rab3_TOMO_p_synapsin_p_dimerise_a = VesicleBind.Create((ves, ves), (synapsin, Rab3_TOMO_p_synapsin_p), (synapsin_dimer, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer.v, None))
        synapsin_Rab3_TOMO_p_synapsin_p_dimerise_b = VesicleBind.Create((ves, ves), (synapsin, Rab3_TOMO_p_synapsin_p), (synapsin_dimer, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer_p.v, None))
        synapsin_Rab3_TOMO_p_synapsin_p_dimerise_c = VesicleBind.Create((ves, ves), (synapsin, Rab3_TOMO_p_synapsin_p), (synapsin_dimer, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer.v, None))
        synapsin_Rab3_TOMO_p_synapsin_p_dimerise_d = VesicleBind.Create((ves, ves), (synapsin, Rab3_TOMO_p_synapsin_p), (synapsin_dimer, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer_p.v, None))
        synapsin_Rab3_TOMO_p_synapsin_p_dimerise_e = VesicleBind.Create((ves, ves), (synapsin, Rab3_TOMO_p_synapsin_p), (synapsin_dimer, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer.v, None))
        synapsin_Rab3_TOMO_p_synapsin_p_dimerise_f = VesicleBind.Create((ves, ves), (synapsin, Rab3_TOMO_p_synapsin_p), (synapsin_dimer, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer_p.v, None))
        
        synapsin_Rab3_TOMO_p_synapsin_p_unbind = VesicleUnbind.Create((ves, ves), (synapsin_dimer, Rab3_TOMO_p_synapsin_dimer_p), (synapsin, Rab3_TOMO_p_synapsin_p), kcst=koff_synapsin_dimer_p, immobilization=NO_EFFECT)
        
        #synapsin_p - synapsin_p -> synapsin_dimer_p - synapsin_dimer_p
        synapsin_p_synapsin_p_dimerise_a = VesicleBind.Create((ves, ves), (synapsin_p, synapsin_p), (synapsin_dimer_p, synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer.v, None))
        synapsin_p_synapsin_p_dimerise_b = VesicleBind.Create((ves, ves), (synapsin_p, synapsin_p), (synapsin_dimer_p, synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer_p.v, None))
        synapsin_p_synapsin_p_dimerise_c = VesicleBind.Create((ves, ves), (synapsin_p, synapsin_p), (synapsin_dimer_p, synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer.v, None))
        synapsin_p_synapsin_p_dimerise_d = VesicleBind.Create((ves, ves), (synapsin_p, synapsin_p), (synapsin_dimer_p, synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer_p.v, None))
        synapsin_p_synapsin_p_dimerise_e = VesicleBind.Create((ves, ves), (synapsin_p, synapsin_p), (synapsin_dimer_p, synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer.v, None))
        synapsin_p_synapsin_p_dimerise_f = VesicleBind.Create((ves, ves), (synapsin_p, synapsin_p), (synapsin_dimer_p, synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer_p.v, None))
        
        synapsin_p_synapsin_p_dimerise_unbind = VesicleUnbind.Create((ves, ves), (synapsin_dimer_p, synapsin_dimer_p), (synapsin_p, synapsin_p), kcst=koff_synapsin_dimer_p, immobilization=NO_EFFECT)
        
        #synapsin_p - Rab3_TOMO_synapsin -> synapsin_dimer_p - Rab3_TOMO_synapsin_dimer
        synapsin_p_Rab3_TOMO_synapsin_dimerise_a = VesicleBind.Create((ves, ves), (synapsin_p, Rab3_TOMO_synapsin), (synapsin_dimer_p, Rab3_TOMO_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer.v, None))
        synapsin_p_Rab3_TOMO_synapsin_dimerise_b = VesicleBind.Create((ves, ves), (synapsin_p, Rab3_TOMO_synapsin), (synapsin_dimer_p, Rab3_TOMO_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer_p.v, None))
        synapsin_p_Rab3_TOMO_synapsin_dimerise_c = VesicleBind.Create((ves, ves), (synapsin_p, Rab3_TOMO_synapsin), (synapsin_dimer_p, Rab3_TOMO_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer.v, None))
        synapsin_p_Rab3_TOMO_synapsin_dimerise_d = VesicleBind.Create((ves, ves), (synapsin_p, Rab3_TOMO_synapsin), (synapsin_dimer_p, Rab3_TOMO_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer_p.v, None))
        synapsin_p_Rab3_TOMO_synapsin_dimerise_e = VesicleBind.Create((ves, ves), (synapsin_p, Rab3_TOMO_synapsin), (synapsin_dimer_p, Rab3_TOMO_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer.v, None))
        synapsin_p_Rab3_TOMO_synapsin_dimerise_f = VesicleBind.Create((ves, ves), (synapsin_p, Rab3_TOMO_synapsin), (synapsin_dimer_p, Rab3_TOMO_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer_p.v, None))
        
        synapsin_p_Rab3_TOMO_synapsin_dimerise_unbind = VesicleUnbind.Create((ves, ves), (synapsin_dimer_p, Rab3_TOMO_synapsin_dimer), (synapsin_p, Rab3_TOMO_synapsin), kcst=koff_synapsin_dimer_p, immobilization=NO_EFFECT)
        
        #synapsin_p - Rab3_TOMO_synapsin_p -> synapsin_dimer_p - Rab3_TOMO_synapsin_dimer_p
        synapsin_p_Rab3_TOMO_synapsin_p_dimerise_a = VesicleBind.Create((ves, ves), (synapsin_p, Rab3_TOMO_synapsin_p), (synapsin_dimer_p, Rab3_TOMO_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer.v, None))
        synapsin_p_Rab3_TOMO_synapsin_p_dimerise_b = VesicleBind.Create((ves, ves), (synapsin_p, Rab3_TOMO_synapsin_p), (synapsin_dimer_p, Rab3_TOMO_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer_p.v, None))
        synapsin_p_Rab3_TOMO_synapsin_p_dimerise_c = VesicleBind.Create((ves, ves), (synapsin_p, Rab3_TOMO_synapsin_p), (synapsin_dimer_p, Rab3_TOMO_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer.v, None))
        synapsin_p_Rab3_TOMO_synapsin_p_dimerise_d = VesicleBind.Create((ves, ves), (synapsin_p, Rab3_TOMO_synapsin_p), (synapsin_dimer_p, Rab3_TOMO_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer_p.v, None))
        synapsin_p_Rab3_TOMO_synapsin_p_dimerise_e = VesicleBind.Create((ves, ves), (synapsin_p, Rab3_TOMO_synapsin_p), (synapsin_dimer_p, Rab3_TOMO_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer.v, None))
        synapsin_p_Rab3_TOMO_synapsin_p_dimerise_f = VesicleBind.Create((ves, ves), (synapsin_p, Rab3_TOMO_synapsin_p), (synapsin_dimer_p, Rab3_TOMO_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer_p.v, None))
        
        synapsin_p_Rab3_TOMO_synapsin_p_unbind = VesicleUnbind.Create((ves, ves), (synapsin_dimer_p, Rab3_TOMO_synapsin_dimer_p), (synapsin_p, Rab3_TOMO_synapsin_p), kcst=koff_synapsin_dimer_p, immobilization=NO_EFFECT)
        
        #synapsin_p - Rab3_TOMO_p_synapsin -> synapsin_dimer_p - Rab3_TOMO_p_synapsin_dimer
        synapsin_p_Rab3_TOMO_p_synapsin_dimerise_a = VesicleBind.Create((ves, ves), (synapsin_p, Rab3_TOMO_p_synapsin), (synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer.v, None))
        synapsin_p_Rab3_TOMO_p_synapsin_dimerise_b = VesicleBind.Create((ves, ves), (synapsin_p, Rab3_TOMO_p_synapsin), (synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer_p.v, None))
        synapsin_p_Rab3_TOMO_p_synapsin_dimerise_c = VesicleBind.Create((ves, ves), (synapsin_p, Rab3_TOMO_p_synapsin), (synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer.v, None))
        synapsin_p_Rab3_TOMO_p_synapsin_dimerise_d = VesicleBind.Create((ves, ves), (synapsin_p, Rab3_TOMO_p_synapsin), (synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer_p.v, None))
        synapsin_p_Rab3_TOMO_p_synapsin_dimerise_e = VesicleBind.Create((ves, ves), (synapsin_p, Rab3_TOMO_p_synapsin), (synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer.v, None))
        synapsin_p_Rab3_TOMO_p_synapsin_dimerise_f = VesicleBind.Create((ves, ves), (synapsin_p, Rab3_TOMO_p_synapsin), (synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer_p.v, None))
        
        synapsin_p_Rab3_TOMO_p_synapsin_unbind = VesicleUnbind.Create((ves, ves), (synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer), (synapsin_p, Rab3_TOMO_p_synapsin), kcst=koff_synapsin_dimer_p, immobilization=NO_EFFECT)
        
        #synapsin_p - Rab3_TOMO_p_synapsin_p -> synapsin_dimer_p - Rab3_TOMO_p_synapsin_dimer_p
        synapsin_p_Rab3_TOMO_p_synapsin_p_dimerise_a = VesicleBind.Create((ves, ves), (synapsin_p, Rab3_TOMO_p_synapsin_p), (synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer.v, None))
        synapsin_p_Rab3_TOMO_p_synapsin_p_dimerise_b = VesicleBind.Create((ves, ves), (synapsin_p, Rab3_TOMO_p_synapsin_p), (synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer_p.v, None))
        synapsin_p_Rab3_TOMO_p_synapsin_p_dimerise_c = VesicleBind.Create((ves, ves), (synapsin_p, Rab3_TOMO_p_synapsin_p), (synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer.v, None))
        synapsin_p_Rab3_TOMO_p_synapsin_p_dimerise_d = VesicleBind.Create((ves, ves), (synapsin_p, Rab3_TOMO_p_synapsin_p), (synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer_p.v, None))
        synapsin_p_Rab3_TOMO_p_synapsin_p_dimerise_e = VesicleBind.Create((ves, ves), (synapsin_p, Rab3_TOMO_p_synapsin_p), (synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer.v, None))
        synapsin_p_Rab3_TOMO_p_synapsin_p_dimerise_f = VesicleBind.Create((ves, ves), (synapsin_p, Rab3_TOMO_p_synapsin_p), (synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer_p.v, None))
        
        synapsin_p_Rab3_TOMO_p_synapsin_p_unbind = VesicleUnbind.Create((ves, ves), (synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer_p), (synapsin_p, Rab3_TOMO_p_synapsin_p), kcst=koff_synapsin_dimer_p, immobilization=NO_EFFECT)
        
        #Rab3_TOMO_synapsin - Rab3_TOMO_synapsin -> Rab3_TOMO_synapsin_dimer - Rab3_TOMO_synapsin_dimer
        Rab3_TOMO_synapsin_Rab3_TOMO_synapsin_dimerise_a = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin, Rab3_TOMO_synapsin), (Rab3_TOMO_synapsin_dimer, Rab3_TOMO_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer.v, None))
        Rab3_TOMO_synapsin_Rab3_TOMO_synapsin_dimerise_b = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin, Rab3_TOMO_synapsin), (Rab3_TOMO_synapsin_dimer, Rab3_TOMO_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer_p.v, None))
        Rab3_TOMO_synapsin_Rab3_TOMO_synapsin_dimerise_c = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin, Rab3_TOMO_synapsin), (Rab3_TOMO_synapsin_dimer, Rab3_TOMO_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer.v, None))
        Rab3_TOMO_synapsin_Rab3_TOMO_synapsin_dimerise_d = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin, Rab3_TOMO_synapsin), (Rab3_TOMO_synapsin_dimer, Rab3_TOMO_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer_p.v, None))
        Rab3_TOMO_synapsin_Rab3_TOMO_synapsin_dimerise_e = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin, Rab3_TOMO_synapsin), (Rab3_TOMO_synapsin_dimer, Rab3_TOMO_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer.v, None))
        Rab3_TOMO_synapsin_Rab3_TOMO_synapsin_dimerise_f = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin, Rab3_TOMO_synapsin), (Rab3_TOMO_synapsin_dimer, Rab3_TOMO_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer_p.v, None))
        
        Rab3_TOMO_synapsin_Rab3_TOMO_synapsin_unbind = VesicleUnbind.Create((ves, ves), (Rab3_TOMO_synapsin_dimer, Rab3_TOMO_synapsin_dimer), (Rab3_TOMO_synapsin, Rab3_TOMO_synapsin), kcst=koff_synapsin_dimer, immobilization=NO_EFFECT)
        
        #Rab3_TOMO_synapsin - Rab3_TOMO_synapsin_p -> Rab3_TOMO_synapsin_dimer - Rab3_TOMO_synapsin_dimer_p
        Rab3_TOMO_synapsin_Rab3_TOMO_synapsin_p_dimerise_a = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin, Rab3_TOMO_synapsin_p), (Rab3_TOMO_synapsin_dimer, Rab3_TOMO_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer.v, None))
        Rab3_TOMO_synapsin_Rab3_TOMO_synapsin_p_dimerise_b = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin, Rab3_TOMO_synapsin_p), (Rab3_TOMO_synapsin_dimer, Rab3_TOMO_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer_p.v, None))
        Rab3_TOMO_synapsin_Rab3_TOMO_synapsin_p_dimerise_c = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin, Rab3_TOMO_synapsin_p), (Rab3_TOMO_synapsin_dimer, Rab3_TOMO_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer.v, None))
        Rab3_TOMO_synapsin_Rab3_TOMO_synapsin_p_dimerise_d = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin, Rab3_TOMO_synapsin_p), (Rab3_TOMO_synapsin_dimer, Rab3_TOMO_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer_p.v, None))
        Rab3_TOMO_synapsin_Rab3_TOMO_synapsin_p_dimerise_e = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin, Rab3_TOMO_synapsin_p), (Rab3_TOMO_synapsin_dimer, Rab3_TOMO_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer.v, None))
        Rab3_TOMO_synapsin_Rab3_TOMO_synapsin_p_dimerise_f = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin, Rab3_TOMO_synapsin_p), (Rab3_TOMO_synapsin_dimer, Rab3_TOMO_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer_p.v, None))
        
        Rab3_TOMO_synapsin_Rab3_TOMO_synapsin_p_unbind = VesicleUnbind.Create((ves, ves), (Rab3_TOMO_synapsin_dimer, Rab3_TOMO_synapsin_dimer_p), (Rab3_TOMO_synapsin, Rab3_TOMO_synapsin_p), kcst=koff_synapsin_dimer_p, immobilization=NO_EFFECT)
        
        #Rab3_TOMO_synapsin - Rab3_TOMO_p_synapsin -> Rab3_TOMO_synapsin_dimer - Rab3_TOMO_p_synapsin_dimer
        Rab3_TOMO_synapsin_Rab3_TOMO_p_synapsin_dimerise_a = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin, Rab3_TOMO_p_synapsin), (Rab3_TOMO_synapsin_dimer, Rab3_TOMO_p_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer.v, None))
        Rab3_TOMO_synapsin_Rab3_TOMO_p_synapsin_dimerise_b = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin, Rab3_TOMO_p_synapsin), (Rab3_TOMO_synapsin_dimer, Rab3_TOMO_p_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer_p.v, None))
        Rab3_TOMO_synapsin_Rab3_TOMO_p_synapsin_dimerise_c = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin, Rab3_TOMO_p_synapsin), (Rab3_TOMO_synapsin_dimer, Rab3_TOMO_p_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer.v, None))
        Rab3_TOMO_synapsin_Rab3_TOMO_p_synapsin_dimerise_d = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin, Rab3_TOMO_p_synapsin), (Rab3_TOMO_synapsin_dimer, Rab3_TOMO_p_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer_p.v, None))
        Rab3_TOMO_synapsin_Rab3_TOMO_p_synapsin_dimerise_e = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin, Rab3_TOMO_p_synapsin), (Rab3_TOMO_synapsin_dimer, Rab3_TOMO_p_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer.v, None))
        Rab3_TOMO_synapsin_Rab3_TOMO_p_synapsin_dimerise_f = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin, Rab3_TOMO_p_synapsin), (Rab3_TOMO_synapsin_dimer, Rab3_TOMO_p_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer_p.v, None))
        
        Rab3_TOMO_synapsin_Rab3_TOMO_p_unbind = VesicleUnbind.Create((ves, ves), (Rab3_TOMO_synapsin_dimer, Rab3_TOMO_p_synapsin_dimer), (Rab3_TOMO_synapsin, Rab3_TOMO_p_synapsin), kcst=koff_synapsin_dimer, immobilization=NO_EFFECT)
        
        #Rab3_TOMO_synapsin - Rab3_TOMO_p_synapsin_p -> Rab3_TOMO_synapsin_dimer - Rab3_TOMO_p_synapsin_dimer_p
        Rab3_TOMO_synapsin_Rab3_TOMO_p_synapsin_p_dimerise_a = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin, Rab3_TOMO_p_synapsin_p), (Rab3_TOMO_synapsin_dimer, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer.v, None))
        Rab3_TOMO_synapsin_Rab3_TOMO_p_synapsin_p_dimerise_b = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin, Rab3_TOMO_p_synapsin_p), (Rab3_TOMO_synapsin_dimer, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer_p.v, None))
        Rab3_TOMO_synapsin_Rab3_TOMO_p_synapsin_p_dimerise_c = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin, Rab3_TOMO_p_synapsin_p), (Rab3_TOMO_synapsin_dimer, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer.v, None))
        Rab3_TOMO_synapsin_Rab3_TOMO_p_synapsin_p_dimerise_d = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin, Rab3_TOMO_p_synapsin_p), (Rab3_TOMO_synapsin_dimer, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer_p.v, None))
        Rab3_TOMO_synapsin_Rab3_TOMO_p_synapsin_p_dimerise_e = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin, Rab3_TOMO_p_synapsin_p), (Rab3_TOMO_synapsin_dimer, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer.v, None))
        Rab3_TOMO_synapsin_Rab3_TOMO_p_synapsin_p_dimerise_f = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin, Rab3_TOMO_p_synapsin_p), (Rab3_TOMO_synapsin_dimer, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer_p.v, None))
        
        Rab3_TOMO_synapsin_Rab3_TOMO_p_synapsin_p_unbind = VesicleUnbind.Create((ves, ves), (Rab3_TOMO_synapsin_dimer, Rab3_TOMO_p_synapsin_dimer_p), (Rab3_TOMO_synapsin, Rab3_TOMO_p_synapsin_p), kcst=koff_synapsin_dimer_p, immobilization=NO_EFFECT)
        
        #Rab3_TOMO_synapsin_p - Rab3_TOMO_synapsin_p -> Rab3_TOMO_synapsin_dimer_p - Rab3_TOMO_synapsin_dimer_p
        Rab3_TOMO_synapsin_p_Rab3_TOMO_synapsin_p_dimerise_a = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin_p, Rab3_TOMO_synapsin_p), (Rab3_TOMO_synapsin_dimer_p, Rab3_TOMO_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer.v, None))
        Rab3_TOMO_synapsin_p_Rab3_TOMO_synapsin_p_dimerise_b = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin_p, Rab3_TOMO_synapsin_p), (Rab3_TOMO_synapsin_dimer_p, Rab3_TOMO_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer_p.v, None))
        Rab3_TOMO_synapsin_p_Rab3_TOMO_synapsin_p_dimerise_c = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin_p, Rab3_TOMO_synapsin_p), (Rab3_TOMO_synapsin_dimer_p, Rab3_TOMO_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer.v, None))
        Rab3_TOMO_synapsin_p_Rab3_TOMO_synapsin_p_dimerise_d = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin_p, Rab3_TOMO_synapsin_p), (Rab3_TOMO_synapsin_dimer_p, Rab3_TOMO_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer_p.v, None))
        Rab3_TOMO_synapsin_p_Rab3_TOMO_synapsin_p_dimerise_e = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin_p, Rab3_TOMO_synapsin_p), (Rab3_TOMO_synapsin_dimer_p, Rab3_TOMO_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer.v, None))
        Rab3_TOMO_synapsin_p_Rab3_TOMO_synapsin_p_dimerise_f = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin_p, Rab3_TOMO_synapsin_p), (Rab3_TOMO_synapsin_dimer_p, Rab3_TOMO_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer_p.v, None))
        
        Rab3_TOMO_synapsin_p_Rab3_TOMO_synapsin_p_unbind = VesicleUnbind.Create((ves, ves), (Rab3_TOMO_synapsin_dimer_p, Rab3_TOMO_synapsin_dimer_p), (Rab3_TOMO_synapsin_p, Rab3_TOMO_synapsin_p), kcst=koff_synapsin_dimer_p, immobilization=NO_EFFECT)
        
        #Rab3_TOMO_synapsin_p - Rab3_TOMO_p_synapsin -> Rab3_TOMO_synapsin_dimer_p - Rab3_TOMO_p_synapsin_dimer
        Rab3_TOMO_synapsin_p_Rab3_TOMO_p_synapsin_dimerise_a = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin_p, Rab3_TOMO_p_synapsin), (Rab3_TOMO_synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer.v, None))
        Rab3_TOMO_synapsin_p_Rab3_TOMO_p_synapsin_dimerise_b = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin_p, Rab3_TOMO_p_synapsin), (Rab3_TOMO_synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer_p.v, None))
        Rab3_TOMO_synapsin_p_Rab3_TOMO_p_synapsin_dimerise_c = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin_p, Rab3_TOMO_p_synapsin), (Rab3_TOMO_synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer.v, None))
        Rab3_TOMO_synapsin_p_Rab3_TOMO_p_synapsin_dimerise_d = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin_p, Rab3_TOMO_p_synapsin), (Rab3_TOMO_synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer_p.v, None))
        Rab3_TOMO_synapsin_p_Rab3_TOMO_p_synapsin_dimerise_e = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin_p, Rab3_TOMO_p_synapsin), (Rab3_TOMO_synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer.v, None))
        Rab3_TOMO_synapsin_p_Rab3_TOMO_p_synapsin_dimerise_f = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin_p, Rab3_TOMO_p_synapsin), (Rab3_TOMO_synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer_p.v, None))
        
        Rab3_TOMO_synapsin_p_Rab3_TOMO_p_synapsin_unbind = VesicleUnbind.Create((ves, ves), (Rab3_TOMO_synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer), (Rab3_TOMO_synapsin_p, Rab3_TOMO_p_synapsin), kcst=koff_synapsin_dimer_p, immobilization=NO_EFFECT)
        
        #Rab3_TOMO_synapsin_p - Rab3_TOMO_p_synapsin_p -> Rab3_TOMO_synapsin_dimer_p - Rab3_TOMO_p_synapsin_dimer_p
        Rab3_TOMO_synapsin_p_Rab3_TOMO_p_synapsin_p_dimerise_a = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin_p, Rab3_TOMO_p_synapsin_p), (Rab3_TOMO_synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer.v, None))
        Rab3_TOMO_synapsin_p_Rab3_TOMO_p_synapsin_p_dimerise_b = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin_p, Rab3_TOMO_p_synapsin_p), (Rab3_TOMO_synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer_p.v, None))
        Rab3_TOMO_synapsin_p_Rab3_TOMO_p_synapsin_p_dimerise_c = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin_p, Rab3_TOMO_p_synapsin_p), (Rab3_TOMO_synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer.v, None))
        Rab3_TOMO_synapsin_p_Rab3_TOMO_p_synapsin_p_dimerise_d = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin_p, Rab3_TOMO_p_synapsin_p), (Rab3_TOMO_synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer_p.v, None))
        Rab3_TOMO_synapsin_p_Rab3_TOMO_p_synapsin_p_dimerise_e = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin_p, Rab3_TOMO_p_synapsin_p), (Rab3_TOMO_synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer.v, None))
        Rab3_TOMO_synapsin_p_Rab3_TOMO_p_synapsin_p_dimerise_f = VesicleBind.Create((ves, ves), (Rab3_TOMO_synapsin_p, Rab3_TOMO_p_synapsin_p), (Rab3_TOMO_synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer_p.v, None))
        
        Rab3_TOMO_synapsin_p_Rab3_TOMO_p_synapsin_p_unbind = VesicleUnbind.Create((ves, ves), (Rab3_TOMO_synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer_p), (Rab3_TOMO_synapsin_p, Rab3_TOMO_p_synapsin_p), kcst=koff_synapsin_dimer_p, immobilization=NO_EFFECT)
        
        #Rab3_TOMO_p_synapsin - Rab3_TOMO_p_synapsin -> Rab3_TOMO_p_synapsin_dimer - Rab3_TOMO_p_synapsin_dimer
        Rab3_TOMO_p_synapsin_Rab3_TOMO_p_synapsin_dimerise_a = VesicleBind.Create((ves, ves), (Rab3_TOMO_p_synapsin, Rab3_TOMO_p_synapsin), (Rab3_TOMO_p_synapsin_dimer, Rab3_TOMO_p_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer.v, None))
        Rab3_TOMO_p_synapsin_Rab3_TOMO_p_synapsin_dimerise_b = VesicleBind.Create((ves, ves), (Rab3_TOMO_p_synapsin, Rab3_TOMO_p_synapsin), (Rab3_TOMO_p_synapsin_dimer, Rab3_TOMO_p_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer_p.v, None))
        Rab3_TOMO_p_synapsin_Rab3_TOMO_p_synapsin_dimerise_c = VesicleBind.Create((ves, ves), (Rab3_TOMO_p_synapsin, Rab3_TOMO_p_synapsin), (Rab3_TOMO_p_synapsin_dimer, Rab3_TOMO_p_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer.v, None))
        Rab3_TOMO_p_synapsin_Rab3_TOMO_p_synapsin_dimerise_d = VesicleBind.Create((ves, ves), (Rab3_TOMO_p_synapsin, Rab3_TOMO_p_synapsin), (Rab3_TOMO_p_synapsin_dimer, Rab3_TOMO_p_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer_p.v, None))
        Rab3_TOMO_p_synapsin_Rab3_TOMO_p_synapsin_dimerise_e = VesicleBind.Create((ves, ves), (Rab3_TOMO_p_synapsin, Rab3_TOMO_p_synapsin), (Rab3_TOMO_p_synapsin_dimer, Rab3_TOMO_p_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer.v, None))
        Rab3_TOMO_p_synapsin_Rab3_TOMO_p_synapsin_dimerise_f = VesicleBind.Create((ves, ves), (Rab3_TOMO_p_synapsin, Rab3_TOMO_p_synapsin), (Rab3_TOMO_p_synapsin_dimer, Rab3_TOMO_p_synapsin_dimer), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer_p.v, None))
        
        Rab3_TOMO_p_synapsin_Rab3_TOMO_p_synapsin_unbind = VesicleUnbind.Create((ves, ves), (Rab3_TOMO_p_synapsin_dimer, Rab3_TOMO_p_synapsin_dimer), (Rab3_TOMO_p_synapsin, Rab3_TOMO_p_synapsin), kcst=koff_synapsin_dimer, immobilization=NO_EFFECT)
        
        #Rab3_TOMO_p_synapsin - Rab3_TOMO_p_synapsin_p -> Rab3_TOMO_p_synapsin_dimer - Rab3_TOMO_p_synapsin_dimer_p
        Rab3_TOMO_p_synapsin_Rab3_TOMO_p_synapsin_p_dimerise_a = VesicleBind.Create((ves, ves), (Rab3_TOMO_p_synapsin, Rab3_TOMO_p_synapsin_p), (Rab3_TOMO_p_synapsin_dimer, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer.v, None))
        Rab3_TOMO_p_synapsin_Rab3_TOMO_p_synapsin_p_dimerise_b = VesicleBind.Create((ves, ves), (Rab3_TOMO_p_synapsin, Rab3_TOMO_p_synapsin_p), (Rab3_TOMO_p_synapsin_dimer, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer_p.v, None))
        Rab3_TOMO_p_synapsin_Rab3_TOMO_p_synapsin_p_dimerise_c = VesicleBind.Create((ves, ves), (Rab3_TOMO_p_synapsin, Rab3_TOMO_p_synapsin_p), (Rab3_TOMO_p_synapsin_dimer, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer.v, None))
        Rab3_TOMO_p_synapsin_Rab3_TOMO_p_synapsin_p_dimerise_d = VesicleBind.Create((ves, ves), (Rab3_TOMO_p_synapsin, Rab3_TOMO_p_synapsin_p), (Rab3_TOMO_p_synapsin_dimer, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer_p.v, None))
        Rab3_TOMO_p_synapsin_Rab3_TOMO_p_synapsin_p_dimerise_e = VesicleBind.Create((ves, ves), (Rab3_TOMO_p_synapsin, Rab3_TOMO_p_synapsin_p), (Rab3_TOMO_p_synapsin_dimer, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer.v, None))
        Rab3_TOMO_p_synapsin_Rab3_TOMO_p_synapsin_p_dimerise_f = VesicleBind.Create((ves, ves), (Rab3_TOMO_p_synapsin, Rab3_TOMO_p_synapsin_p), (Rab3_TOMO_p_synapsin_dimer, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer_p.v, None))
        
        Rab3_TOMO_p_synapsin_Rab3_TOMO_p_synapsin_p_unbind = VesicleUnbind.Create((ves, ves), (Rab3_TOMO_p_synapsin_dimer, Rab3_TOMO_p_synapsin_dimer_p), (Rab3_TOMO_p_synapsin, Rab3_TOMO_p_synapsin_p), kcst=koff_synapsin_dimer_p, immobilization=NO_EFFECT)
        
        #Rab3_TOMO_p_synapsin_p - Rab3_TOMO_p_synapsin_p -> Rab3_TOMO_p_synapsin_dimer_p - Rab3_TOMO_p_synapsin_dimer_p
        Rab3_TOMO_p_synapsin_p_Rab3_TOMO_p_synapsin_p_dimerise_a = VesicleBind.Create((ves, ves), (Rab3_TOMO_p_synapsin_p, Rab3_TOMO_p_synapsin_p), (Rab3_TOMO_p_synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer.v, None))
        Rab3_TOMO_p_synapsin_p_Rab3_TOMO_p_synapsin_p_dimerise_b = VesicleBind.Create((ves, ves), (Rab3_TOMO_p_synapsin_p, Rab3_TOMO_p_synapsin_p), (Rab3_TOMO_p_synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(synapsin_dimer_p.v, None))
        Rab3_TOMO_p_synapsin_p_Rab3_TOMO_p_synapsin_p_dimerise_c = VesicleBind.Create((ves, ves), (Rab3_TOMO_p_synapsin_p, Rab3_TOMO_p_synapsin_p), (Rab3_TOMO_p_synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer.v, None))
        Rab3_TOMO_p_synapsin_p_Rab3_TOMO_p_synapsin_p_dimerise_d = VesicleBind.Create((ves, ves), (Rab3_TOMO_p_synapsin_p, Rab3_TOMO_p_synapsin_p), (Rab3_TOMO_p_synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_synapsin_dimer_p.v, None))
        Rab3_TOMO_p_synapsin_p_Rab3_TOMO_p_synapsin_p_dimerise_e = VesicleBind.Create((ves, ves), (Rab3_TOMO_p_synapsin_p, Rab3_TOMO_p_synapsin_p), (Rab3_TOMO_p_synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer.v, None))
        Rab3_TOMO_p_synapsin_p_Rab3_TOMO_p_synapsin_p_dimerise_f = VesicleBind.Create((ves, ves), (Rab3_TOMO_p_synapsin_p, Rab3_TOMO_p_synapsin_p), (Rab3_TOMO_p_synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer_p), 0, 1e-07, kcst=kon_synapsin_dimer, immobilization=NO_EFFECT, deps=(Rab3_TOMO_p_synapsin_dimer_p.v, None))
        
        Rab3_TOMO_p_synapsin_p_Rab3_TOMO_p_synapsin_p_unbind = VesicleUnbind.Create((ves, ves), (Rab3_TOMO_p_synapsin_dimer_p, Rab3_TOMO_p_synapsin_dimer_p), (Rab3_TOMO_p_synapsin_p, Rab3_TOMO_p_synapsin_p), kcst=koff_synapsin_dimer_p, immobilization=NO_EFFECT)

    with memsys:
        diff_PMCA_P0 = Diffusion.Create(PMCA_P0, DCST_MEM)
        diff_PMCA_P1 = Diffusion.Create(PMCA_P1, DCST_MEM)

        # Ca + PMCA <->  CaPMCA -> PMCA
        Ca.i + PMCA_P0.s <r[1]> PMCA_P1.s >r[2]> PMCA_P0.s
        r[1].K = k1_pmca, k2_pmca
        r[2].K = k3_pmca

        # PMCA leak
        PMCA_P0.s >r[1]> Ca.i + PMCA_P0.s
        r[1].K = kl_pmca

        diff_syt = Diffusion.Create(syt, 2.5e-13)
        sdiff_glu = Diffusion.Create(glu, 0.0)

        # CaP channel
        if not cachan_static:
            Diffusion(CaPchan[...], 0.01e-12)
        with CaPchan[...]:
            CaPc.s <r[1]> CaPo.s
            r[1].K = alpha_cap, beta_cap
        OC_CaP = GHKCurr.Create(
            CaPchan[CaPo, CaPo, CaPo], Ca, CaP_P,
            computeflux=True,
            virtual_oconc=Ca_oconc,
        )

        # Things will get transported by endocytosis- so even though they will be clamped to 0 in some places need to define them in patches
        diff_RIM = Diffusion.Create(RIM, 0.0)
        sdiff_M13 = Diffusion.Create(M13, 0.0)
        diff_RIM_M13 = Diffusion.Create(RIM_M13, 0.0)
        sdiff_Rab3 = Diffusion.Create(Rab3, 2.5e-13)
        diff_RIM_M13_Rab3 = Diffusion.Create(RIM_M13_Rab3, 0.0)
        diff_SYX = Diffusion.Create(SYX, 2.429e-13)
        diff_SYX_M18 = Diffusion.Create(SYX_M18, 2.429e-13)
        sdiff_M18 = Diffusion.Create(M18, 0.0)
        diff_RIM_M13_Rab3_SYX_M18 = Diffusion.Create(RIM_M13_Rab3_SYX_M18, 0.0)
        diff_RIM_M13_Rab3_SYX_M18_SNP25  = Diffusion.Create(RIM_M13_Rab3_SYX_M18_SNP25, 0.0)
        diff_SNP25 = Diffusion.Create(SNP25, 6.51e-13)
        diff_SYB = Diffusion.Create(SYB, 2.5e-13)
        diff_SNARE = Diffusion.Create(SNARE, 0.0)
        sdiff_CXN = Diffusion.Create(CXN, 0.0)

        #######RIM binds and activates Munc13 homodimer
        kon_rim_m13 = 10e06
        koff_rim_m13 = 0

        kon_syx_m18 = 5e06 #Burkhardt 2008
        koff_syx_m18 = 0.0011 #Burkhardt 2008

        #reactions
        M13.i + RIM.s <r[1]> RIM_M13.s
        r[1].K = kon_rim_m13, koff_rim_m13

        ##Blocking of Syntaxin-1 by MUNC18a
        M18.i + SYX.s <r[1]> SYX_M18.s
        r[1].K = kon_syx_m18, koff_syx_m18

        kon_snare_snap = 1.7e05
        koff_snare_snap = 0.26
        kon_nsf_snap = 20e06
        kcat_nsf = 0.116 #Vivona 2013
        # #initial loss of CXN, Ca,  syt, etc from all SNARE variants
        k_rim_diss = 1000e12

        SNARE.s >r[1]> Rab3.i + RIM_M13.s + cisSNARE.s
        r[1].K = k_rim_diss
        SNARE_syt_CXN_Ca3_bCa2.s >r[1]> CXN.i + 5 * Ca.i + RIM_M13.s + Rab3.s + cisSNARE.s + syt.s
        r[1].K = k_rim_diss
        SNARE_syt_CXN_Ca3_bCa.s >r[1]> CXN.i + 4 * Ca.i + RIM_M13.s + Rab3.s + cisSNARE.s + syt.s
        r[1].K = k_rim_diss
        SNARE_syt_CXN_Ca2_bCa2.s >r[1]> CXN.i + 4 * Ca.i + RIM_M13.s + Rab3.s + cisSNARE.s + syt.s
        r[1].K = k_rim_diss
        SNARE_syt_CXN_Ca2_bCa.s >r[1]> CXN.i + 3 * Ca.i + RIM_M13.s + Rab3.s + cisSNARE.s + syt.s
        r[1].K = k_rim_diss
        SNARE_syt_CXN_Ca_bCa2.s >r[1]> CXN.i + 3 * Ca.i + RIM_M13.s + Rab3.s + cisSNARE.s + syt.s
        r[1].K = k_rim_diss
        SNARE_syt_CXN_Ca_bCa.s >r[1]> CXN.i + 2 * Ca.i + RIM_M13.s + Rab3.s + cisSNARE.s + syt.s
        r[1].K = k_rim_diss
        SNARE_syt_CXN_bCa2.s >r[1]> CXN.i + 2 * Ca.i + RIM_M13.s + Rab3.s + cisSNARE.s + syt.s
        r[1].K = k_rim_diss
        SNARE_syt_CXN_bCa.s >r[1]> CXN.i + Ca.i + RIM_M13.s + Rab3.s + cisSNARE.s + syt.s
        r[1].K = k_rim_diss
        SNARE_syt_CXN_Ca3.s >r[1]> CXN.i + 3 * Ca.i + RIM_M13.s + Rab3.s + cisSNARE.s + syt.s
        r[1].K = k_rim_diss
        SNARE_syt_CXN_Ca2.s >r[1]> CXN.i + 2 * Ca.i + RIM_M13.s + Rab3.s + cisSNARE.s + syt.s
        r[1].K = k_rim_diss
        SNARE_syt_CXN_Ca.s >r[1]> CXN.i + Ca.i + RIM_M13.s + Rab3.s + cisSNARE.s + syt.s
        r[1].K = k_rim_diss
        SNARE_syt_CXN.s >r[1]> CXN.i + RIM_M13.s + Rab3.s + cisSNARE.s + syt.s
        r[1].K = k_rim_diss
        SNARE_syt.s >r[1]> RIM_M13.s + Rab3.s + cisSNARE.s + syt.s
        r[1].K = k_rim_diss
        
        # # #disassembly of SNARE by NSF/aSNAP
        aSNAP.i + cisSNARE.s >r[1]> SNARE_aSNAP.s
        r[1].K = kon_snare_snap
        SNARE_aSNAP.s >r[1]> aSNAP.i + SNARE.s
        r[1].K = koff_snare_snap
        NSF.i + SNARE_aSNAP.s >r[1]> SNARE_aSNAP_NSF.s
        r[1].K = kon_nsf_snap
        SNARE_aSNAP_NSF.s >r[1]> M18.i + NSF.i + aSNAP.i + SNARE_DISS.s + SNP25.s + SYB.s + SYX.s
        r[1].K = kcat_nsf

        sdiff_DYNpp = Diffusion.Create(DYNpp, 3.2464e-12)
        sdiff_DYNp = Diffusion.Create(DYNp, 3.2464e-12)
        sdiff_DYN = Diffusion.Create(DYN, 3.2464e-12)
        sdiff_CaN_DYNpp = Diffusion.Create(CaN_DYNpp, 3.2464e-12)
        sdiff_CaN_DYNp = Diffusion.Create(CaN_DYNp, 3.2464e-12)

        sdiff_DYN_SYN1 = Diffusion.Create(DYN_SYN1, 3.2464e-12)

        kon_dyn_mem = 0
        koff_dyn_mem = 0
        DYN_SYN1.i <r[1]> DYN_SYN1.s
        r[1].K = kon_dyn_mem, koff_dyn_mem
        DYNp.i <r[1]> DYNp.s
        r[1].K = kon_dyn_mem, koff_dyn_mem
        DYN.i <r[1]> DYN.s
        r[1].K = kon_dyn_mem, koff_dyn_mem

        sdiff_DYN_SYN1_CDK5 = Diffusion.Create(DYN_SYN1_CDK5, DCST_MEM)

        #phosphorylation of dynamin in the membrane
        CDK5.i + DYN_SYN1.s <r[1]> DYN_SYN1_CDK5.s
        r[1].K = kon_cdk_dyn, koff_cdk_dyn
        DYN_SYN1_CDK5.s >r[1]> CDK5.i + DYNp.i + SYN1.i
        r[1].K = kcat_cdk
        
        CDK5.s + DYN.s <r[1]> DYN_CDK5.s
        r[1].K = kon_cdk_dyn, koff_cdk_dyn
        DYN_CDK5.s >r[1]> CDK5.s + DYNp.s
        r[1].K = kcat_cdk

        sdiff_AP2 = Diffusion.Create(AP2, 0.0)
        sdiff_RabAd = Diffusion.Create(RabAd, 0.0)
        sdiff_DYN_AD = Diffusion.Create(DYN_AD, 0.0)
        sdiff_AP180 = Diffusion.Create(AP180, 0.0)

        sdiff_synapsin = Diffusion.Create(synapsin, 0.0)
        sdiff_synapsin_p = Diffusion.Create(synapsin_p, 0.0)
        sdiff_synapsin_pka = Diffusion.Create(PKA_synapsin, 0.0)
        sdiff_synapsin_pp2a = Diffusion.Create(PP2A_synapsin_p, 0.0)
        sdiff_Rab3_TOMO_synapsin_actin = Diffusion.Create(Rab3_TOMO_synapsin_actin, 0.0)
        sdiff_Rab3_TOMO_p_synapsin_actin = Diffusion.Create(Rab3_TOMO_p_synapsin_actin, 0.0)
        sdiff_Rab3_TOMO_synapsin = Diffusion.Create(Rab3_TOMO_synapsin, 0.0)
        sdiff_Rab3_TOMO_synapsin_p_ves = Diffusion.Create(Rab3_TOMO_synapsin_p, 0.0)
        sdiff_Rab3_TOMO_p_synapsin = Diffusion.Create(Rab3_TOMO_p_synapsin, 0.0)
        sdiff_Rab3_TOMO_p_synapsin_p = Diffusion.Create(Rab3_TOMO_p_synapsin_p, 0.0)
        sdiff_Rab3_TOMO_p = Diffusion.Create(Rab3_TOMO_p, 0.0)
        sdiff_Rab3_TOMO = Diffusion.Create(Rab3_TOMO, 0.0)
        sdiff_Rab3_TOMO_CDK5 = Diffusion.Create(Rab3_TOMO_CDK5, 0.0)
        sdiff_Rab3_TOMO_p_CaN_CaM = Diffusion.Create(Rab3_TOMO_p_CaN_CaM, 0.0)
        sdiff_Rab3_TOMO_synapsin_PKA = Diffusion.Create(Rab3_TOMO_synapsin_PKA, 0.0)
        sdiff_Rab3_TOMO_p_synapsin_PKA = Diffusion.Create(Rab3_TOMO_p_synapsin_PKA, 0.0)
        sdiff_Rab3_TOMO_synapsin_CDK5 = Diffusion.Create(Rab3_TOMO_synapsin_CDK5, 0.0)
        sdiff_Rab3_TOMO_synapsin_p_CDK5 = Diffusion.Create(Rab3_TOMO_synapsin_p_CDK5, 0.0)
        sdiff_Rab3_TOMO_p_synapsin_CaN_CaM = Diffusion.Create(Rab3_TOMO_p_synapsin_CaN_CaM, 0.0)
        sdiff_Rab3_TOMO_p_synapsin_p_CaN_CaM = Diffusion.Create(Rab3_TOMO_p_synapsin_p_CaN_CaM, 0.0)
        sdiff_RIM_M13_Rab3_SYX_M18_SNP25_TOMO = Diffusion.Create(RIM_M13_Rab3_SYX_M18_SNP25_TOMO, 0)
        sdiff_RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p = Diffusion.Create(RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p, 0)
        sdiff_RIM_M13_Rab3_SYX_M18_SNP25_TOMOx = Diffusion.Create(RIM_M13_Rab3_SYX_M18_SNP25_TOMOx, 0)
        sdiff_RIM_M13_Rab3_SYX_M18_SNP25_TOMOx_p = Diffusion.Create(RIM_M13_Rab3_SYX_M18_SNP25_TOMOx_p, 0)
        sdiff_RIM_M13_Rab3_SYX_M18_SNP25_TOMO_P = Diffusion.Create(RIM_M13_Rab3_SYX_M18_SNP25_TOMO_P, 0)
        sdiff_RIM_M13_Rab3_SYX_M18_SNP25_TOMOx_P = Diffusion.Create(RIM_M13_Rab3_SYX_M18_SNP25_TOMOx_P, 0)
        sdiff_RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p_P = Diffusion.Create(RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p_P, 0)
        sdiff_RIM_M13_Rab3_SYX_M18_SNP25_TOMOx_p_P = Diffusion.Create(RIM_M13_Rab3_SYX_M18_SNP25_TOMOx_p_P, 0)
        sdiff_RIM_M13_Rab3_SYX_M18_SNP25_TOMO_PKA = Diffusion.Create(RIM_M13_Rab3_SYX_M18_SNP25_TOMO_PKA, 0)
        sdiff_RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p_PKA = Diffusion.Create(RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p_PKA, 0)
        sdiff_RIM_M13_Rab3_SYX_M18_SNP25_TOMO_P_PP2A = Diffusion.Create(RIM_M13_Rab3_SYX_M18_SNP25_TOMO_P_PP2A, 0)
        sdiff_RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p_P_PP2A = Diffusion.Create(RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p_P_PP2A, 0)
        sdiff_CXN_site = Diffusion.Create(ves_CXN_site, 0)
        sdiff_rab3_site = Diffusion.Create(ves_rab3_site, 0)
        sdiff_M13_site = Diffusion.Create(ves_m13_site, 0)
        sdiff_M18_site = Diffusion.Create(ves_m18_site, 0)
        sdiff_syndap_site = Diffusion.Create(ves_syndap_site, 0)
        sdiff_SYN1 = Diffusion.Create(SYN1, 0)
        sdiff_asnap_site = Diffusion.Create(ves_asnap_site, 0)
        sdiff_asnap = Diffusion.Create(aSNAP, 0)
        sdiff_nsf_site = Diffusion.Create(ves_nsf_site, 0)
        sdiff_NSF = Diffusion.Create(NSF, 0)

        #####Remove vesicle-protein sites and other proteins from plasma membrane after vesicle fusion
        ves_syn1_site.s >r[1]> None
        r[1].K = 1000000.0
        ves_CXN_site.s >r[1]> None
        r[1].K = 1000000.0
        ves_m13_site.s >r[1]> None
        r[1].K = 1000000.0
        ves_m18_site.s >r[1]> None
        r[1].K = 1000000.0
        ves_syndap_site.s >r[1]> None
        r[1].K = 1000000.0
        ves_asnap_site.s >r[1]> None
        r[1].K = 1000000.0
        ves_nsf_site.s >r[1]> None
        r[1].K = 1000000.0
        ves_rab3_site.s >r[1]> None
        r[1].K = 1000000.0
        NSF.s >r[1]> NSF.i
        r[1].K = 1000000.0
        
        synapsin.s >r[1]> synapsin.i
        r[1].K = 1000000.0
        synapsin_p.s >r[1]> synapsin.i
        r[1].K = 1000000.0
        PKA_synapsin.s >r[1]> PKA_synapsin.i
        r[1].K = 1000000.0
        PP2A_synapsin_p.s >r[1]> PP2A_synapsin_p.i
        r[1].K = 1000000.0
        
        TOMO_p.s >r[1]> TOMO_p.i
        r[1].K = 1000000.0
        Rab3_TOMO_p.s >r[1]> Rab3.i + TOMO_p.i
        r[1].K = 1000000.0
        Rab3_TOMO.s >r[1]> Rab3.i + TOMO.i
        r[1].K = 1000000.0
        Rab3_TOMO_p_synapsin.s >r[1]> Rab3.i + TOMO_p.i + synapsin.i
        r[1].K = 1000000.0
        Rab3_TOMO_synapsin.s >r[1]> Rab3.i + TOMO.i + synapsin.i
        r[1].K = 1000000.0
        Rab3_TOMO_p_synapsin_p.s >r[1]> Rab3.i + TOMO_p.i + synapsin_p.i
        r[1].K = 1000000.0
        Rab3_TOMO_synapsin_p.s >r[1]> Rab3.i + TOMO.i + synapsin_p.i
        r[1].K = 1000000.0
        
        
        actin.s >r[1]> actin.i
        r[1].K = 1000000.0
        synapsin_actin.s >r[1]> actin.i + synapsin.i
        r[1].K = 1000000.0
        Rab3_TOMO_synapsin_actin.s >r[1]> Rab3.i + TOMO.i + actin.i + synapsin.i
        r[1].K = 1000000.0
        Rab3_TOMO_p_synapsin_actin.s >r[1]> Rab3.i + TOMO_p.i + actin.i + synapsin.i
        r[1].K = 1000000.0
        Rab3_TOMO_synapsin_PKA.s >r[1]> PKA.i + Rab3.i + TOMO.i + synapsin.i
        r[1].K = 1000000.0
        Rab3_TOMO_p_synapsin_PKA.s >r[1]> PKA.i + Rab3.i + TOMO_p.i + synapsin.i
        r[1].K = 1000000.0
        synapsin_actin_p.s >r[1]> actin.i + synapsin.i
        r[1].K = 1000000.0
        Rab3_TOMO_p_CaN_CaM.s >r[1]> CaN_CaM.i + Rab3.i + TOMO_p.i
        r[1].K = 1000000.0
        Rab3_TOMO_p_synapsin_CaN_CaM.s >r[1]> CaN_CaM.i + Rab3.i + TOMO_p.i + synapsin.i
        r[1].K = 1000000.0
        Rab3_TOMO_p_synapsin_p_CaN_CaM.s >r[1]> CaN_CaM.i + Rab3.i + TOMO_p.i + synapsin_p.i
        r[1].K = 1000000.0
        Rab3_TOMO_CDK5.s >r[1]> CDK5.i + Rab3.i + TOMO.i
        r[1].K = 1000000.0
        Rab3_TOMO_synapsin_CDK5.s >r[1]> CDK5.i + Rab3.i + TOMO.i + synapsin.i
        r[1].K = 1000000.0
        Rab3_TOMO_synapsin_p_CDK5.s >r[1]> CDK5.i + Rab3.i + TOMO.i + synapsin_p.i
        r[1].K = 1000000.0
        Rab3_TOMO_p_synapsin_p_PP2A.s >r[1]> PP2A.i + Rab3.i + TOMO_p.i + synapsin_p.i
        r[1].K = 1000000.0
        Rab3_TOMO_synapsin_p_PP2A.s >r[1]> PP2A.i + Rab3.i + TOMO.i + synapsin_p.i
        r[1].K = 1000000.0
        
        RIM_M13_Rab3_SYX_M18_SNP25_TOMO.s >r[1]> M18.i + Rab3.i + TOMO.i + RIM_M13.s + SNP25.s + SYX.s
        r[1].K = 1000000.0
        RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p.s >r[1]> M18.i + Rab3.i + TOMO_p.i + RIM_M13.s + SNP25.s + SYX.s
        r[1].K = 1000000.0
        RIM_M13_Rab3_SYX_M18_SNP25_TOMOx.s >r[1]> M18.i + Rab3.i + TOMO.i + RIM_M13.s + SNP25.s + SYX.s
        r[1].K = 1000000.0
        RIM_M13_Rab3_SYX_M18_SNP25_TOMOx_p.s >r[1]> M18.i + Rab3.i + TOMO_p.i + RIM_M13.s + SNP25.s + SYX.s
        r[1].K = 1000000.0
        RIM_M13_Rab3_SYX_M18_SNP25_TOMO_P.s >r[1]> M18.i + Rab3.i + TOMO.i + RIM_M13.s + SNP25.s + SYX.s
        r[1].K = 1000000.0
        RIM_M13_Rab3_SYX_M18_SNP25_TOMOx_P.s >r[1]> M18.i + Rab3.i + TOMO.i + RIM_M13.s + SNP25.s + SYX.s
        r[1].K = 1000000.0
        RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p_P.s >r[1]> M18.i + Rab3.i + TOMO_p.i + RIM_M13.s + SNP25.s + SYX.s
        r[1].K = 1000000.0
        RIM_M13_Rab3_SYX_M18_SNP25_TOMOx_p_P.s >r[1]> M18.i + Rab3.i + TOMO_p.i + RIM_M13.s + SNP25.s + SYX.s
        r[1].K = 1000000.0
        RIM_M13_Rab3_SYX_M18_SNP25_TOMO_PKA.s >r[1]> M18.i + PKA.i + Rab3.i + TOMO.i + RIM_M13.s + SNP25.s + SYX.s
        r[1].K = 1000000.0
        RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p_PKA.s >r[1]> M18.i + PKA.i + Rab3.i + TOMO_p.i + RIM_M13.s + SNP25.s + SYX.s
        r[1].K = 1000000.0
        RIM_M13_Rab3_SYX_M18_SNP25_TOMO_P_PP2A.s >r[1]> M18.i + PP2A.i + Rab3.i + TOMO.i + RIM_M13.s + SNP25.s + SYX.s
        r[1].K = 1000000.0
        RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p_P_PP2A.s >r[1]> M18.i + PP2A.i + Rab3.i + TOMO_p.i + RIM_M13.s + SNP25.s + SYX.s
        r[1].K = 1000000.0



    with ERsys:
        diff_SERCA = Diffusion.Create(SERCA, DCST_MEM)
        diff_CaSERCA = Diffusion.Create(CaSERCA, DCST_MEM)
        diff_Ca2SERCA = Diffusion.Create(Ca2SERCA, DCST_MEM)

        # Ca + SERCA <->  Ca1SERCA +Ca <->  Ca2SERCA  ->  SERCA
        SERCA.s + Ca.o <r[1]> CaSERCA.s
        CaSERCA.s + Ca.o <r[2]> Ca2SERCA.s >r[3]> 2 * Ca.i + SERCA.s
        r[1].K = 17147000000.0, 8426.3
        r[2].K = 17147000000.0, 8426.3
        r[3].K = 250.0

    with syn_sys:
        #Diffusion of glutamate in the cleft
        diff_glu = Diffusion.Create(glu, DCST_CYT)
        diff_REL_IND = Diffusion.Create(REL_IND, 0)

    
    
    #######Synaptotagmin Model#######
    #Parameters
    kon_SNARE_syt_CXN_ca1 = 1000e06
    koff_SNARE_syt_CXN_ca1 = 1270
    koff_SNARE_syt_CXN_ca1_s = 1000 #Millet 2002
    
    kon_SNARE_syt_CXN_ca2 = 1000e06
    koff_SNARE_syt_CXN_ca2 = 227670
    koff_SNARE_syt_CXN_ca2_s = 2000
    
    kon_SNARE_syt_CXN_ca3 = 1000e06
    koff_SNARE_syt_CXN_ca3 = 12370
    koff_SNARE_syt_CXN_ca3_s = 5000
    
    kon_SNARE_syt_CXN_bca1 = 1000e06
    koff_SNARE_syt_CXN_bca1 = 50000
    
    kon_SNARE_syt_CXN_bca2 = 1000e06
    koff_SNARE_syt_CXN_bca2 = 25780
    
    kon_SNARE_syt_CXN_mem  = 2500
    koff_SNARE_syt_CXN_mem = 700

    with ves_ssys:
        #vesicle Diffusion of Synaptotagmin
        vdiff_syt = Diffusion.Create(syt, DCST_MEM_VES)
        #Reactions
        Ca.o + SNARE_syt_CXN.v <r[1]> SNARE_syt_CXN_Ca.v
        r[1].K = kon_SNARE_syt_CXN_ca1, koff_SNARE_syt_CXN_ca1
        
        Ca.o + SNARE_syt_CXN_Ca.v <r[1]> SNARE_syt_CXN_Ca2.v
        r[1].K = kon_SNARE_syt_CXN_ca2, koff_SNARE_syt_CXN_ca2
        
        Ca.o + SNARE_syt_CXN_Ca2.v <r[1]> SNARE_syt_CXN_Ca3.v
        r[1].K = kon_SNARE_syt_CXN_ca3, koff_SNARE_syt_CXN_ca3
        
        Ca.o + SNARE_syt_CXN.v <r[1]> SNARE_syt_CXN_bCa.v
        r[1].K = kon_SNARE_syt_CXN_bca1, koff_SNARE_syt_CXN_bca1
        
        Ca.o + SNARE_syt_CXN_bCa.v <r[1]> SNARE_syt_CXN_bCa2.v
        r[1].K = kon_SNARE_syt_CXN_bca2, koff_SNARE_syt_CXN_bca2
        
        Ca.o + SNARE_syt_CXN_bCa.v <r[1]> SNARE_syt_CXN_Ca_bCa.v
        r[1].K = kon_SNARE_syt_CXN_ca1, koff_SNARE_syt_CXN_ca1
        
        Ca.o + SNARE_syt_CXN_Ca_bCa.v <r[1]> SNARE_syt_CXN_Ca2_bCa.v
        r[1].K = kon_SNARE_syt_CXN_ca2, koff_SNARE_syt_CXN_ca2
        
        Ca.o + SNARE_syt_CXN_Ca2_bCa.v <r[1]> SNARE_syt_CXN_Ca3_bCa.v
        r[1].K = kon_SNARE_syt_CXN_ca3, koff_SNARE_syt_CXN_ca3
        
        Ca.o + SNARE_syt_CXN_Ca.v <r[1]> SNARE_syt_CXN_Ca_bCa.v
        r[1].K = kon_SNARE_syt_CXN_bca1, koff_SNARE_syt_CXN_bca1
        
        Ca.o + SNARE_syt_CXN_Ca_bCa.v <r[1]> SNARE_syt_CXN_Ca_bCa2.v
        r[1].K = kon_SNARE_syt_CXN_bca2, koff_SNARE_syt_CXN_bca2
        
        Ca.o + SNARE_syt_CXN_Ca2.v <r[1]> SNARE_syt_CXN_Ca2_bCa.v
        r[1].K = kon_SNARE_syt_CXN_bca1, koff_SNARE_syt_CXN_bca1
        
        Ca.o + SNARE_syt_CXN_Ca2_bCa.v <r[1]> SNARE_syt_CXN_Ca2_bCa2.v
        r[1].K = kon_SNARE_syt_CXN_bca2, koff_SNARE_syt_CXN_bca2
        
        Ca.o + SNARE_syt_CXN_Ca3.v <r[1]> SNARE_syt_CXN_Ca3_bCa.v
        r[1].K = kon_SNARE_syt_CXN_bca1, koff_SNARE_syt_CXN_bca1
        
        Ca.o + SNARE_syt_CXN_Ca3_bCa.v <r[1]> SNARE_syt_CXN_Ca3_bCa2.v
        r[1].K = kon_SNARE_syt_CXN_bca2, koff_SNARE_syt_CXN_bca2
        
        Ca.o + SNARE_syt_CXN_bCa2.v <r[1]> SNARE_syt_CXN_Ca_bCa2.v
        r[1].K = kon_SNARE_syt_CXN_ca1, koff_SNARE_syt_CXN_ca1_s
        
        Ca.o + SNARE_syt_CXN_Ca_bCa2.v <r[1]> SNARE_syt_CXN_Ca2_bCa2.v
        r[1].K = kon_SNARE_syt_CXN_ca2, koff_SNARE_syt_CXN_ca2_s
        
        Ca.o + SNARE_syt_CXN_Ca2_bCa2.v <r[1]> SNARE_syt_CXN_Ca3_bCa2.v
        r[1].K = kon_SNARE_syt_CXN_ca3, koff_SNARE_syt_CXN_ca3_s

        kon_m13_syx = 20e06 #Mezer 2004
        koff_m13_syx = 2.6 #Mezer 2004
        kon_syx_snp25 = 10e06
        koff_syx_snp25 = 1.26
        kon_syb_syx = 2.35e06
        koff_syb_syx = 0.0047
        #######Vesicle docks by Rab3 via RIM_M13 interaction
        kon_rim_rab3 = 1000e9
        koff_rim_rab3 = 0#15

        kon_syt_snare = 57e06 #32e06 #Mezer 2004
        koff_syt_snare = 0.24 #9.4 #Mezer 2004

        kon_cxn_snare = 30e06 #7e06 #Mezer 2004
        koff_cxn_snare = 0.33 #Mezer 2004

        vdiff_Rab3 = Diffusion.Create(Rab3, 5e-13) #estimated from similar proteins
        #vesicle Diffusion of Synaptobrevin
        vdiff_syb = Diffusion.Create(SYB, 3.245e-13)

        #reactions
        mdist=10e-09

        RIM_M13.s + Rab3.v <r[1]> RIM_M13_Rab3.v
        r[1].K = kon_rim_rab3, koff_rim_rab3
        r[1].Immobilization = IMMOBILIZING, MOBILIZING
        r[1].MaxDistance = mdist, None

        ##Munc13 (active form in complex with RIM)  binds and opens syntaxin-1 (which will allow it to bind SNAP25) (Lai 2017)
        SYX_M18.s + RIM_M13_Rab3.v <r[1]> RIM_M13_Rab3_SYX_M18.v
        r[1].K = kon_m13_syx, koff_m13_syx
        
        #######SNAP25 is incorporated into the SNARE complex (binds to SYX) (Lai 2017)
        SNP25.s + RIM_M13_Rab3_SYX_M18.v <r[1]> RIM_M13_Rab3_SYX_M18_SNP25.v
        r[1].K = kon_syx_snp25, koff_syx_snp25
        
        #######Synaptobrevin is incorporated into the SNARE complex (binds to SYX_SNP25 complex)
        RIM_M13_Rab3_SYX_M18_SNP25.v + SYB.v <r[1]> SNARE.v
        r[1].K = kon_syb_syx, koff_syb_syx

        ##Binding of syt to SNARE
        SNARE.v + syt.v <r[1]> SNARE_syt.v
        r[1].K = kon_syt_snare, koff_syt_snare
        
        ##Binding of CXN to SNARE
        CXN.o + SNARE_syt.v <r[1]> SNARE_syt_CXN.v
        r[1].K = kon_cxn_snare, koff_cxn_snare

        ##VESICLE FUSION EVENT (EXOCYTOSIS)
        fusion2 = Exocytosis.Create(3000.0, deps=2 * SNARE_syt_CXN_Ca3_bCa2.v)
        fusion3 = Exocytosis.Create(30000.0, deps=3 * SNARE_syt_CXN_Ca3_bCa2.v)
        fusion4 = Exocytosis.Create(300000.0, deps=4 * SNARE_syt_CXN_Ca3_bCa2.v)

        #removal of DYN_SYN1 from newly endocytosed vesicle
        DYN_SYN1.v >r[1]> DYNp.o + SYN1.o
        r[1].K = 1000000000.0

        diff_synapsin_mem = Diffusion.Create(synapsin, DCST_MEM_VES)
        diff_synapsin_p_mem = Diffusion.Create(synapsin_p, DCST_MEM_VES)
        diff_synapsin_pka_mem = Diffusion.Create(PKA_synapsin, DCST_MEM_VES)
        diff_synapsin_pp2a_ves = Diffusion.Create(PP2A_synapsin_p, DCST_MEM_VES)

        actin.o + synapsin.v <r[1]> synapsin_actin.v
        r[1].K = kon_synapsin_actin, koff_synapsin_actin
        r[1].Immobilization = IMMOBILIZING, MOBILIZING
        actin.o + synapsin_p.v <r[1]> synapsin_actin_p.v
        r[1].K = kon_synapsin_actin, koff_synapsin_actin_p
        r[1].Immobilization = IMMOBILIZING, MOBILIZING

        diff_Rab3_TOMO_synapsin_ves = Diffusion.Create(Rab3_TOMO_synapsin, DCST_MEM_VES)
        diff_Rab3_TOMO_synapsin_p_ves = Diffusion.Create(Rab3_TOMO_synapsin_p, DCST_MEM_VES)
        diff_Rab3_TOMO_p_synapsin_ves = Diffusion.Create(Rab3_TOMO_p_synapsin, DCST_MEM_VES)
        diff_Rab3_TOMO_p_synapsin_p_ves = Diffusion.Create(Rab3_TOMO_p_synapsin_p, DCST_MEM_VES)
        diff_Rab3_TOMO_p_ves = Diffusion.Create(Rab3_TOMO_p, DCST_MEM_VES)
        diff_Rab3_TOMO_ves = Diffusion.Create(Rab3_TOMO, DCST_MEM_VES)
        diff_Rab3_TOMO_CDK5_ves = Diffusion.Create(Rab3_TOMO_CDK5, 5e-13)
        diff_Rab3_TOMO_p_CaN_CaM_ves = Diffusion.Create(Rab3_TOMO_p_CaN_CaM, DCST_MEM_VES)
        diff_Rab3_TOMO_synapsin_PKA_ves = Diffusion.Create(Rab3_TOMO_synapsin_PKA, DCST_MEM_VES)
        diff_Rab3_TOMO_p_synapsin_PKA_ves = Diffusion.Create(Rab3_TOMO_p_synapsin_PKA, DCST_MEM_VES)
        diff_Rab3_TOMO_synapsin_CDK5_ves = Diffusion.Create(Rab3_TOMO_synapsin_CDK5, DCST_MEM_VES)
        diff_Rab3_TOMO_synapsin_p_CDK5_ves = Diffusion.Create(Rab3_TOMO_synapsin_p_CDK5, DCST_MEM_VES)
        diff_Rab3_TOMO_p_synapsin_CaN_CaM_ves = Diffusion.Create(Rab3_TOMO_p_synapsin_CaN_CaM, DCST_MEM_VES)
        diff_Rab3_TOMO_p_synapsin_p_CaN_CaM_ves = Diffusion.Create(Rab3_TOMO_p_synapsin_p_CaN_CaM, DCST_MEM_VES)

        CDK5.o + Rab3_TOMO.v <r[1]> Rab3_TOMO_CDK5.v
        r[1].K = kon_cdk_tomo, koff_cdk_tomo
        Rab3_TOMO_CDK5.v >r[1]> CDK5.o + Rab3_TOMO_p.v
        r[1].K = kcat_cdk_tomo

        CaN_CaM.o + Rab3_TOMO_p.v <r[1]> Rab3_TOMO_p_CaN_CaM.v
        r[1].K = kon_can_tomo, koff_can_tomo
        Rab3_TOMO_p_CaN_CaM.v >r[1]> CaN_CaM.o + Rab3_TOMO.v
        r[1].K = kcat_can_tomo
        TOMO_p.o + Rab3.v <r[1]> Rab3_TOMO_p.v
        r[1].K = kon_rab3_TOMO_p, koff_rab3_TOMO_p
        Rab3_TOMO.v >r[1]> TOMO.o + Rab3.v
        r[1].K = koff_rab3_TOMO_p
        Rab3_TOMO_p.v + synapsin.v <r[1]> Rab3_TOMO_p_synapsin.v
        r[1].K = kon_TOMO_p_synapsin, koff_TOMO_p_synapsin
        Rab3_TOMO_p.v + synapsin_p.v <r[1]> Rab3_TOMO_p_synapsin_p.v
        r[1].K = kon_TOMO_p_synapsin, koff_TOMO_p_synapsin
        Rab3_TOMO.v + synapsin.v <r[1]> Rab3_TOMO_synapsin.v
        r[1].K = kon_TOMO_synapsin, koff_TOMO_synapsin
        Rab3_TOMO.v + synapsin_p.v <r[1]> Rab3_TOMO_synapsin_p.v
        r[1].K = kon_TOMO_synapsin, koff_TOMO_synapsin

        #free vesicle synapsin phosphorylation
        PKA.o + synapsin.v <r[1]> PKA_synapsin.v
        r[1].K = kon_PKA_syn, koff_PKA_syn
        PKA_synapsin.v >r[1]> PKA.o + synapsin_p.v
        r[1].K = kcat_PKA_syn
        
        PKA.o + Rab3_TOMO_synapsin.v <r[1]> Rab3_TOMO_synapsin_PKA.v
        r[1].K = kon_PKA_syn, koff_PKA_syn
        Rab3_TOMO_synapsin_PKA.v >r[1]> PKA.o + Rab3_TOMO_synapsin_p.v
        r[1].K = kcat_PKA_syn
        
        PKA.o + Rab3_TOMO_p_synapsin.v <r[1]> Rab3_TOMO_p_synapsin_PKA.v
        r[1].K = kon_PKA_syn, koff_PKA_syn
        Rab3_TOMO_p_synapsin_PKA.v >r[1]> PKA.o + Rab3_TOMO_p_synapsin_p.v
        r[1].K = kcat_PKA_syn

        #dephos of vesicle synapsin_p
        PP2A.o + synapsin_p.v <r[1]> PP2A_synapsin_p.v
        r[1].K = kon_pp2a_synapsin, koff_pp2a_synapsin
        PP2A_synapsin_p.v >r[1]> PP2A.o + synapsin.v
        r[1].K = kcat_pp2a_syn
        
        PP2A.o + Rab3_TOMO_p_synapsin_p.v <r[1]> Rab3_TOMO_p_synapsin_p_PP2A.v
        r[1].K = kon_pp2a_synapsin, koff_pp2a_synapsin
        Rab3_TOMO_p_synapsin_p_PP2A.v >r[1]> PP2A.o + Rab3_TOMO_p_synapsin.v
        r[1].K = kcat_pp2a_syn
        
        PP2A.o + Rab3_TOMO_synapsin_p.v <r[1]> Rab3_TOMO_synapsin_p_PP2A.v
        r[1].K = kon_pp2a_synapsin, koff_pp2a_synapsin
        Rab3_TOMO_synapsin_p_PP2A.v >r[1]> PP2A.o + Rab3_TOMO_synapsin.v
        r[1].K = kcat_pp2a_syn

        CDK5.o + Rab3_TOMO_synapsin.v <r[1]> Rab3_TOMO_synapsin_CDK5.v
        r[1].K = kon_cdk_tomo, koff_cdk_tomo
        Rab3_TOMO_synapsin_CDK5.v >r[1]> CDK5.o + Rab3_TOMO_p_synapsin.v
        r[1].K = kcat_cdk_tomo
        
        CDK5.o + Rab3_TOMO_synapsin_p.v <r[1]> Rab3_TOMO_synapsin_p_CDK5.v
        r[1].K = kon_cdk_tomo, koff_cdk_tomo
        Rab3_TOMO_synapsin_p_CDK5.v >r[1]> CDK5.o + Rab3_TOMO_p_synapsin_p.v
        r[1].K = kcat_cdk_tomo

        CaN_CaM.o + Rab3_TOMO_p_synapsin.v <r[1]> Rab3_TOMO_p_synapsin_CaN_CaM.v
        r[1].K = kon_can_tomo, koff_can_tomo
        Rab3_TOMO_p_synapsin_CaN_CaM.v >r[1]> CaN_CaM.o + Rab3_TOMO_synapsin.v
        r[1].K = kcat_can_tomo
        
        CaN_CaM.o + Rab3_TOMO_p_synapsin_p.v <r[1]> Rab3_TOMO_p_synapsin_p_CaN_CaM.v
        r[1].K = kon_can_tomo, koff_can_tomo
        Rab3_TOMO_p_synapsin_p_CaN_CaM.v >r[1]> CaN_CaM.o + Rab3_TOMO_synapsin_p.v
        r[1].K = kcat_can_tomo

        #####Rab3-TOMO binds to the synapsin dimers#####
        Rab3_TOMO_p.v + synapsin_dimer.v <r[1]> Rab3_TOMO_p_synapsin_dimer.v
        r[1].K = kon_TOMO_p_synapsin, koff_TOMO_p_synapsin
        Rab3_TOMO_p.v + synapsin_dimer_p.v <r[1]> Rab3_TOMO_p_synapsin_dimer_p.v
        r[1].K = kon_TOMO_p_synapsin, koff_TOMO_p_synapsin
        Rab3_TOMO.v + synapsin_dimer.v <r[1]> Rab3_TOMO_synapsin_dimer.v
        r[1].K = kon_TOMO_synapsin, koff_TOMO_synapsin
        Rab3_TOMO.v + synapsin_dimer_p.v <r[1]> Rab3_TOMO_synapsin_dimer_p.v
        r[1].K = kon_TOMO_synapsin, koff_TOMO_synapsin

        #PKA phos.
        PKA.o + synapsin_dimer.v <r[1]> PKA_synapsin_dimer.v >r[2]> PKA.o + synapsin_dimer_p.v
        r[1].K = kon_PKA_syn, koff_PKA_syn
        r[2].K = kcat_PKA_syn
        
        PKA.o + Rab3_TOMO_synapsin_dimer.v <r[1]> Rab3_TOMO_synapsin_dimer_PKA.v >r[2]> PKA.o + Rab3_TOMO_synapsin_dimer_p.v
        r[1].K = kon_PKA_syn, koff_PKA_syn
        r[2].K = kcat_PKA_syn
        
        PKA.o + Rab3_TOMO_p_synapsin_dimer.v <r[1]> Rab3_TOMO_p_synapsin_dimer_PKA.v >r[2]> PKA.o + Rab3_TOMO_p_synapsin_dimer_p.v
        r[1].K = kon_PKA_syn, koff_PKA_syn
        r[2].K = kcat_PKA_syn
        
        PP2A.o + synapsin_dimer_p.v <r[1]> PP2A_synapsin_dimer_p.v >r[2]> PP2A.o + synapsin_dimer.v
        r[1].K = kon_pp2a_synapsin, koff_pp2a_synapsin
        r[2].K = kcat_pp2a_syn
        
        PP2A.o + Rab3_TOMO_synapsin_dimer_p.v <r[1]> Rab3_TOMO_synapsin_dimer_p_PP2A.v >r[2]> PP2A.o + Rab3_TOMO_synapsin_dimer.v
        r[1].K = kon_pp2a_synapsin, koff_pp2a_synapsin
        r[2].K = kcat_pp2a_syn
        
        PP2A.o + Rab3_TOMO_p_synapsin_dimer_p.v <r[1]> Rab3_TOMO_p_synapsin_dimer_p_PP2A.v >r[2]> PP2A.o + Rab3_TOMO_p_synapsin_dimer.v
        r[1].K = kon_pp2a_synapsin, koff_pp2a_synapsin
        r[2].K = kcat_pp2a_syn
        
        CDK5.o + Rab3_TOMO_synapsin_dimer.v <r[1]> Rab3_TOMO_synapsin_dimer_CDK5.v >r[2]> CDK5.o + Rab3_TOMO_p_synapsin_dimer.v
        r[1].K = kon_cdk_tomo, koff_cdk_tomo
        r[2].K = kcat_cdk_tomo
        
        CDK5.o + Rab3_TOMO_synapsin_dimer_p.v <r[1]> Rab3_TOMO_synapsin_dimer_p_CDK5.v >r[2]> CDK5.o + Rab3_TOMO_p_synapsin_dimer_p.v
        r[1].K = kon_cdk_tomo, koff_cdk_tomo
        r[2].K = kcat_cdk_tomo

        CaN_CaM.o + Rab3_TOMO_p_synapsin_dimer.v <r[1]> Rab3_TOMO_p_synapsin_dimer_CaN_CaM.v >r[2]> CaN_CaM.o + Rab3_TOMO_synapsin_dimer.v
        r[1].K = kon_can_tomo, koff_can_tomo
        r[2].K = kcat_can_tomo
        
        CaN_CaM.o + Rab3_TOMO_p_synapsin_dimer_p.v <r[1]> Rab3_TOMO_p_synapsin_dimer_p_CaN_CaM.v >r[2]> CaN_CaM.o + Rab3_TOMO_synapsin_dimer_p.v
        r[1].K = kon_can_tomo, koff_can_tomo
        r[2].K = kcat_can_tomo

        #binding of TOMO to SNARE (independent of CDK5 phos state, indicated by small p)
        #######TOMO is incorporated into the SNARE complex (binds to SYX_SNP25 complex)
        TOMO.o + RIM_M13_Rab3_SYX_M18_SNP25.v <r[1]> RIM_M13_Rab3_SYX_M18_SNP25_TOMO.v
        r[1].K = kon_tomo_syx, koff_tomo_syx
        
        TOMO_p.o + RIM_M13_Rab3_SYX_M18_SNP25.v <r[1]> RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p.v
        r[1].K = kon_tomo_syx, koff_tomo_syx
        
        #switching of TOMO to SYB-displaceable state (unphosphorylated, slow)
        RIM_M13_Rab3_SYX_M18_SNP25_TOMO.v <r[1]> RIM_M13_Rab3_SYX_M18_SNP25_TOMOx.v
        r[1].K = kflip_tomo_f, kflip_tomo_r
        
        RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p.v <r[1]> RIM_M13_Rab3_SYX_M18_SNP25_TOMOx_p.v
        r[1].K = kflip_tomo_f, kflip_tomo_r
        
        #switching of TOMO to SYB-displaceable state (phosphorylated by PKA indicated by big P, fast)
        RIM_M13_Rab3_SYX_M18_SNP25_TOMO_P.v <r[1]> RIM_M13_Rab3_SYX_M18_SNP25_TOMOx_P.v
        r[1].K = kflip_tomo_P_f, kflip_tomo_r
        
        RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p.v >r[1]> RIM_M13_Rab3_SYX_M18_SNP25_TOMOx_p.v
        r[1].K = kflip_tomo_P_f
        RIM_M13_Rab3_SYX_M18_SNP25_TOMOx_p_P.v >r[1]> RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p_P.v
        r[1].K = kflip_tomo_r
        
        #displacement of TOMO by SYB
        RIM_M13_Rab3_SYX_M18_SNP25_TOMOx.v + SYB.v >r[1]> SNARE.v + TOMO.v
        r[1].K = kdisp_tomo_syx
        RIM_M13_Rab3_SYX_M18_SNP25_TOMOx_p.v + SYB.v >r[1]> SNARE.v + TOMO_p.v
        r[1].K = kdisp_tomo_syx
        
        #displacement of TOMO_P by SYB
        RIM_M13_Rab3_SYX_M18_SNP25_TOMOx_P.v + SYB.v >r[1]> SNARE.v + TOMO.v
        r[1].K = kdisp_tomo_syx
        RIM_M13_Rab3_SYX_M18_SNP25_TOMOx_p_P.v + SYB.v >r[1]> SNARE.v + TOMO_p.v
        r[1].K = kdisp_tomo_syx
        
        #phosphorylation of TOMO by PKA
        PKA.o + RIM_M13_Rab3_SYX_M18_SNP25_TOMO.v <r[1]> RIM_M13_Rab3_SYX_M18_SNP25_TOMO_PKA.v
        r[1].K = kon_tomo_pka, koff_tomo_pka
        RIM_M13_Rab3_SYX_M18_SNP25_TOMO_PKA.v >r[1]> PKA.o + RIM_M13_Rab3_SYX_M18_SNP25_TOMO_P.v
        r[1].K = kcat_pka_tomo
        
        PKA.o + RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p.v <r[1]> RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p_PKA.v
        r[1].K = kon_tomo_pka, koff_tomo_pka
        RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p_PKA.v >r[1]> PKA.o + RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p_P.v
        r[1].K = kcat_pka_tomo
        
        #dephosphorylation of TOMO by PP2A
        PP2A.o + RIM_M13_Rab3_SYX_M18_SNP25_TOMO_P.v <r[1]> RIM_M13_Rab3_SYX_M18_SNP25_TOMO_P_PP2A.v
        r[1].K = kon_tomo_PP2A, koff_tomo_PP2A
        RIM_M13_Rab3_SYX_M18_SNP25_TOMO_P_PP2A.v >r[1]> PP2A.o + RIM_M13_Rab3_SYX_M18_SNP25_TOMO.v
        r[1].K = kcat_PP2A_tomo
        
        PP2A.o + RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p_P.v <r[1]> RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p_P_PP2A.v
        r[1].K = kon_tomo_PP2A, koff_tomo_PP2A
        RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p_P_PP2A.v >r[1]> PP2A.o + RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p.v
        r[1].K = kcat_PP2A_tomo
        
        #####BINDING OF PROTEINS TO VESICLES#######
        ######Binding of synapsin to vesicles
        kon_synapsin_ves = 10.7e06
        koff_synapsin_ves = 0.123
        koff_synapsin_p_ves = 123
        synapsin.o + ves_syn1_site.v <r[1]> synapsin.v
        r[1].K = kon_synapsin_ves, koff_synapsin_ves
        synapsin_p.v >r[1]> synapsin_p.o + ves_syn1_site.v
        r[1].K = koff_synapsin_p_ves

        CXN.o + ves_CXN_site.v <r[1]> CXN.v
        r[1].K = kon_cxn_ves, koff_cxn_ves
        Rab3.o + ves_rab3_site.v <r[1]> Rab3.v
        r[1].K = kon_rab3_ves, koff_rab3_ves
        M13.o + ves_m13_site.v <r[1]> M13.v
        r[1].K = kon_m13_ves, koff_m13_ves
        M18.o + ves_m18_site.v <r[1]> M18.v
        r[1].K = kon_m18_ves, koff_m18_ves
        SYN1.o + ves_syndap_site.v <r[1]> SYN1.v
        r[1].K = kon_syn1_ves, koff_syn1_ves
        aSNAP.o + ves_asnap_site.v <r[1]> aSNAP.v
        r[1].K = kon_asnap_ves, koff_asnap_ves
        NSF.o + ves_nsf_site.v <r[1]> NSF.v
        r[1].K = kon_nsf_ves, koff_nsf_ves

    with raft_sys:
        #raft dephosphorylation
        CaN_CaM.i + DYNpp.r <r[1]> CaN_DYNpp.r
        r[1].K = kon_can_dyn, koff_can_dyn
        CaN_DYNpp.r >r[1]> CaN_CaM.i + DYNp.r
        r[1].K = kcat_can
        
        CaN_CaM.i + DYNp.r <r[1]> CaN_DYNp.r
        r[1].K = kon_can_dyn, koff_can_dyn
        CaN_DYNp.r >r[1]> CaN_CaM.i + DYN.r
        r[1].K = kcat_can

        #Raft binding of dynamin to syndapin 1 (SYN1)
        SYN1.i + DYN.r <r[1]> DYN_SYN1.r
        r[1].K = kon_dyn_syn, koff_dyn_syn

        SYB.s + AP180.r >r[1]> SYB.r
        r[1].K = kon_syb_ap180
        syt.s + AP2.r >r[1]> syt.r
        r[1].K = kon_syt_ap2
        Rab3.s + RabAd.r >r[1]> Rab3.r
        r[1].K = kon_rab3_rabad

        DYN_SYN1.s + DYN_AD.r <r[1]> DYN_SYN1.r
        r[1].K = kon_dyn_raft, koff_dyn_raft
        DYNp.s + DYN_AD.r <r[1]> DYNp.r
        r[1].K = kon_dyn_raft, koff_dyn_raft
        DYN.s + DYN_AD.r <r[1]> DYN.r
        r[1].K = kon_dyn_raft, koff_dyn_raft

        ##formation of new vesicle by endocytosis
        raftendo = RaftEndocytosis.Create(ves, 1000.0, deps=DYN_SYN1.r + 10 * Rab3.r + 26 * syt.r + 68 * SYB.r, inside=True)
        
####################################
#FIND ALL SPECIES WITHOUT MEMBRANE DIFFUSION RULE

vsysSpecs = set(spec.getID() for spec in vsys._getStepsObjects()[0].getAllSpecs())
memsysSpecs = set(spec.getID() for spec in memsys._getStepsObjects()[0].getAllSpecs())
ves_ssysSpecs = set(spec.getID() for spec in ves_ssys._getStepsObjects()[0].getAllSpecs())
raft_sysSpecs = set(spec.getID() for spec in raft_sys._getStepsObjects()[0].getAllSpecs())
print(f'TOTAL CYTOSOL SPECIES: {len(vsysSpecs)}')
print(f'TOTAL MEMBRANE SPECIES: {len(memsysSpecs)}')
print(f'TOTAL VES SPECIES: {len(ves_ssysSpecs)}')

#find intersect of vesicle species and surface species
ves_surf_spec_diff = ves_ssysSpecs - memsysSpecs
print(f'There are {len(ves_surf_spec_diff)} ves species not in the surface system')
print('MISSING SURF SPECIES:', ves_surf_spec_diff)

print(f'TOTAL RAFT SPECIES: {len(raft_sysSpecs)}')
raft_surf_spec_diff = raft_sysSpecs - memsysSpecs
print(f'There are {len(raft_surf_spec_diff)} raft species not in the surface system')
print('MISSING RAFT SPECIES:', raft_surf_spec_diff)


####Geometry Definition####
scale = 1e-9
scale_100 = 1e-7

mesh = TetMesh.Load(MESH_PATH)

print("MESH LOADED")

with mesh:
    #find min and max y-values of outer_tets bounding box
    cyt_y_min = min(tet.center.y for tet in mesh.outer.tets)
    cyt_y_max = max(tet.center.y for tet in mesh.outer.tets)

    print("cyt_y_min", cyt_y_min)
    print("cyt_y_max", cyt_y_max)

    print('FOUND OUTER TETS')

    cyt_tets = TetList(tet for tet in mesh.spine.tets if cyt_y_min < tet.center.y < cyt_y_max)
    axon_tets = mesh.spine.tets - cyt_tets

    print('ALL tets in mesh:', len(mesh.tets))
    print('Total cytosol tets:', len(mesh.spine.tets))
    print('Total axonal tets:', len(axon_tets))
    print('Total bouton tets:', len(cyt_tets))
    print('Sum bouton + axomal:', len(axon_tets)+len(cyt_tets))
    print('Sum bouton, axon, ER, outer:', len(mesh.spine.tets)+len(mesh.er.tets)+len(mesh.outer.tets))

    ERmemb_tris = cyt_tets.surface & mesh.er.tets.surface
    memb_tris = cyt_tets.surface & mesh.outer.tets.surface
    print('Number of triangles on membrane =', len(memb_tris))

    #open actin_tets list
    print('OPENING ACTIN TETS...')
    with open(os.path.join(MODEL_PATH, 'actin_tets_215nm.txt'), 'rb') as f:
        actin_tets = pickle.load(f)
    print('ACTIN TETS LOADED:', len(actin_tets))

    ######FIND AZONE TRIANGLES######
    #find membrane triangles within ellipse enclosing coordinates
    #ellipse foci
    focus1 = np.array([2036.7e-9, 2982.6e-9])
    focus2 = np.array([2185.8e-9, 3280.8e-9])
    #find memb triangles closest to foci
    focus1_tri = min(memb_tris, key=lambda tri: np.linalg.norm(tri.center[:2] - focus1))
    focus2_tri = min(memb_tris, key=lambda tri: np.linalg.norm(tri.center[:2] - focus2))
    print("Ellipse focus 1 tri:", focus1_tri)
    print("Ellipse focus 2 tri:", focus2_tri)

    radius_maj = 380e-9#349.8e-9
    azone_tris = TriList(tri for tri in memb_tris if np.linalg.norm(tri.center - focus1_tri.center) + np.linalg.norm(tri.center - focus2_tri.center) < radius_maj)
    print('Number of triangles in active zone =', len(azone_tris))
    print("active zone area =", azone_tris.Area)

    ######FIND ENDOCYTIC ZONE TRIANGLES######
    radius_maj_endo = 550e-9 #550e-9#349.8e-9
    endo_tris = TriList(tri for tri in memb_tris if radius_maj <= np.linalg.norm(tri.center - focus1_tri.center) + np.linalg.norm(tri.center - focus2_tri.center) < radius_maj_endo)
    print('Number of triangles in endo zone:', len(endo_tris))
    print("endo zone area =", endo_tris.Area)
    print("membrane area =", memb_tris.Area)
    print("ER vol =", mesh.er.tets.Vol)

    #create compartments
    #cytosol includes all tets not in the recycling vesicle AZ cluster
    axon = Compartment.Create(axon_tets)
    cytosol = Compartment.Create(cyt_tets, vsys)
    cytER = Compartment.Create(mesh.er.tets, cytERsys)
    cleft = Compartment.Create(mesh.outer.tets, syn_sys)

    vac1 = Compartment.Create(mesh.vec1.tets)
    vac2 = Compartment.Create(mesh.vec2.tets)
    vac3 = Compartment.Create(mesh.vec3.tets)
    vac4 = Compartment.Create(mesh.vec4.tets)
    vac5 = Compartment.Create(mesh.vec5.tets)
    mito = Compartment.Create(mesh.mit.tets)

    print("cytosol vol =", cytosol.Vol)
    print("axon vol =", axon.Vol)

    ERmemb = Patch.Create(ERmemb_tris, cytER, cytosol, ERsys)
    memb = Patch.Create(memb_tris, cytosol, cleft, memsys)

    cyt_tris = mesh.surface & cytosol.surface
    cyt_verts = VertList(set(vert for tet in cyt_tets for vert in tet.verts))
    print("Number of cyt_verts =", len(cyt_verts))

    # Create excitable membrane
    membrane = Membrane.Create([memb])
    print("MEMBRANE CREATED")

    #find tetrahedra in restricted diffusion "tethering" zone around Azone.
    tether_zone = TetList(tet for tet in cyt_tets if any(np.linalg.norm(tri.center - tet.center) < 100e-09 for tri in azone_tris))
    print("Tetrahedra in tether zone:", len(tether_zone))

    dock_tris = TriList(dock_tris)
    #create dictionary of dock tris and surrounding triangles (to capture escaped RIM)
    all_surr_tris = TriList()
    azone_and_endo_tris = azone_tris + endo_tris
    for tri in dock_tris:
        surr = TriList(t for t in azone_and_endo_tris if np.linalg.norm(t.center - tri.center) < 45e-9)
        all_surr_tris += surr
        print(f'Dock tri {tri} has {len(surr)} surrounding tris.')
        print(f'Dock tri {tri} has {surr.indices} surrounding tris.')
    all_surr_tris_set = set(all_surr_tris)
    print('TRI SURROUNDING HAD DUPLICATES:', len(all_surr_tris) != len(all_surr_tris_set))
    print(len(all_surr_tris), 'vs', len(all_surr_tris_set))
    print('DUPLICATES:', [t for t, v in collections.Counter(all_surr_tris).items() if v > 1])

    #add each surrounding triangle to closest dock triangle
    dock_tri_to_surround_tris = {}
    for tri in all_surr_tris_set:
        if tri not in dock_tris:
            dock_tri = min(dock_tris, key=lambda dt: np.linalg.norm(dt.center - tri.center))
            dock_tri_to_surround_tris.setdefault(dock_tri, TriList()).append(tri)

    surr_tri_count = 0
    for d, surr in dock_tri_to_surround_tris.items():
        print(f'Dock tris {d} surrond tris: {surr.indices}')
        surr_tri_count += len(surr)
    print('TOTAL SURR TRIS:', surr_tri_count)

    #create dictionary of dock tris and "tether path" ends (tri:[[x0,y0],[x1,y1]])
    dock_tri_to_tether_path = {}
    #find norms for all dock_tris
    print('FINDING TETHER PATH POINTS')
    for tri in dock_tris:
        print('DOCK TRI COORD:', tri.center)
        outside = tri.center - tri.norm * 100e-9
        inside = tri.center + tri.norm * 100e-9
        # Make sure it points in cytosol
        if mesh.tets[inside] not in cyt_tets:
            print('NORM NOT IN CYTOSOL.... INVERTING')
            outside, inside = inside, outside
        dock_tri_to_tether_path[tri] = (outside, inside)
    print(f'FOUND ALL {dock_tri_to_tether_path} TETHER PATH POINTS')

    min_dist_dock = min(np.linalg.norm(tri1.center - tri2.center) for i, tri1 in enumerate(dock_tris) for tri2 in dock_tris[i+1:])
    print(f'DOCK MIN DIST: {min_dist_dock * 1e09} nm')

    print("active zone area = ", azone_tris.Area)
    print("endo zone area = ", endo_tris.Area)
    mem_area = cyt_tris.Area
    print("membrane area = ", mem_area)
    print("ER vol =", cytER.Vol)


#################################### SIMULATION DATA FUNCTIONS ETC ####################################


def getTOMOcount(rs):
    TOMO_vesSpecs = [
        'Rab3_TOMO', 'Rab3_TOMO_synapsin', 'Rab3_TOMO_synapsin_actin', 'Rab3_TOMO_synapsin_PKA', 'Rab3_TOMO_synapsin_p',
        'Rab3_TOMO_p_CaN_CaM', 'Rab3_TOMO_CDK5', 'Rab3_TOMO_synapsin_CDK5', 'Rab3_TOMO_synapsin_p_PP2A',
        'Rab3_TOMO_synapsin_dimer_p_PP2A', 'Rab3_TOMO_synapsin_dimer_PKA', 'Rab3_TOMO_synapsin_dimer',
        'Rab3_TOMO_synapsin_dimer_p', 'Rab3_TOMO_synapsin_dimer_CDK5', 'Rab3_TOMO_synapsin_dimer_p_CDK5',
    ]
    cyt_TOMO_count = rs.SUM(rs.cytosol.LIST('TOMO', 'TOMO_CDK5').Count)
    ves_TOMO_count = rs.SUM(rs.cytosol.ves('surf').LIST(*TOMO_vesSpecs).Count)
    return cyt_TOMO_count + ves_TOMO_count

def getVesicleLinkSpecs(vesref):
    linkspecs = [
        'synapsin_dimer', 'synapsin_dimer_p', 'Rab3_TOMO_synapsin_dimer', 'Rab3_TOMO_synapsin_dimer_p',
        'Rab3_TOMO_p_synapsin_dimer', 'Rab3_TOMO_p_synapsin_dimer_p', 'PKA_synapsin_dimer', 'Rab3_TOMO_synapsin_dimer_PKA',
        'Rab3_TOMO_p_synapsin_dimer_PKA', 'PP2A_synapsin_dimer_p', 'Rab3_TOMO_synapsin_dimer_p_PP2A',
        'Rab3_TOMO_p_synapsin_dimer_p_PP2A', 'Rab3_TOMO_synapsin_dimer_CDK5', 'Rab3_TOMO_synapsin_dimer_p_CDK5',
        'Rab3_TOMO_p_synapsin_dimer_CaN_CaM', 'Rab3_TOMO_p_synapsin_dimer_p_CaN_CaM',
    ]
    return sum(vesref('surf').LIST(*linkspecs).Count)

def getPhosphoSynapsins(rs):
    syn_phos_ves_specs = [
        'synapsin_p',
        'Rab3_TOMO_synapsin_actin_p',
        'Rab3_TOMO_p_synapsin_actin_p',
        'Rab3_TOMO_p_synapsin_p',
        'Rab3_TOMO_synapsin_p',
        'Rab3_TOMO_p_synapsin_p_PP2A',
        'Rab3_TOMO_synapsin_p_PP2A',
        'Rab3_TOMO_synapsin_p_CDK5',
        'Rab3_TOMO_p_synapsin_p_CaN_CaM',
        'synapsin_dimer_p',
        'Rab3_TOMO_synapsin_dimer_p',
        'Rab3_TOMO_p_synapsin_dimer_p',
        'PP2A_synapsin_dimer_p',
        'Rab3_TOMO_synapsin_dimer_p_PP2A',
        'Rab3_TOMO_p_synapsin_dimer_p_PP2A',
        'Rab3_TOMO_synapsin_dimer_p_CDK5',
        'Rab3_TOMO_p_synapsin_dimer_p_CaN_CaM',
    ]
    cyt_syn_phos = rs.SUM(rs.cytosol.LIST(synapsin_p, synapsin_actin_p).Count)
    ves_syn_phos = rs.SUM(rs.cytosol.ves('surf').LIST(*syn_phos_ves_specs).Count)
    return cyt_syn_phos + ves_syn_phos

rng = RNG('mt19937', 512, rng_seed)

print ("Creating MPI solver...")
use_efield = True
sim = Simulation('TetVesicle', syt_model, mesh, rng, MPI.EF_DV_PETSC if use_efield else MPI.EF_NONE)
print ("MPI solver created.")

def generateDockVesicle(sim, tri, dock_tri_to_tether_path, dock_positions):
    vesref = sim.cytosol.addVesicle('ves')
    outside, inside = dock_tri_to_tether_path[tri]
    pos = tri.center + (inside - tri.center) / np.linalg.norm(inside - tri.center) * (ves.Diameter / 2 + 1e-9)
    vesref.Pos = pos
    # test vesicle actually moved
    if np.linalg.norm(vesref.Pos - pos) < 1e-9:
        dock_positions.append(pos)
        print("Added dock vesicle")
    else:
        print("VESICLE DOCK ERROR!")
    return vesref

FUSED_VES = []
fusion_count = 0

# Data recording
rs = ResultSelector(sim)

custom_values = CustomResults(sim, [int] * 9 + [list] * 4)
custom_values.labels = [
    'len(nonrec_ves_clust)',
    'len(nonrec_ves_free)',
    'len(docked_ves)',
    'len(init_vesicles_used)',
    'len(init_rec_vesicles_used)',
    'len(init_nonrec_vesicles_used)',
    'len(new_vesicles_used)',
    'len(init_rec_vesicles_used) + len(init_nonrec_vesicles_used) + len(new_vesicles_used)',
    'total_link_specs',
    'VES_RAB3_ALL',
    'VES_RAB3_REC',
    'VES_RAB3_RES',
    'rel_probs',
]

recorded_values = (
       rs.cytosol.LIST(CaM[N2, C2], AC18_CaM, cAMP, PKA, Ca).Count
    << rs.cleft.REL_IND.Count \
    << rs.SUM(rs.cytosol.LIST(CaN_CaM, CaN_DYNpp, CaN_DYNp).Count) # count active Calcineurin all forms in cytosol
    << getTOMOcount(rs) # count all dephosphorylated TOMO
    << rs.SUM(rs.cytosol.LIST(TOMO, TOMO_p).Count) # count all free cytosolic TOMO
    << rs.SUM(rs.cytosol.ves('surf').LIST(
        'RIM_M13_Rab3_SYX_M18_SNP25_TOMO',
        'RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p',
        'RIM_M13_Rab3_SYX_M18_SNP25_TOMO_P',
        'RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p_P',
        'RIM_M13_Rab3_SYX_M18_SNP25_TOMO_PKA',
        'RIM_M13_Rab3_SYX_M18_SNP25_TOMO_p_PKA',
    ).Count) # Count TOMO bound to Syx in non-displaceable state
    << rs.SUM(rs.cytosol.ves('surf').LIST(
                    'RIM_M13_Rab3_SYX_M18_SNP25_TOMOx',
                    'RIM_M13_Rab3_SYX_M18_SNP25_TOMOx_p',
                    'RIM_M13_Rab3_SYX_M18_SNP25_TOMOx_P',
                    'RIM_M13_Rab3_SYX_M18_SNP25_TOMOx_p_P',
    ).Count) # Count TOMO bound to Syx in displaceable state
    << rs.raftendo.Extent # count recycled vesicles
    << rs.SUM(rs.memb.raft.LIST(DYNpp, CaN_DYNpp).Count)
    << rs.SUM(rs.memb.raft.LIST(DYNp, CaN_DYNp).Count)
    << rs.memb.raft.DYN.Count
    << rs.SUM(rs.memb.raft.LIST(DYN_SYN1, DYN_SYN1_CDK5).Count) # count dynamin in rafts
    << rs.memb.raft.SNARE_DISS.Count # count disassembled SNAREs in Rafts
    << rs.memb.SYB.Count #count SYB in membrane
    << rs.memb.raft.SYB.Count #count SYB in Rafts
    << rs.memb.raft.Count #count Rafts
    << rs.cytosol.ves('surf').Rab3.Count #count free Rab3 on vesicles
    << getPhosphoSynapsins(rs) #count phosphorylated synapsin
)

if save_blender_data:
    blender_selectors = [
        rs.TETS(cytosol.tets).Ca.Count,
        rs.cytosol.VESICLES().Pos,
        rs.cytosol.VESICLES()('surf').POINTSPECS([Rab3, SYB, syt]).PosSpherical,
        rs.cytosol.VESICLES()('in').glu.Count,
        rs.cytosol.VESICLES()('surf').LINKSPECS().Pos,
        rs.cytosol.VESICLES()('surf').LINKSPECS().LinkedTo,
        rs.memb.RAFTS().Pos,
        rs.memb.RAFTS().LIST(Rab3, SYB, syt, DYNp).Count,
        rs.ALL(Exocytosis, RaftEndocytosis).Events,
    ]
    sim.toSave(recorded_values, *blender_selectors, dt=DT)
else:
    sim.toSave(recorded_values, dt=DT)

sim.toSave(custom_values)

os.makedirs(DATA_PATH, exist_ok=True)
hdfPrefix = os.path.join(DATA_PATH, f'{DATA_ID}_{JOB_ID}_{JOB_INDEX}_results')
with XDMFHandler(hdfPrefix) as hdf:
    sim.toDB(hdf, f'{DATA_ID}_{JOB_ID}_{JOB_INDEX}')

    for j in range(NITER):
        sim.newRun()

        add_tethers_time = 0.7
        pka_act_time = 50
        camkii_deact_time = 50
        rel_prob_test_time = 0.05
        print("Running iteration", j + 1)

        #create list of pulse times
        pulse_times = np.round(np.arange(5000, INT*1000, 1000 / RATE))
        pulse_times_off = pulse_times+1
        check_release_times_pre = pulse_times-1
        check_release_times = pulse_times+5

        if checkpoint_mode in ['restore', 'both']:
            # First restore the simulation to get the previous time
            sim.restore(checkpoint_path)

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

            # If the events already happened, prevent them from happening again
            if sim.Time > add_tethers_time:
                add_tethers_time = math.inf
            if sim.Time > pka_act_time:
                pka_act_time = math.inf
            if sim.Time > camkii_deact_time:
                camkii_deact_time = math.inf
            if sim.Time > rel_prob_test_time:
                rel_prob_test_time = math.inf

            # Restore the state of script variables
            pkl_path = checkpoint_path + '_script_state.pkl'
            with open(pkl_path, 'rb') as f:
                def restore_vesraft_list(f, cls=VesicleReference):
                    return [cls(sim, vesTpe, vesIdx) for vesTpe, vesIdx in pickle.load(f)]

                dock_positions = pickle.load(f)
                dock_ves_init_tets = {
                    VesicleReference(sim, vesTpe, vesIdx): mesh.tets[tetIdx]
                    for (vesTpe, vesIdx), tetIdx in pickle.load(f).items()
                }
                (nonrec_ves, nonrec_ves_init, rec_ves, recycled_ves_free, recycled_ves_all, init_vesicle_ind, init_vesicles_used, 
                    init_rec_vesicles_used, init_nonrec_vesicles_used, new_vesicles_used, fused_ind_all, REC_FUSIONS, RES_FUSIONS,
                    rec_ves_init, all_ves_pre_ap) = [
                        restore_vesraft_list(f) for i in range(15)
                    ]
                new_rafts = restore_vesraft_list(f, cls=RaftReference)
                filled_rafts = restore_vesraft_list(f, cls=RaftReference)
                r2v_r = restore_vesraft_list(f, cls=RaftReference)
                r2v_v = restore_vesraft_list(f)
                raft_to_vesicle = {r:v for r, v in zip(r2v_r, r2v_v)}
                (fusion_events_tot, rel_ind_count_pre, rel_ind_count_post, total_raft_endo,
                    total_time, action_potentials_fired, ap_n, rel_probs, recycled_vesicles_tot,
                    raft_count_init, synapsin_total, rel_ind_count_pre_ap) = pickle.load(f)
                rel_probs = np.append(rel_probs, np.zeros(len(pulse_times)))
        else:
            # First time setup (no restore from checkpoint)

            if use_efield:
                print("E-FIELD IS ACTIVE")
                sim.EfieldDT = EF_DT
                sim.membrane.Potential = init_pot
                sim.membrane.VolRes = Ra
                sim.membrane.Capac = memb_capac
                sim.VERTS(cyt_verts).VClamped = True
            else:
                print("E-FIELD IS INACTIVE")

            # #################ADD DOCK VESICLES#####################
            dock_ves_init = []
            dock_positions = []
            for tri in dock_tris:
                vesref = generateDockVesicle(sim, tri, dock_tri_to_tether_path, dock_positions)
                dock_ves_init.append(vesref)
                print('Added ves:', vesref.idx)
                vesref('surf').Rab3.Count = 10
                vesref('surf').SYB.Count = 69
                vesref('surf').syt.Count = 27
                vesref('in').REL_IND.Count = 1

            print('Dock vesicles added:', len(dock_positions))
            print('Dock triangles:', dock_tris.indices)

            #find tetrahedron for each dock vesicle
            dock_ves_init_tets = {vesref: mesh.tets[vesref.Pos] for vesref in dock_ves_init}
            #set initial diffusion rate for dock vesicles to zero
            print('SET DOCK TETS DIFF TO ZERO:', list(dock_ves_init_tets.values()))
            sim.TETS(dock_ves_init_tets.values()).ves.Dcst = 0

            ####################################################################
            # non-recycling vesicles
            ####################################################################
            # add non-recycling vesicles to bouton
            ves_n = 250
            nonrec_ves = []
            for v in range(ves_n):
                vesref = sim.cytosol.addVesicle(ves)
                nonrec_ves.append(vesref)
                #add synapsin to non-recycling vesicles
                #Add vesicle membrane species
                vesref('surf').ves_syn1_site.Count = 10#75
                vesref('surf').ves_CXN_site.Count = 1000
                vesref('surf').ves_m13_site.Count = 1000
                vesref('surf').ves_m18_site.Count = 1000
                vesref('surf').ves_syndap_site.Count = 1000
                vesref('surf').ves_asnap_site.Count = 1000
                vesref('surf').ves_nsf_site.Count = 1000
                #adding Rab3 directly for release prob tests
                vesref('surf').Rab3.Count = 10
                vesref('surf').SYB.Count = 69
                vesref('surf').syt.Count = 27
                vesref('in').REL_IND.Count = 1
            print("Added synapsin to non-recycling vesicles")
            nonrec_ves_init = copy.copy(nonrec_ves)

            #add recycling vesicles without synapsin
            rec_ves_n = 50
            rec_ves = []
            for r in range(rec_ves_n):
                vesref = sim.cytosol.addVesicle(ves)
                rec_ves.append(vesref)
                #Add vesicle membrane species
                vesref('surf').ves_CXN_site.Count = 1000
                vesref('surf').ves_m13_site.Count = 1000
                vesref('surf').ves_m18_site.Count = 1000
                vesref('surf').ves_syndap_site.Count = 1000
                vesref('surf').ves_asnap_site.Count = 1000
                vesref('surf').ves_nsf_site.Count = 1000
                #adding Rab3 directly for release prob tests
                vesref('surf').Rab3.Count = 10 #Takamori 2006
                vesref('surf').SYB.Count = 69
                vesref('surf').syt.Count = 27
                #Add glutamate to each vesicle
                vesref('in').REL_IND.Count = 1
                # vesref('in').glu.Count = 400
            print('Added recycling vesicles:', len(rec_ves))
            #save recycling vesicles at init
            rec_ves_init = rec_ves + dock_ves_init
            print('Rec ves ind: ', rec_ves_init)

            RAFT_N = 40
            raft_count = 0
            while raft_count < RAFT_N :
                raftref = sim.TRI(random.choice(endo_tris)).addRaft(raft)
                if raftref.exists():
                    raftref.Rab3.Count = 10 # Takamori 2006
                    raftref.SYB.Count = 69
                    raftref.syt.Count = 27
                    raftref.DYNp.Count = 1
                    raftref.ves_CXN_site.Count = 1000
                    raftref.ves_m13_site.Count = 1000
                    raftref.ves_m18_site.Count = 1000
                    raftref.ves_syndap_site.Count = 1000
                    raftref.ves_asnap_site.Count = 1000
                    raftref.ves_nsf_site.Count = 1000
                    raft_count += 1
            raft_count_init = RAFT_N

            # populate docking sites with RIM
            sim.TRIS(dock_tris).RIM_M13.Count = 4

            if cachan_static:
                # add 4 Ca channel to each dock triangle
                sim.TRIS(dock_tris).CaPchan[CaPc, CaPc, CaPc].Count = 4
            else:
                # add 15 Ca channels randomly over active zone
                indices = random.sample(azone_tris.indices)
                import mpi4py.MPI
                indices = mpi4py.MPI.COMM_WORLD.bcast(indices, root=0)
                sim.TRIS(indices).CaPchan[CaPc, CaPc, CaPc].Count = 1

            print(f'Added {sim.memb.RIM_M13.Count} RIM to membrane')
            print(f'Added {sim.memb.CaPchan[CaPc, CaPc, CaPc].Count} CaP_m0 to membrane')

            #scale membrane proteins by 1.6 (area of mesh membrane/cylinder membrane)
            mem_factor = 1.6

            sim.cytosol.Ca.Conc = 0.045e-6
            # total CB in rat HC: 1.98e-6
            sim.cytosol.CBhi.Conc = 0.99e-6 # i.e. 1/2 of total CB molarity (1.98*10e-6)
            sim.cytosol.CBlo.Conc = 0.99e-6 # i.e. 1/2 of total CB molarity (1.98*10e-6)

            # total CR in rat HC: 2.47e-6
            # WE MODEL TWO IDENTICAL PAIRS OF COOPERATIVE BINDING SITES AS ONE SPECIES.
            # SO WE SET ITS CONCENTRATION TO 4/5 OF THE TOTAL CONCENTRATION OF CR. THE REMAINING 1/5 IS THE INDEPENDENT SITE/
            sim.cytosol.CRTT.Conc = 0.1976e-6  # 4/5 of total CR concentration (2 pairs of cooperative sites)
            sim.cytosol.CRind.Conc = 0.494e-6  # 1/5 of total CR

            # total PV in rat HC: 4.55e-6
            sim.cytosol.PV.Conc = 4.55e-6

            sim.cytosol.PV_Ca.Conc = 8.4e-6
            sim.cytosol.MgPV.Conc = 30.45e-6

            #Add AC18 and PKAinact
            sim.cytosol.AC18.Conc = 0.2e-6#0.5e-6
            sim.cytosol.R2C2.Conc = 1e-6#0.5e-6

            SERCA_ro = 1000*1e12 # Bartol 2015. Computational reconstitution...
            PMCA_ro = 180*1e12 # Bartol 2015. Computational reconstitution...
            SERCA_count = ERmemb.Area * SERCA_ro
            PMCA_count = mem_area * PMCA_ro
            sim.ERmemb.SERCA.Count = SERCA_count
            sim.memb.PMCA_P0.Count = PMCA_count
            print('Added SERCA:', sim.ERmemb.SERCA.Count)
            print('Added PMCA:', sim.memb.PMCA_P0.Count)

            sim.cytER.Ca.Conc = 150e-6
            sim.cytER.Ca.Clamped = True # clamped means the conc won't change as simulation runs.

            sim.cytosol.M18.Conc = 28.4e-6
            print('Added M18 to cytosol:', sim.cytosol.M18.Count)
            sim.cytosol.M13.Conc = 10.36e-6
            print('Added M13 to cytosol:', sim.cytosol.M13.Count)
            sim.cytosol.CXN.Conc = 16.59e-6
            print('Added complexin to cytosol:', sim.cytosol.CXN.Count)

            sim.cytosol.Rab3.Conc =  0#125.85e-6
            print('Added Rab3 to cytosol:', sim.cytosol.Rab3.Count)
            sim.cytosol.aSNAP.Conc = 7.68e-6
            print('Added aSNAP to cytosol:', sim.cytosol.aSNAP.Count)
            sim.cytosol.NSF.Conc = 27.14e-6
            print('Added NSF to cytosol:', sim.cytosol.NSF.Count)
            sim.cytosol.synapsin.Conc = 156e-6
            synapsin_total = sim.cytosol.synapsin.Count
            print('Added synapsins to cytosol:', synapsin_total)

            #add PP2A to cytocol
            sim.cytosol.PP2A.Conc = 0.5e-6

            sim.memb.SYX.Count = 20096*mem_factor
            sim.memb.SNP25.Count = 26686*mem_factor

            if cam_present:
                sim.cytosol.CaM[N0, C0].Conc = 6e-05 #Wilhelm 2014
                print("CALMODULIN PRESENT")
            else:
                sim.cytosol.CaM[N0, C0].Conc = 0
                print("CALMODULIN NOT PRESENT")
            sim.cytosol.CaN.Conc = 5e-06
            sim.cytosol.DYNp.Count = 0 #20 #plus 40 in initial rafts
            sim.cytosol.CDK5.Conc = 1e-06
            sim.cytosol.Rab3.Conc = 0.00012585

            sim.cytosol.SYN1.Conc = 2.137e-05
            print('Added Syndapin1 to cytosol:', sim.cytosol.SYN1.Count)

            #add TOMOSYN
            sim.cytosol.TOMO_p.Count = tomo_init

            sim.TETS(actin_tets).actin.Count = 1

            fusion_events_tot = 0
            #check release indicator in cleft pre and post run
            rel_ind_count_pre = 0
            rel_ind_count_post = 0
            #create empty list to contain recycled vescicles
            recycled_ves_free = []
            recycled_ves_all = []
            #dictionary of Rafts to Vesicle indices (connecting endocytosis)
            raft_to_vesicle = {}
            #record total raft endocytosis events
            total_raft_endo = 0
            total_time = 0

            rel_ind_count_pre_ap = sim.cleft.REL_IND.Count
            all_ves_pre_ap = VesicleList(sim.cytosol.ves)

            #####SCALE DIFFUSION OF SYNTAXIN AND SNAP25#####################
            sim.TRIS(azone_tris).diff_SYX.D = 0.2429e-12*snare_diff_scale
            sim.TRIS(azone_tris).diff_SYX_M18.D = 0.2429e-12*snare_diff_scale
            sim.TRIS(azone_tris).diff_SNP25.D = 0.651e-12*snare_diff_scale

            action_potentials_fired = 0
            #Save release probs
            ap_n = len(pulse_times)
            rel_probs = np.zeros(ap_n) #10 APs
            #Save release dock sites
            #rel_tris = [] #10 APs
            #count fusions from "recycling" vs "reserve" pool
            REC_FUSIONS = []
            RES_FUSIONS = []
            #get all initial vesicle indices
            init_vesicle_ind = VesicleList(sim.cytosol.ves)
            #build list for vesicles (f)used from initial pools
            init_vesicles_used = []
            init_rec_vesicles_used = []
            init_nonrec_vesicles_used = []
            new_vesicles_used = []
            fused_ind_all = []
            #record total recycled vesicles
            recycled_vesicles_tot = 0
            #list of new rafts
            new_rafts = []
            filled_rafts = []

            #create array of time points
            tpnts = np.arange(0.0, INT, DT)
            #find number of time points
            ntpnts = tpnts.shape[0]

        print('Number of timepoints =', ntpnts)

        if MPI.rank == 0 and exportParam:
            from steps.utils import ExportParameters
            ExportParameters(sim, 'Vesicle_Cycle', method='pdf',
                hideColumns=['Defined in', 'valence Units'],
                unitsToSimplify=[
                    'uM', 'uM^-1 s^-1', 'mV', 'um^2 s^-1', 'pS', 'um', 'um^2'
                ], numPrecision=5,
            )

        for i, t in enumerate(tpnts):

            if t >= add_tethers_time:
                #create tether paths for each dock triangle
                print("CREATING TETHER PATHS...")
                print('ADDING VESICLES TO TETHER PATHS...')
                tether_pull_rate_m_s = 200.0 / 60 * 1e-06 # um / min -> m / s
                print(f'TETHER RATE: {tether_pull_rate_m_s} m/s')
                tether_paths = []
                for tri in dock_tris:
                    path = sim.addVesiclePath(f'actin_{tri.idx}')
                    pos1, pos2 = dock_tri_to_tether_path[tri]
                    p2 = path.addPoint(pos2)
                    p1 = path.addPoint(pos1)
                    print(f'{path.idx} PATH LENGTH: {np.linalg.norm(pos1 - pos2) * 1e9} nm')
                    #connect path points
                    path.addBranch(p2,  {p1: 1.0})
                    path.addVesicle(ves, tether_pull_rate_m_s, dependencies=8 * Rab3.v, stoch_stepsize=2e-9)
                    print('Added vesicle to path:', path.idx)
                print('VESICLES ADDED TO TETHER PATHS')
                print("TETHER PATHS CREATED")
                add_tethers_time = math.inf

            print(t)
            if t >= pka_act_time:
                pka_act_time = math.inf
                sim.cytosol.PKA.Count = 60
            if t >= camkii_deact_time:
                camkii_deact_time = math.inf
            #get all vesicles and rafts before the next time step to test for new vesicles by endoctosis
            all_ves_prerun = VesicleList(sim.cytosol.ves)
            raft_prerun = RaftList(sim.memb.raft)
            raft_positions_prerun = {raftref: np.array(raftref.Pos) for raftref in raft_prerun}
            #get number of exocytosis events pre-run
            exo_extent_pre = sum(sim.ALL(Exocytosis).Extent)
            #count release indicators in cleft pre run
            rel_ind_count_pre = sim.cleft.REL_IND.Count
            print('Exo_extent_pre =', rel_ind_count_pre)

            ##########RUN TIME STEP###########
            start = time.time()
            sim.run(t)
            end = time.time()

            #count release indicators in cleft post run
            rel_ind_count_post = sim.cleft.REL_IND.Count
            print("rel_ind_count_post:", rel_ind_count_post)
            exo_extent_post = sum(sim.ALL(Exocytosis).Extent)
            print('Exo_extent_post = ', rel_ind_count_post)
            print("ACTION POTENTIALS FIRED: ", action_potentials_fired)
            fusion_events = rel_ind_count_post - rel_ind_count_pre
            print("fusion_events:", fusion_events)
            fusion_events_tot += fusion_events
            print('fusion_events_tot =', fusion_events_tot)

            time_took = end - start
            print('Took:', time_took)
            if i != 0:
                total_time += time_took
                mean_time = total_time / i
                print('Mean time:', mean_time)
                print('Time remaining (hours):', (ntpnts - i) * mean_time / 3600)

            #move any docked vesicles to their docking site
            #get list of all vesicles (some may have appeared by endocytosis)
            #get list of newly formed vesicles
            all_ves = VesicleList(sim.cytosol.ves)
            new_ves = [v for v in all_ves if v not in all_ves_prerun]
            for v in new_ves:
                #add REL_IND to vesicle
                v('in').REL_IND.Count = 1
            #check that all docked_ves still exist (i.e. haven't fused)
            for v in all_ves_prerun:
                if v not in all_ves:
                    FUSED_VES.append(v)

            if use_efield:
                if round(t * 1000) in pulse_times:
                    sim.membrane.Potential = 40e-3
                    print('Action potential ON')
                if round(t * 1000) in pulse_times_off:
                    sim.membrane.Potential = init_pot
                    print('Action potential OFF')
                if round(t * 1000) in check_release_times_pre:
                    rel_ind_count_pre_ap = sim.cleft.REL_IND.Count
                    all_ves_pre_ap = VesicleList(sim.cytosol.ves)
                if round(t * 1000) in check_release_times:
                    all_ves_post_ap = VesicleList(sim.cytosol.ves)
                    rel_ind_count_post_ap = sim.cleft.REL_IND.Count
                    fusions = rel_ind_count_post_ap - rel_ind_count_pre_ap
                    if fusions != 0:
                        rel_probs[action_potentials_fired] = 1
                        #get all current vesicles to check for fusion
                        fused_ind = []
                        for v in all_ves_pre_ap:
                            if (v not in all_ves_post_ap):
                                fused_ind.append(v)
                                fused_ind_all.append(v)
                                #add fused vesicle to used vesicle lists
                        print('FUSED VESICLES:', fused_ind)
                        for v in fused_ind:
                            if v in init_vesicle_ind:
                                init_vesicles_used.append(v)
                                print(v, 'is an init vesicle type')
                            if v in rec_ves_init:
                                init_rec_vesicles_used.append(v)
                                print(v, 'is a recycling vesicle type')
                            if v in nonrec_ves_init:
                                init_nonrec_vesicles_used.append(v)
                                print(v, 'is a reserve vesicle type')
                            if v not in init_vesicle_ind:
                                new_vesicles_used.append(v)
                                print(v, 'is a new vesicle type')
                            print(v, 'is not in the list', nonrec_ves_init)
                            print('init_nonrec_vesicles_used', init_nonrec_vesicles_used)
                                #add indices to ves_docktris
                        #create a new Rafts for each fused vesicle
                        count = 0
                        attempts = 0
                        print('CREATING NEW RAFT')
                        while count < fusions:
                            rand_tri = random.choice(endo_tris)
                            print('RAND TRI SELECTED:', rand_tri)
                            raftref = sim.TRI(rand_tri).addRaft(raft)
                            print('INDEX: ', raftref)
                            if (raftref.exists()):
                                print('ADDING SPECS TO NEW RAFT')
                                raftref.AP180.Count = 69
                                raftref.AP2.Count = 27
                                raftref.Rab3.Count = 10
                                raftref.DYN_AD.Count = 1
                                raftref.ves_CXN_site.Count = 1000
                                raftref.ves_m13_site.Count = 1000
                                raftref.ves_m18_site.Count = 1000
                                raftref.ves_syndap_site.Count = 1000
                                raftref.ves_asnap_site.Count = 1000
                                raftref.ves_nsf_site.Count = 1000
                                count += 1
                                new_rafts.append(raftref)
                            else:
                                attempts+=1
                                print('RAND TRI SELECTION ATTEMPTS:', attempts)
                                continue
                    else:
                        rel_probs[action_potentials_fired] = 0
                    action_potentials_fired += 1

            # generate list of current docked vesicles
            docked_ves = []
            for v in all_ves:
                #check immobility
                if v.Immobility != 0: # ONLY DOCKED VESICLES ARE IMMOBILE IN rec_ves_ind
                    surf_actin_n = sum(v('surf').LIST(synapsin_actin, Rab3_TOMO_p_synapsin_actin, Rab3_TOMO_synapsin_actin).Count)
                    if surf_actin_n == 0: #Check not a non-recycling vesicle bound to actin
                        ves_pos = v.Pos
                        min_dist_ind_2 = np.argmin([np.linalg.norm(tri.center - ves_pos) for tri in dock_tris])
                        dock_pos = dock_positions[min_dist_ind_2]
                        # only move vesicle if closer enough to dock site (i.e. not actin bound)
                        if np.linalg.norm(ves_pos - dock_pos) < 100e-9:
                            print("Attempted new position:", dock_pos)
                            v.Pos = dock_pos
                            new_pos = v.Pos
                            # test vesicle has actually been moved to dock site (i.e. not blocked by another vesicle)
                            if np.linalg.norm(new_pos - dock_pos) < 1e-09:
                                docked_ves.append(v)
                            print('Ves current pos:', ves_pos)
                            print('Ves new pos:', new_pos)
            print('Dock ves ind:', docked_ves)

            #remove dock ves from nonrec_ves (if docked)
            for d in docked_ves:
                if d in nonrec_ves:
                    nonrec_ves.remove(d)
                if d in rec_ves:
                    rec_ves.remove(d)
            for d in FUSED_VES:
                if d in nonrec_ves:
                    nonrec_ves.remove(d)
                if d in nonrec_ves_init:
                    if d not in RES_FUSIONS:
                        RES_FUSIONS.append(d)
                if d in rec_ves:
                    rec_ves.remove(d)
                if d in rec_ves_init:
                    if d not in REC_FUSIONS:
                        REC_FUSIONS.append(d)

            # create dictionary of docked vesicles and their dock sites (ves -> dock_tri)
            ves_docktris = {v: min(dock_tris, key=lambda t: np.linalg.norm(t.center - v.Pos)) for v in docked_ves}

            #find Raft each new vesicle came from.
            rafts = RaftList(sim.memb.raft)
            rafts_endocytosed = []
            for r in raft_prerun:
                if r not in rafts:
                    rafts_endocytosed.append(r)
            total_raft_endo += len(rafts_endocytosed)
            print('Rafts endocytosed:', total_raft_endo)
            for v in new_ves:
                if len(rafts_endocytosed) > 0:
                    closest_raft = min(rafts_endocytosed, key=lambda r: np.linalg.norm(raft_positions_prerun[r] - v.Pos))
                    raft_to_vesicle[closest_raft] = v
                    rafts_endocytosed.remove(closest_raft)
                    print(f'Raft {closest_raft} formed vesicle {v}')

            #######RELEASE PROBABILITY TEST####################################################
            if t >= rel_prob_test_time:
                rel_prob_test_time = math.inf
                # add SNARE_SYT_CXN complexes to each docked vesicle
                for v in docked_ves:
                    print('Priming vesicle:', v)
                    v('surf').SNARE_syt_CXN.Count = 4
                    # remove RIM and related structures from vesicle dock sites
                    v('surf').RIM_M13_Rab3.Count = 0
                    v('surf').RIM_M13_Rab3_SYX_M18.Count = 0
                    v('surf').RIM_M13_Rab3_SYX_M18_SNP25.Count = 0
                    v('surf').SNARE.Count = 0
                    v('surf').SNARE_syt.Count = 0
                    v('surf').SNARE_syt_CXN_Ca.Count = 0
                    v('surf').SNARE_syt_CXN_Ca2.Count = 0
                    v('surf').SNARE_syt_CXN_Ca3.Count = 0
                    v('surf').SNARE_syt_CXN_bCa.Count = 0
                    v('surf').SNARE_syt_CXN_bCa2.Count = 0
                    v('surf').SNARE_syt_CXN_Ca_bCa.Count = 0
                    v('surf').SNARE_syt_CXN_Ca_bCa2.Count = 0
                    v('surf').SNARE_syt_CXN_Ca2_bCa.Count = 0
                    v('surf').SNARE_syt_CXN_Ca2_bCa2.Count = 0
                    v('surf').SNARE_syt_CXN_Ca3_bCa.Count = 0
                    v('surf').SNARE_syt_CXN_Ca3_bCa2.Count = 0
                # set normal diffusion rate for dock vesicles
                for tet in dock_ves_init_tets.values():
                    sim.TET(tet).ves.Dcst = ves_diff_k
                    print('SET DOCK TET DIFF TO ZERO:', tet)

            print('Dock ves ind:', docked_ves)
            ###################################################################################

            #shift all SNARE complexes to pole of vesicle (closest to membrane)
            for v in docked_ves:
                xd, yd, zd = ves_docktris[v].center - v.Pos
                theta = np.arctan2(yd, xd)
                phi = np.arccos(zd / ves_radius)
                v('surf').SNARE.PosSpherical = [theta, phi]
                v('surf').SNARE_syt.PosSpherical = [theta, phi]
                v('surf').SNARE_syt_CXN.PosSpherical = [theta, phi]
                v('surf').SNARE_syt_CXN_Ca.PosSpherical = [theta, phi]
                v('surf').SNARE_syt_CXN_Ca2.PosSpherical = [theta, phi]
                v('surf').SNARE_syt_CXN_Ca3.PosSpherical = [theta, phi]
                v('surf').SNARE_syt_CXN_bCa.PosSpherical = [theta, phi]
                v('surf').SNARE_syt_CXN_bCa2.PosSpherical = [theta, phi]
                v('surf').SNARE_syt_CXN_Ca_bCa.PosSpherical = [theta, phi]
                v('surf').SNARE_syt_CXN_Ca_bCa2.PosSpherical = [theta, phi]
                v('surf').SNARE_syt_CXN_Ca2_bCa.PosSpherical = [theta, phi]
                v('surf').SNARE_syt_CXN_Ca3_bCa.PosSpherical = [theta, phi]
                v('surf').SNARE_syt_CXN_Ca3_bCa2.PosSpherical = [theta, phi]
                v('surf').RIM_M13_Rab3.PosSpherical = [theta, phi]
                v('surf').RIM_M13_Rab3_SYX_M18.PosSpherical = [theta, phi]
                v('surf').RIM_M13_Rab3_SYX_M18_SNP25.PosSpherical = [theta, phi]

            # find all recycling vesicles not docked
            rec_ves_free = [v for v in all_ves if (v not in docked_ves) and (v not in nonrec_ves)]

            #separate non-recycling vesicles into free and clustered
            nonrec_ves_free = []
            nonrec_ves_clust = []
            actin_bound_ves = []
            total_link_specs = 0
            for v in nonrec_ves:
                dimer_count = getVesicleLinkSpecs(v)
                total_link_specs += dimer_count
                if dimer_count < 1:
                    nonrec_ves_free.append(v)
                else:
                    nonrec_ves_clust.append(v)
                ves_actin_count = sum(v('surf').LIST(synapsin_actin, Rab3_TOMO_p_synapsin_actin).Count)
                if ves_actin_count != 0:
                    actin_bound_ves.append(v)
            #find recycled vesicles
            recycled_ves_free += [v for v in new_ves if (v not in docked_ves) and (v not in FUSED_VES)]
            #check no recycled vesicles have docked
            rec_docked = [v for v in recycled_ves_free if (v in docked_ves) or (v in FUSED_VES)]
            for v in rec_docked:
                recycled_ves_free.remove(v)
            #remove recycled vesicles in rec_ves_free
            ves_free_to_remove = [v for v in rec_ves_free if v in recycled_ves_free]
            for v in ves_free_to_remove:
                rec_ves_free.remove(v)

            #concatenate rec_ves_free and nonrec_ves_free to get all free vesicles
            all_ves_free = rec_ves_free + nonrec_ves_free

            print('Number of docked vesicles =', len(docked_ves))
            print('Number of recycled vesicles =', len(recycled_ves_free))
            #sanity check
            print('Total vesicles:', len(all_ves))
            ves_sum = docked_ves + all_ves_free + recycled_ves_free
            print('Summed vesicles:', len(ves_sum) + len(nonrec_ves_clust))

            print('rec_ves_free:', len(rec_ves_free))
            print('nonrec_ves_free:', len(nonrec_ves_free))
            print('all_ves_free:', len(all_ves_free))

            rafts_n = len(RaftList(sim.memb.raft))
            print("N rafts:", rafts_n)

            print('Number of SYB in rafts:', sim.memb.raft.SYB.Count)

            print('SYB in membrane =', sim.memb.SYB.Count)

            #####calculate SYX and SNAP25 in AZONE####
            print('SYX in AZone:', sum(sim.TRIS(azone_tris).LIST(SYX, SYX_M18).Count))
            print('SNAP25 in AZone:', sum(sim.TRIS(azone_tris).SNP25.Count))

            endo_events = raft_count_init - rafts_n
            print('Total endo events:', endo_events)
            print('Ca4CaM count =', sim.cytosol.CaM[N2, C2].Count)
            print('CaN_CaM count =', sim.cytosol.CaN_CaM.Count)
            print('DYNpp count =', sim.cytosol.DYNpp.Count)
            print('DYN count =', sim.cytosol.DYN.Count)
            print('DYN MEM count =', sum(sim.memb.LIST(DYN_SYN1, DYN_SYN1_CDK5).Count))

            synapsin_n = sum(sim.cytosol.LIST(synapsin, synapsin_actin).Count)
            print('Bound synapsin:', synapsin_total - synapsin_n)
            synapsin_cluster_perc = (1 - (synapsin_n / synapsin_total)) * 100
            print('Percentage SYNAPSIN in cluster:', synapsin_cluster_perc)
            print('TOTAL CXN in CYTOSOL:', sim.cytosol.CXN.Count)
            print('TOTAL CXN in AXON:', sim.axon.CXN.Count)

            print("FUSION COUNT:", fusion_count)
            print('TOTAL Ca in CYTOSOL:', 1e06 * sim.cytosol.Ca.Conc)
            print('TOTAL CaM[N2, C2] in CYTOSOL:', sim.cytosol.CaM[N2, C2].Count)
            print('TOTAL AC18 in CYTOSOL:', sim.cytosol.AC18.Count)
            print('TOTAL AC18_CaM in CYTOSOL:', sim.cytosol.AC18_CaM.Count)
            print('TOTAL cAMP in CYTOSOL:', sim.cytosol.cAMP.Count)
            print('TOTAL R2C2 in CYTOSOL:', sim.cytosol.R2C2.Count)
            print('TOTAL R2C2_cAMP in CYTOSOL:', sim.cytosol.R2C2_cAMP.Count)
            print('TOTAL PKA in CYTOSOL:', sim.cytosol.PKA.Count)

            print('Number of vesicles bound to actin =', len(actin_bound_ves))
            syn_dimers_n = sum(sim.cytosol.ves('surf').ALL(LinkSpecies).Count)
            print('Number of synapsin dimers =', syn_dimers_n)

            print('NON-CLUSTERED VESICLES:', len(nonrec_ves_free))
            print('CLUSTERED VESICLES:', len(nonrec_ves_clust))
            print('TOTAL PHOSPHO SYNAPSIN in CYTOSOL:', sim.cytosol.synapsin_p.Count)

            if MPI.rank == 0:
                sd_val = recorded_values.data[-1, -1, 16]
                print('Number of disassembled cis-SNAREs:', sd_val)

            print('Fused vesicles: ', sorted(fused_ind_all, key=lambda vr: vr.idx))

            if MPI.rank == 0:
                print('Total phosho synapsins: ', recorded_values.data[-1, -1, 21])
            print('Total rafts:', sim.memb.raft.Count)

            total_dynamin_cytosol = sum(sim.cytosol.LIST(DYN, DYNp, DYN_SYN1).Count)
            print('Total dynamin in cytosol:', total_dynamin_cytosol)
            total_dynamin_membrane = sum(sim.memb.LIST(DYN, DYNp, DYN_SYN1).Count)
            print('Total dynamin in membrane:', total_dynamin_membrane)

            # Dynamin in rafts
            DYNpp_raft = sum(sim.memb.raft.LIST(DYNpp, CaN_DYNpp).Count)
            DYNp_raft = sum(sim.memb.raft.LIST(DYNp, CaN_DYNp).Count)
            DYN_raft = rs.memb.raft.DYN.Count
            DYN_SYN1_raft = sum(sim.memb.raft.LIST(DYN_SYN1, DYN_SYN1_CDK5).Count)
            total_dynamin_rafts = DYNpp_raft + DYNp_raft + DYN_raft + DYN_SYN1_raft
            print('Total dynamin in rafts:', total_dynamin_rafts)
            print('DYNpp in rafts:', DYNpp_raft)
            print('DYNp in rafts:', DYNp_raft)
            print('DYN in rafts:', DYN_raft)
            print('DYN_SYN1 in rafts:', DYN_SYN1_raft)

            # check dynamin saturation and add to raft as necessary (only one per time step in chronological order)
            print('Total new rafts: ', len(new_rafts))
            print('Total filled rafts: ', len(filled_rafts))
            if total_dynamin_rafts < 60:
                for raftref in new_rafts:
                    ves_protein_count = raftref.SYB.Count + raftref.syt.Count
                    if ves_protein_count > 93:
                        raftref.DYNp.Count = 1
                        raftref.DYN_AD.Count = 0
                        new_rafts.remove(raftref)
                        filled_rafts.append(raftref)
                        break

            #maintain RIM_M13 in dock tris
            dock_tris_occupied = list(ves_docktris.values())
            for tri in dock_tris:
                if tri not in dock_tris_occupied:
                    sim.TRI(tri).RIM_M13.Count = 4
            for tri in azone_and_endo_tris:
                if tri not in dock_tris:
                    sim.TRI(tri).RIM_M13.Count = 0

            #separate vesicles into RECYCLING and RESERVE POOLS
            #create lists of Rab3 counts
            RAB3_dict = sim.VESICLES(sim.cytosol.ves('surf')).Rab3.Count
            VES_RAB3_ALL = list(RAB3_dict.values())
            VES_RAB3_REC = [RAB3_dict[v.idx] for v in rec_ves]
            VES_RAB3_RES = [RAB3_dict[v.idx] for v in nonrec_ves]
            custom_values.save([
                len(nonrec_ves_clust),
                len(nonrec_ves_free),
                len(docked_ves),
                len(init_vesicles_used),
                len(init_rec_vesicles_used),
                len(init_nonrec_vesicles_used),
                len(new_vesicles_used),
                len(init_rec_vesicles_used) + len(init_nonrec_vesicles_used) + len(new_vesicles_used),
                total_link_specs,
                VES_RAB3_ALL,
                VES_RAB3_REC,
                VES_RAB3_RES,
                list(rel_probs),
            ])

            if LOG_FILE is not None:
                LOG_FILE.flush()

    if checkpoint_mode in ['checkpoint', 'both']:
        if checkpoint_mode == 'both':
            checkpoint_path = checkpoint_path_2
        # Checkpoint the state of the solver
        sim.checkpoint(checkpoint_path)

        # Also checkpoint the state of the script variables
        def checkpoint_vesraft_list(f, lst):
            pickle.dump([(vesref._type, vesref.idx) for vesref in lst], f)

        pkl_path = checkpoint_path + '_script_state.pkl'
        with open(pkl_path, 'wb') as f:
            pickle.dump(dock_positions, f)
            pickle.dump({(vesref._type, vesref.idx): tet.idx for vesref, tet in dock_ves_init_tets.items()}, f)
            veslists = [
                nonrec_ves, nonrec_ves_init, rec_ves, recycled_ves_free, recycled_ves_all, init_vesicle_ind, init_vesicles_used,
                init_rec_vesicles_used, init_nonrec_vesicles_used, new_vesicles_used, fused_ind_all, REC_FUSIONS, RES_FUSIONS,
                rec_ves_init, all_ves_pre_ap,
            ]
            for lst in veslists:
                checkpoint_vesraft_list(f, lst)

            checkpoint_vesraft_list(f, new_rafts)
            checkpoint_vesraft_list(f, filled_rafts)
            checkpoint_vesraft_list(f, raft_to_vesicle.keys())
            checkpoint_vesraft_list(f, raft_to_vesicle.values())
            pickle.dump((
                fusion_events_tot, rel_ind_count_pre, rel_ind_count_post, total_raft_endo,
                total_time, action_potentials_fired, ap_n, rel_probs, recycled_vesicles_tot,
                raft_count_init, synapsin_total, rel_ind_count_pre_ap,
            ), f)

if LOG_FILE is not None:
    LOG_FILE.close()
