Exp1_DopaEffect15MAR2024.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db2bcb7a-16b6-448a-b16f-410ab7cb46f7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Start time:  2024-03-15 13:26:28.453214\n",
      "\n",
      "Creating network of 3 cell populations on 1 hosts...\n",
      "  Number of cells on node 0: 30 \n",
      "  Done; cell creation time = 0.01 s.\n",
      "Making connections...\n",
      "  Number of connections on node 0: 870 \n",
      "  Number of synaptic contacts on node 0: 1450 \n",
      "  Done; cell connection time = 0.11 s.\n",
      "Adding stims...\n",
      "  Number of stims on node 0: 100 \n",
      "  Done; cell stims creation time = 0.01 s.\n",
      "Recording 60 traces of 2 types on node 0\n",
      "\n",
      "Running simulation using NEURON for 30000.0 ms...\n",
      "  Done; run time = 224.63 s; real-time ratio: 0.13.\n",
      "\n",
      "Gathering data...\n",
      "  Done; gather time = 0.20 s.\n",
      "\n",
      "Analyzing...\n",
      "  Cells: 30\n",
      "  Connections: 900 (30.00 per cell)\n",
      "  Synaptic contacts: 1540 (51.33 per cell)\n",
      "  Spikes: 6295 (6.99 Hz)\n",
      "  Simulated time: 30.0 s; 1 workers\n",
      "  Run time: 224.63 s\n",
      "Saving output as patient_0_data.pkl ... \n",
      "Finished saving!\n",
      "  Done; saving time = 3.42 s.\n",
      "Plotting recorded cell traces ... cell\n",
      "Preparing spike data...\n",
      "Plotting raster...\n",
      "Unable to run analyzeRaster from sim.plotting and sim.analysis. Reason: module 'netpyne.analysis' has no attribute 'analyzeRaster'\n",
      "  Done; plotting time = 4.91 s\n",
      "\n",
      "Total time = 233.34 s\n",
      "\n",
      "Start time:  2024-03-15 13:30:21.805697\n",
      "\n",
      "Creating network of 3 cell populations on 1 hosts...\n",
      "  Number of cells on node 0: 30 \n",
      "  Done; cell creation time = 0.00 s.\n",
      "Making connections...\n",
      "  Number of connections on node 0: 870 \n",
      "  Number of synaptic contacts on node 0: 1450 \n",
      "  Done; cell connection time = 0.10 s.\n",
      "Adding stims...\n",
      "  Number of stims on node 0: 100 \n",
      "  Done; cell stims creation time = 0.02 s.\n",
      "Recording 60 traces of 2 types on node 0\n",
      "\n",
      "Running simulation using NEURON for 30000.0 ms...\n",
      "  Done; run time = 249.01 s; real-time ratio: 0.12.\n",
      "\n",
      "Gathering data...\n",
      "  Done; gather time = 0.29 s.\n",
      "\n",
      "Analyzing...\n",
      "  Cells: 30\n",
      "  Connections: 900 (30.00 per cell)\n",
      "  Synaptic contacts: 1540 (51.33 per cell)\n",
      "  Spikes: 6295 (6.99 Hz)\n",
      "  Simulated time: 30.0 s; 1 workers\n",
      "  Run time: 249.01 s\n",
      "Saving output as patient_0_data.pkl ... \n",
      "Finished saving!\n",
      "  Done; saving time = 3.54 s.\n",
      "Plotting recorded cell traces ... cell\n",
      "Preparing spike data...\n",
      "Plotting raster...\n",
      "Unable to run analyzeRaster from sim.plotting and sim.analysis. Reason: module 'netpyne.analysis' has no attribute 'analyzeRaster'\n",
      "  Done; plotting time = 5.13 s\n",
      "\n",
      "Total time = 258.14 s\n",
      "Patient 0:\n",
      "Burst Duration starting at 5 seconds: 668 ms\n",
      "Burst Frequency: 4e-05 Hz\n",
      "\n",
      "Start time:  2024-03-15 13:34:39.948796\n",
      "\n",
      "Creating network of 3 cell populations on 1 hosts...\n",
      "  Number of cells on node 0: 30 \n",
      "  Done; cell creation time = 0.00 s.\n",
      "Making connections...\n",
      "  Number of connections on node 0: 870 \n",
      "  Number of synaptic contacts on node 0: 1450 \n",
      "  Done; cell connection time = 0.10 s.\n",
      "Adding stims...\n",
      "  Number of stims on node 0: 100 \n",
      "  Done; cell stims creation time = 0.01 s.\n",
      "Recording 60 traces of 2 types on node 0\n",
      "\n",
      "Running simulation using NEURON for 30000.0 ms...\n",
      "  Done; run time = 231.06 s; real-time ratio: 0.13.\n",
      "\n",
      "Gathering data...\n",
      "  Done; gather time = 0.28 s.\n",
      "\n",
      "Analyzing...\n",
      "  Cells: 30\n",
      "  Connections: 900 (30.00 per cell)\n",
      "  Synaptic contacts: 1540 (51.33 per cell)\n",
      "  Spikes: 6452 (7.17 Hz)\n",
      "  Simulated time: 30.0 s; 1 workers\n",
      "  Run time: 231.06 s\n",
      "Saving output as patient_1_data.pkl ... \n",
      "Finished saving!\n",
      "  Done; saving time = 3.44 s.\n",
      "Plotting recorded cell traces ... cell\n",
      "Preparing spike data...\n",
      "Plotting raster...\n",
      "Unable to run analyzeRaster from sim.plotting and sim.analysis. Reason: module 'netpyne.analysis' has no attribute 'analyzeRaster'\n",
      "  Done; plotting time = 4.86 s\n",
      "\n",
      "Total time = 239.81 s\n",
      "\n",
      "Start time:  2024-03-15 13:38:39.766654\n",
      "\n",
      "Creating network of 3 cell populations on 1 hosts...\n",
      "  Number of cells on node 0: 30 \n",
      "  Done; cell creation time = 0.01 s.\n",
      "Making connections...\n",
      "  Number of connections on node 0: 870 \n",
      "  Number of synaptic contacts on node 0: 1450 \n",
      "  Done; cell connection time = 0.10 s.\n",
      "Adding stims...\n",
      "  Number of stims on node 0: 100 \n",
      "  Done; cell stims creation time = 0.00 s.\n",
      "Recording 60 traces of 2 types on node 0\n",
      "\n",
      "Running simulation using NEURON for 30000.0 ms...\n",
      "  Done; run time = 230.77 s; real-time ratio: 0.13.\n",
      "\n",
      "Gathering data...\n",
      "  Done; gather time = 0.26 s.\n",
      "\n",
      "Analyzing...\n",
      "  Cells: 30\n",
      "  Connections: 900 (30.00 per cell)\n",
      "  Synaptic contacts: 1540 (51.33 per cell)\n",
      "  Spikes: 6452 (7.17 Hz)\n",
      "  Simulated time: 30.0 s; 1 workers\n",
      "  Run time: 230.77 s\n",
      "Saving output as patient_1_data.pkl ... \n",
      "Finished saving!\n",
      "  Done; saving time = 3.49 s.\n",
      "Plotting recorded cell traces ... cell\n",
      "Preparing spike data...\n",
      "Plotting raster...\n",
      "Unable to run analyzeRaster from sim.plotting and sim.analysis. Reason: module 'netpyne.analysis' has no attribute 'analyzeRaster'\n",
      "  Done; plotting time = 4.54 s\n",
      "\n",
      "Total time = 239.24 s\n",
      "Patient 1:\n",
      "Burst Duration starting at 5 seconds: 634 ms\n",
      "Burst Frequency: 4e-05 Hz\n",
      "\n",
      "Start time:  2024-03-15 13:42:39.010329\n",
      "\n",
      "Creating network of 3 cell populations on 1 hosts...\n",
      "  Number of cells on node 0: 30 \n",
      "  Done; cell creation time = 0.19 s.\n",
      "Making connections...\n",
      "  Number of connections on node 0: 870 \n",
      "  Number of synaptic contacts on node 0: 1450 \n",
      "  Done; cell connection time = 0.09 s.\n",
      "Adding stims...\n",
      "  Number of stims on node 0: 100 \n",
      "  Done; cell stims creation time = 0.02 s.\n",
      "Recording 60 traces of 2 types on node 0\n",
      "\n",
      "Running simulation using NEURON for 30000.0 ms...\n",
      "  Done; run time = 229.14 s; real-time ratio: 0.13.\n",
      "\n",
      "Gathering data...\n",
      "  Done; gather time = 0.29 s.\n",
      "\n",
      "Analyzing...\n",
      "  Cells: 30\n",
      "  Connections: 900 (30.00 per cell)\n",
      "  Synaptic contacts: 1540 (51.33 per cell)\n",
      "  Spikes: 6012 (6.68 Hz)\n",
      "  Simulated time: 30.0 s; 1 workers\n",
      "  Run time: 229.14 s\n",
      "Saving output as patient_2_data.pkl ... \n",
      "Finished saving!\n",
      "  Done; saving time = 3.46 s.\n",
      "Plotting recorded cell traces ... cell\n",
      "Preparing spike data...\n",
      "Plotting raster...\n",
      "Unable to run analyzeRaster from sim.plotting and sim.analysis. Reason: module 'netpyne.analysis' has no attribute 'analyzeRaster'\n",
      "  Done; plotting time = 4.46 s\n",
      "\n",
      "Total time = 237.69 s\n",
      "\n",
      "Start time:  2024-03-15 13:46:36.705698\n",
      "\n",
      "Creating network of 3 cell populations on 1 hosts...\n",
      "  Number of cells on node 0: 30 \n",
      "  Done; cell creation time = 0.00 s.\n",
      "Making connections...\n",
      "  Number of connections on node 0: 870 \n",
      "  Number of synaptic contacts on node 0: 1450 \n",
      "  Done; cell connection time = 0.09 s.\n",
      "Adding stims...\n",
      "  Number of stims on node 0: 100 \n",
      "  Done; cell stims creation time = 0.01 s.\n",
      "Recording 60 traces of 2 types on node 0\n",
      "\n",
      "Running simulation using NEURON for 30000.0 ms...\n",
      "  Done; run time = 231.94 s; real-time ratio: 0.13.\n",
      "\n",
      "Gathering data...\n",
      "  Done; gather time = 0.29 s.\n",
      "\n",
      "Analyzing...\n",
      "  Cells: 30\n",
      "  Connections: 900 (30.00 per cell)\n",
      "  Synaptic contacts: 1540 (51.33 per cell)\n",
      "  Spikes: 6012 (6.68 Hz)\n",
      "  Simulated time: 30.0 s; 1 workers\n",
      "  Run time: 231.94 s\n",
      "Saving output as patient_2_data.pkl ... \n",
      "Finished saving!\n",
      "  Done; saving time = 3.46 s.\n",
      "Plotting recorded cell traces ... cell\n",
      "Preparing spike data...\n",
      "Plotting raster...\n",
      "Unable to run analyzeRaster from sim.plotting and sim.analysis. Reason: module 'netpyne.analysis' has no attribute 'analyzeRaster'\n",
      "  Done; plotting time = 4.77 s\n",
      "\n",
      "Total time = 240.62 s\n",
      "Patient 2:\n",
      "Burst Duration starting at 5 seconds: 558 ms\n",
      "Burst Frequency: 4e-05 Hz\n",
      "\n",
      "Start time:  2024-03-15 13:50:37.335177\n",
      "\n",
      "Creating network of 3 cell populations on 1 hosts...\n",
      "  Number of cells on node 0: 30 \n",
      "  Done; cell creation time = 0.01 s.\n",
      "Making connections...\n",
      "  Number of connections on node 0: 870 \n",
      "  Number of synaptic contacts on node 0: 1450 \n",
      "  Done; cell connection time = 0.11 s.\n",
      "Adding stims...\n",
      "  Number of stims on node 0: 100 \n",
      "  Done; cell stims creation time = 0.01 s.\n",
      "Recording 60 traces of 2 types on node 0\n",
      "\n",
      "Running simulation using NEURON for 30000.0 ms...\n",
      "  Done; run time = 236.73 s; real-time ratio: 0.13.\n",
      "\n",
      "Gathering data...\n",
      "  Done; gather time = 0.50 s.\n",
      "\n",
      "Analyzing...\n",
      "  Cells: 30\n",
      "  Connections: 900 (30.00 per cell)\n",
      "  Synaptic contacts: 1540 (51.33 per cell)\n",
      "  Spikes: 5757 (6.40 Hz)\n",
      "  Simulated time: 30.0 s; 1 workers\n",
      "  Run time: 236.73 s\n",
      "Saving output as patient_3_data.pkl ... \n",
      "Finished saving!\n",
      "  Done; saving time = 4.42 s.\n",
      "Plotting recorded cell traces ... cell\n",
      "Preparing spike data...\n",
      "Plotting raster...\n",
      "Unable to run analyzeRaster from sim.plotting and sim.analysis. Reason: module 'netpyne.analysis' has no attribute 'analyzeRaster'\n",
      "  Done; plotting time = 4.82 s\n",
      "\n",
      "Total time = 246.66 s\n",
      "\n",
      "Start time:  2024-03-15 13:54:44.000504\n",
      "\n",
      "Creating network of 3 cell populations on 1 hosts...\n",
      "  Number of cells on node 0: 30 \n",
      "  Done; cell creation time = 0.01 s.\n",
      "Making connections...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  Number of connections on node 0: 870 \n",
      "  Number of synaptic contacts on node 0: 1450 \n",
      "  Done; cell connection time = 0.10 s.\n",
      "Adding stims...\n",
      "  Number of stims on node 0: 100 \n",
      "  Done; cell stims creation time = 0.01 s.\n",
      "Recording 60 traces of 2 types on node 0\n",
      "\n",
      "Running simulation using NEURON for 30000.0 ms...\n",
      "  Done; run time = 377.52 s; real-time ratio: 0.08.\n",
      "\n",
      "Gathering data...\n",
      "  Done; gather time = 0.38 s.\n",
      "\n",
      "Analyzing...\n",
      "  Cells: 30\n",
      "  Connections: 900 (30.00 per cell)\n",
      "  Synaptic contacts: 1540 (51.33 per cell)\n",
      "  Spikes: 5757 (6.40 Hz)\n",
      "  Simulated time: 30.0 s; 1 workers\n",
      "  Run time: 377.52 s\n",
      "Saving output as patient_3_data.pkl ... \n",
      "Finished saving!\n",
      "  Done; saving time = 4.87 s.\n",
      "Plotting recorded cell traces ... cell\n",
      "Preparing spike data...\n",
      "Plotting raster...\n"
     ]
    }
   ],
   "source": [
    "### Using DOPA effect as BL since this appears to be a normal healthy DA response SMS 15FEB2024\n",
    "import os\n",
    "from netpyne import specs, sim\n",
    "\n",
    "\n",
    "# Create a folder to save the results if it doesn't exist\n",
    "results_dir = 'results'\n",
    "if not os.path.exists(results_dir):\n",
    "    os.makedirs(results_dir)\n",
    "\n",
    "### Set variables for Individul Patient Loop\n",
    "\n",
    "## equations designed as % changes from baseline = 1\n",
    "\n",
    "## Variables below should increase at similar rates as a fnx of amph. mediated dopa release\n",
    "Nap_pyr_soma=[1.1, 1.1, 1.1, 1.1, 1.1, 1.1, 1.1];\n",
    "Hva_pyr_soma=[1.1, 1.1, 1.1, 1.1, 1.1, 1.1, 1.1];\n",
    "Iks_pyr_soma=[1.1, 1.1, 1.1, 1.1, 1.1, 1.1, 1.1];\n",
    "Nap_pyr_bas=[1.1, 1.1, 1.1, 1.1, 1.1, 1.1, 1.1];\n",
    "Hva_pyr_bas=[1.1, 1.1, 1.1, 1.1, 1.1, 1.1, 1.1];\n",
    "Iks_pyr_bas=[1.1, 1.1, 1.1, 1.1, 1.1, 1.1, 1.1];\n",
    "Nap_pyr_trunk=[1.1, 1.1, 1.1, 1.1, 1.1, 1.1, 1.1];\n",
    "Hva_pyr_trunk=[1.1, 1.1, 1.1, 1.1, 1.1, 1.1, 1.1];\n",
    "Iks_pyr_trunk=[1.1, 1.1, 1.1, 1.1, 1.1, 1.1, 1.1];\n",
    "Hva_pyr_apical=[1.1, 1.1, 1.1, 1.1, 1.1, 1.1, 1.1];\n",
    "Iks_pyr_apical=[1.1, 1.1, 1.1, 1.1, 1.1, 1.1, 1.1];\n",
    "\n",
    "#Connecticity of the circuit should change as a fnx of different patient (ie. schizo lower synatpic connectivity)\n",
    "exex=[1, 1, 1, 1, 1, 1];\n",
    "exin=[1, 1, 1, 1, 1, 1];\n",
    "inex=[1, 1, 1, 1, 1, 1];\n",
    "exexInterAssem=[1, 1, 1, 1, 1, 1];\n",
    "inin=[1, 1, 1, 1, 1, 1];\n",
    "\n",
    "#Strength of AMPA and NMDA conductions will be altered according to Goldman-Rakic 1998 publication\n",
    "AMPAgEE=[.85, .85, .85, .85, .85, .85, .85];\n",
    "NMDAgEE=[.85, .85, .85, .85, .85, .85, .85];\n",
    "AMPAgEE_assem=[.85, .85, .85, .85, .85, .85, .85];\n",
    "NMDAgEE_assem=[.85, .85, .85, .85, .85, .85, .85];\n",
    "AMPAgEI=[1.1, 1.2, 1.3, 1.3, 1.3, 1.3, 1.3];\n",
    "NMDAgEI=[1.1, 1.2, 1.3, 1.3, 1.3, 1.3, 1.3];\n",
    "GABAgIE=[1, 1, 1.1, 1.2, 1.3, 2, 3];\n",
    "GABAgII=[1, 1, 1.1, 1.2, 1.3, 2, 3];\n",
    "\n",
    "\n",
    "# Define the analyzeRaster function outside the loop\n",
    "def analyzeRaster(simData, patient):\n",
    "    spike_times = simData['spkt']\n",
    "    spike_cells = simData['spkid']\n",
    "    duration = 0\n",
    "    burst_count = 0\n",
    "    in_burst = False\n",
    "\n",
    "    for i in range(len(spike_times)):\n",
    "        if spike_times[i] >= 5000:  # Start analysis at 5 seconds\n",
    "            if not in_burst:\n",
    "                in_burst = True\n",
    "                burst_count += 1\n",
    "            duration += 1\n",
    "        else:\n",
    "            if in_burst:\n",
    "                in_burst = False\n",
    "\n",
    "    # Calculate burst frequency\n",
    "    burst_frequency = burst_count / (simConfig.duration - 5000)  # Duration in ms\n",
    "\n",
    "    # Save the data for each patient\n",
    "    output_dir = 'burstdata'\n",
    "    if not os.path.exists(output_dir):\n",
    "        os.makedirs(output_dir)\n",
    "    with open(os.path.join(output_dir, f'burstdata_patient_{patient}.txt'), 'w') as file:\n",
    "        file.write(f'Patient {patient}:\\n')\n",
    "        file.write(f'Burst Duration starting at 5 seconds: {duration} ms\\n')\n",
    "        file.write(f'Burst Frequency: {burst_frequency} Hz\\n')\n",
    "    \n",
    "    # Print the burst duration and frequency\n",
    "    print(f'Patient {patient}:')\n",
    "    print(f'Burst Duration starting at 5 seconds: {duration} ms')\n",
    "    print(f'Burst Frequency: {burst_frequency} Hz')\n",
    "    \n",
    "    # Return the burst duration and frequency\n",
    "    return duration, burst_frequency\n",
    "\n",
    "\n",
    "for patient in (0,1,2,3,4,5):\n",
    "\n",
    "        Dt = 0.05\n",
    "        Tfinal = 30000\n",
    "\n",
    "        # Network parameters\n",
    "        netParams = specs.NetParams()  # object of class NetParams to store the network parameters\n",
    "        netParams.defaultThreshold = 0.0\n",
    "\n",
    "        ###############################################################################\n",
    "        ## Cell types\n",
    "        ###############################################################################\n",
    "        # Pyramidal cell\n",
    "        PyrCell = {'secs':{}, 'globals': {'ko0_k_ion': 3.82, 'ki0_k_ion': 160, 'cao0_ca_ion': 2, 'cai0_ca_ion': 50e-6}}\n",
    "        # soma\n",
    "        PyrCell['secs']['soma'] = {'geom': {}, 'mechs': {}}\n",
    "        PyrCell['secs']['soma']['geom'] = {'diam': 23, 'L': 23, 'cm': 1.2, 'Ra': 150, 'nseg': 1}\n",
    "        PyrCell['secs']['soma']['mechs'] = {\n",
    "            'pas': {'g': 3.33e-5, 'e': -70},\n",
    "            'kdyn': {},\n",
    "            'Naf': {'gnafbar': 0.086},#same\n",
    "            'Nap': {'gnapbar': [0.00154*Nap_pyr_soma[patient]]},#reduced to 75% of baseline #DopaEffect\n",
    "            'Hva': {'ghvabar': [0.000272*Hva_pyr_soma[patient]]},#reduced to 80% of baseline (as in original model) #DopaEffect\n",
    "            'kdr': {'gkdrbar': 0.0338},#same\n",
    "            'IKs': {'gKsbar': [0.000105*Iks_pyr_soma[patient]]},# reducted by 75% of baseline (50% of baseline used in original mode) #DopaEffect\n",
    "            'iC': {'gkcbar': 0.0022},#same\n",
    "            'cadyn': {'CAF': 385.948e-9, 'tca': 250}}#same\n",
    "\n",
    "        # basal dendrite\n",
    "        PyrCell['secs']['Bdend'] = {'geom': {}, 'mechs': {}}\n",
    "        PyrCell['secs']['Bdend']['geom'] = {'diam': 16.0, 'L': 150.0, 'cm': 2.3, 'Ra': 150, 'nseg': 1}\n",
    "        PyrCell['secs']['Bdend']['topol'] = {'parentSec': 'soma', 'parentX': 0, 'childX': 0}\n",
    "        PyrCell['secs']['Bdend']['mechs'] = {\n",
    "            'pas': {'g': 6.39e-5, 'e': -70},\n",
    "            'kdyn': {},\n",
    "            'Naf': {'gnafbar': 0.028},\n",
    "            'Nap': {'gnapbar': [0.0007*Nap_pyr_bas[patient]]},\n",
    "            'Hva': {'ghvabar': [0.00056*Hva_pyr_bas[patient]]},\n",
    "            'kdr': {'gkdrbar': 0.0092},\n",
    "            'IKs': {'gKsbar': [0.00018*Iks_pyr_bas[patient]]},\n",
    "            'iC': {'gkcbar': 0.0038},\n",
    "            'cadyn': {'CAF': 964.87e-9, 'tca': 120}}\n",
    "\n",
    "        # apical proximal dendrite\n",
    "        PyrCell['secs']['Adend_proximal'] = {'geom': {}, 'mechs': {}}\n",
    "        PyrCell['secs']['Adend_proximal']['geom'] = {'diam': 2.6, 'L': 400.0, 'cm': 2.3, 'Ra': 150, 'nseg': 1}\n",
    "        PyrCell['secs']['Adend_proximal']['topol'] = {'parentSec': 'soma', 'parentX': 1, 'childX': 0}\n",
    "        PyrCell['secs']['Adend_proximal']['mechs'] = {\n",
    "            'pas': {'g': 6.39e-5, 'e': -70},\n",
    "            'kdyn': {},\n",
    "            'Naf': {'gnafbar': 0.028},\n",
    "            'Nap': {'gnapbar': [0.0007*Nap_pyr_trunk[patient]]},\n",
    "            'Hva': {'ghvabar': [0.00056*Hva_pyr_trunk[patient]]},\n",
    "            'kdr': {'gkdrbar': 0.0092},\n",
    "            'IKs': {'gKsbar': [0.00018*Iks_pyr_trunk[patient]]},\n",
    "            'iC': {'gkcbar': 0.0038},\n",
    "            'cadyn': {'CAF': 964.87e-9, 'tca': 120}}\n",
    "\n",
    "        # apical distal dendrite\n",
    "        PyrCell['secs']['Adend_distal'] = {'geom': {}, 'mechs': {}}\n",
    "        PyrCell['secs']['Adend_distal']['geom'] = {'diam': 2.6, 'L': 400.0, 'cm': 2.3, 'Ra': 150, 'nseg': 1}\n",
    "        PyrCell['secs']['Adend_distal']['topol'] = {'parentSec': 'Adend_proximal', 'parentX': 1, 'childX': 0}\n",
    "        PyrCell['secs']['Adend_distal']['mechs'] = {\n",
    "            'pas': {'g': 6.39e-5, 'e': -70},\n",
    "            'kdyn': {},\n",
    "            'Naf': {'gnafbar': 0.028},\n",
    "            'Hva': {'ghvabar': [0.00017*Hva_pyr_apical[patient]]},\n",
    "            'kdr': {'gkdrbar': 0.0092},\n",
    "            'IKs': {'gKsbar': [0.00018*Iks_pyr_apical[patient]]},\n",
    "            'iC': {'gkcbar': 0.0022},\n",
    "            'cadyn': {'CAF': 964.87e-9, 'tca': 80}}\n",
    "\n",
    "        netParams.cellParams['PyrCell'] = PyrCell\n",
    "\n",
    "        ###############################################################################\n",
    "        # Inhibitory cell\n",
    "        InhCell = {'secs':{}}\n",
    "        # soma\n",
    "        InhCell['secs']['soma'] = {'geom': {}, 'mechs': {}}\n",
    "        InhCell['secs']['soma']['geom'] = {'diam': 15, 'L': 15, 'cm': 1.0, 'Ra': 150, 'nseg': 1}\n",
    "        InhCell['secs']['soma']['mechs'] = {\n",
    "            'pas': {'g': 1.0e-5, 'e': -68},\n",
    "            'kdyn': {},\n",
    "            'Naf': {'gnafbar': 0.100},\n",
    "            'kdr': {'gkdrbar': 0.040}}\n",
    "\n",
    "        # dendrite\n",
    "        InhCell['secs']['dend'] = {'geom': {}, 'mechs': {}}\n",
    "        InhCell['secs']['dend']['geom'] = {'diam': 10, 'L': 150, 'cm': 1.92, 'Ra': 150, 'nseg': 1}\n",
    "        InhCell['secs']['dend']['topol'] = {'parentSec': 'soma', 'parentX': 1, 'childX': 0}\n",
    "        InhCell['secs']['dend']['mechs'] = {\n",
    "            'pas': {'g': 1.92e-5, 'e': -68},\n",
    "            'kdyn': {},\n",
    "            'Naf': {'gnafbar': 0.020},\n",
    "            'kdr': {'gkdrbar': 0.008}}\n",
    "\n",
    "        netParams.cellParams['InhCell'] = InhCell\n",
    "\n",
    "\n",
    "        ###############################################################################\n",
    "        ## Synaptic mechs\n",
    "        ###############################################################################\n",
    "\n",
    "        netParams.synMechParams['AMPA'] = {'mod': 'Exp2Syn', 'tau1': 0.5, 'tau2': 2.5, 'e': 0}\n",
    "        netParams.synMechParams['NMDA'] = {'mod': 'Exp2SynNMDA', 'tau1': 10, 'tau2': 250.0, 'e': 0}\n",
    "        netParams.synMechParams['GABA'] = {'mod': 'ExpSyn', 'tau': 1.0, 'e': -75}\n",
    "\n",
    "\n",
    "        ###############################################################################\n",
    "        ## Populations\n",
    "        ###############################################################################\n",
    "        # Structuring excitatory assemblies\n",
    "        Nassemblies = 2\n",
    "        Nexc = 10\n",
    "        Nexc_total = Nassemblies * Nexc       # without overlap between assemblies\n",
    "        for na in range(Nassemblies):\n",
    "            ### Populations\n",
    "            netParams.popParams['Pyr_Assembly'+str(na+1)] = {'numCells': Nexc, 'cellType': 'PyrCell'}\n",
    "\n",
    "        # Inhibitory\n",
    "        Ninh = 10                             # global inhibitory neurons \n",
    "        netParams.popParams['InhPop'] = {'numCells': Ninh, 'cellType': 'InhCell'}\n",
    "\n",
    "\n",
    "        ###############################################################################\n",
    "        ## Connections\n",
    "        ###############################################################################\n",
    "        for na in range(Nassemblies):\n",
    "            ## Intra-assembly\n",
    "            # pyr-to-pyr AMPA/NMDA, to basal dendrite\n",
    "            netParams.connParams['PYR->PYR(intra)_AMPA/NMDA_Assembly'+str(na+1)] = {\n",
    "                'preConds': {'pop': 'Pyr_Assembly'+str(na+1)}, \n",
    "                'postConds': {'pop': 'Pyr_Assembly'+str(na+1)},\n",
    "                'connFunc': 'fullConn',\n",
    "                'weight': [2.0e-3*AMPAgEE[patient], 0.56e-3*NMDAgEE[patient]],\n",
    "                'delay': 'uniform(2,4)',\n",
    "                'sec': ['Bdend','Adend_proximal'],   # synapses sampled at these regions\n",
    "                'loc': 0.5,\n",
    "                'synMech': ['AMPA','NMDA']}\n",
    "\n",
    "            # pyr-to-inh AMPA/NMDA, to dendrite\n",
    "            netParams.connParams['PYR->INH_AMPA/NMDA_Assembly'+str(na+1)] = {\n",
    "                'preConds': {'pop': 'Pyr_Assembly'+str(na+1)}, \n",
    "                'postConds': {'pop': 'InhPop'},\n",
    "                'connFunc': 'fullConn',\n",
    "                'weight': [0.6*value for value in [2.0e-3*AMPAgEI[patient], 0.56e-3*NMDAgEI[patient]]],\n",
    "                'delay': 'uniform(2,4)',\n",
    "                'sec': 'dend',\n",
    "                'loc': 0.5,\n",
    "                'synMech': ['AMPA','NMDA']}\n",
    "            \n",
    "            # inh-to-pyr GABA, to soma\n",
    "            netParams.connParams['INH->PYR_GABA_Assembly'+str(na+1)] = {\n",
    "                'preConds': {'pop': 'InhPop'},\n",
    "                'postConds': {'pop': 'Pyr_Assembly'+str(na+1)}, \n",
    "                'connFunc': 'fullConn',\n",
    "                'weight': 7.8e-3*GABAgIE[patient],\n",
    "                'delay': 'uniform(2,4)',\n",
    "                'sec': 'soma',\n",
    "                'loc': 0.5,\n",
    "                'synMech': 'GABA'}\n",
    "            \n",
    "            ## Inter-assembly\n",
    "            for nb in range(Nassemblies):\n",
    "                if nb != na:\n",
    "                    # pyr-to-pyr AMPA/NMDA, to basal dendrite\n",
    "                    netParams.connParams['PYR->PYR(inter)_AMPA/NMDA_Assembly'+str(na+1)+'_to_Assembly'+str(nb+1)] = {\n",
    "                        'preConds': {'pop': 'Pyr_Assembly'+str(na+1)}, \n",
    "                        'postConds': {'pop': 'Pyr_Assembly'+str(nb+1)},\n",
    "                        'connFunc': 'fullConn',\n",
    "                        'weight': [0.1*value for value in [2.0e-3*AMPAgEE_assem[patient], 0.56e-3*NMDAgEE_assem[patient]]],\n",
    "                        'delay': 'uniform(2,4)',\n",
    "                        'sec': ['Bdend','Adend_proximal'],   # synapses sampled at these regions\n",
    "                        'loc': 0.5,\n",
    "                        'synMech': ['AMPA','NMDA']}\n",
    "                \n",
    "        # inh-to-inh GABA, to soma\n",
    "        netParams.connParams['INH->INH_GABA'] = {\n",
    "            'preConds': {'pop': 'InhPop'},\n",
    "            'postConds': {'pop': 'InhPop'}, \n",
    "            'connFunc': 'fullConn',\n",
    "            'weight': 0.65e-3*GABAgII[patient],\n",
    "            'delay': 'uniform(2,4)',\n",
    "            'sec': 'soma',\n",
    "            'loc': 0.5,\n",
    "            'synMech': 'GABA'}\n",
    "                    \n",
    "\n",
    "        ###############################################################################\n",
    "        ## Background stimulation\n",
    "        ###############################################################################\n",
    "        ### Source Stimulation - Spontaneous firing\n",
    "        netParams.stimSourceParams['NetStim_PYR_exc']  = {'type': 'NetStim', 'interval': 1.0, 'start': 0, 'noise': 1, 'seed': 1}\n",
    "        netParams.stimSourceParams['NetStim_PYR_exc2'] = {'type': 'NetStim', 'interval': 5.0, 'start': 0, 'noise': 1, 'seed': 2}\n",
    "        netParams.stimSourceParams['NetStim_PYR_inh']  = {'type': 'NetStim', 'interval': 1.0, 'start': 0, 'noise': 1, 'seed': 3}\n",
    "        netParams.stimSourceParams['NetStim_INH_exc']  = {'type': 'NetStim', 'interval': 1.0, 'start': 0, 'noise': 1, 'seed': 4}\n",
    "        netParams.stimSourceParams['NetStim_INH_exc2'] = {'type': 'NetStim', 'interval': 5.0, 'start': 0, 'noise': 1, 'seed': 5}\n",
    "        netParams.stimSourceParams['NetStim_INH_inh']  = {'type': 'NetStim', 'interval': 1.0, 'start': 0, 'noise': 1, 'seed': 6}\n",
    "\n",
    "        ### Target Stimulation - Spontaneous firing\n",
    "        for na in range(Nassemblies):\n",
    "            # exc to PYR\n",
    "            netParams.stimTargetParams['NetStim_PYR_exc->PYR_Assembly'+str(na+1)] = {\n",
    "                'source': 'NetStim_PYR_exc',\n",
    "                'conds': {'pop': 'Pyr_Assembly'+str(na+1)},\n",
    "                'sec': 'soma',\n",
    "                'loc': 0.5,\n",
    "                'weight': 1.0e-3,\n",
    "                'delay': 2*Dt,\n",
    "                'synMech': 'AMPA'}\n",
    "\n",
    "            netParams.stimTargetParams['NetStim_PYR_exc->PYR_Assembly_NMDA'+str(na+1)] = {\n",
    "                'source': 'NetStim_PYR_exc2',\n",
    "                'conds': {'pop': 'Pyr_Assembly'+str(na+1)},\n",
    "                'sec': 'soma',\n",
    "                'loc': 0.5,\n",
    "                'weight': 0.1e-3,\n",
    "                'delay': 2*Dt,\n",
    "                'synMech': 'NMDA'}\n",
    "\n",
    "            # inh to PYR\n",
    "            netParams.stimTargetParams['NetStim_PYR_inh->PYR_Assembly'+str(na+1)] = {\n",
    "                'source': 'NetStim_PYR_inh',\n",
    "                'conds': {'pop': 'Pyr_Assembly'+str(na+1)},\n",
    "                'sec': 'soma',\n",
    "                'loc': 0.5,\n",
    "                'weight': 7.5e-3,\n",
    "                'delay': 2*Dt,\n",
    "                'synMech': 'GABA'}\n",
    "\n",
    "        # exc to INH\n",
    "        netParams.stimTargetParams['NetStim_INH_exc->INH'] = {\n",
    "            'source': 'NetStim_INH_exc',\n",
    "            'conds': {'pop': 'InhPop'},\n",
    "            'sec': 'soma',\n",
    "            'loc': 0.5,\n",
    "            'weight': 0.55e-3,\n",
    "            'delay': 2*Dt,\n",
    "            'synMech': 'AMPA'}\n",
    "\n",
    "        netParams.stimTargetParams['NetStim_INH_exc->INH_NMDA'] = {\n",
    "            'source': 'NetStim_INH_exc2',\n",
    "            'conds': {'pop': 'InhPop'},\n",
    "            'sec': 'soma',\n",
    "            'loc': 0.5,\n",
    "            'weight': 0.055e-3,\n",
    "            'delay': 2*Dt,\n",
    "            'synMech': 'NMDA'}\n",
    "\n",
    "        # inh to INH\n",
    "        netParams.stimTargetParams['NetStim_INH_inh->INH'] = {\n",
    "            'source': 'NetStim_INH_inh',\n",
    "            'conds': {'pop': 'InhPop'},\n",
    "            'sec': 'soma',\n",
    "            'loc': 0.5,\n",
    "            'weight': 1.25e-3,\n",
    "            'delay': 2*Dt,\n",
    "            'synMech': 'GABA'}\n",
    "\n",
    "\n",
    "        ###############################################################################\n",
    "        ## Assembly recall\n",
    "        ###############################################################################\n",
    "        netParams.stimSourceParams['IClamp'] =  {'type': 'IClamp', 'del': Tfinal/6, 'dur': 100, 'amp': 1.00}\n",
    "        netParams.stimTargetParams['IClamp->Pyr_Assembly1'] = {\n",
    "                'source': 'IClamp',\n",
    "                'sec': 'soma',\n",
    "                'loc': 0.5,\n",
    "                'conds': {'pop':'Pyr_Assembly1'}}\n",
    "\n",
    "\n",
    "    \n",
    "        \n",
    "        ###############################################################################\n",
    "        # Simulation options\n",
    "        ###############################################################################\n",
    "        simConfig = specs.SimConfig()\n",
    "\n",
    "        simConfig.duration = Tfinal\n",
    "        simConfig.dt = Dt\n",
    "        simConfig.hParams = {'v_init': -66,'celsius': 23}\n",
    "        simConfig.verbose = False\n",
    "        simConfig.recordTraces = {'V_soma':{'sec':'soma','loc':0.5,'var':'v'}}\n",
    "        #                          'ko_soma': {'sec': 'soma', 'loc': 0.5, 'var': 'ko'},\n",
    "        #                          'cai_basal':{'sec':'Bdend','loc':0.5,'var':'cai'},\n",
    "        #                          'ko_apical': {'sec': 'Adend_distal', 'loc': 0.5, 'var': 'ko'},\n",
    "        #                          'synapse_inmda':{'sec':'soma','loc':0.5,'synMech':'NMDA','var':'i'},}  # Dict with traces to record\n",
    "        simConfig.recordStim = True\n",
    "        simConfig.recordStep = 0.1        # Step size in ms to save data (eg. V traces, LFP, etc)\n",
    "        simConfig.filename = '00'         # Set file output name\n",
    "        #simConfig.saveJson = True\n",
    "\n",
    "        #simConfig.analysis['plotTraces'] = {'include': ['all'], 'saveFig': True}\n",
    "        #simConfig.analysis['plotRaster'] = {'saveFig': True}\n",
    "        \n",
    "        # Update simConfig to save figures with patient flags\n",
    "        simConfig.analysis['plotTraces'] = {'include': ['all'], 'saveFig': 'traces_pat{0}.png'.format(patient)}\n",
    "        simConfig.analysis['plotRaster'] = {'saveFig': 'raster_pat{0}.png'.format(patient)}\n",
    "       \n",
    "\n",
    "\n",
    "        # Add this block of code inside the loop\n",
    "        simConfig.recordTraces = {'V_soma':{'sec':'soma','loc':0.5,'var':'v'}}\n",
    "        simConfig.recordStim = True\n",
    "        simConfig.recordStep = 0.1\n",
    "        simConfig.recordSpikes = 'all'  # Record spike times for all cells\n",
    "\n",
    "        def analyzeRaster(simData, patient):\n",
    "            spike_times = simData['spkt']\n",
    "            spike_cells = simData['spkid']\n",
    "            duration = 0\n",
    "            burst_count = 0\n",
    "            in_burst = False\n",
    "\n",
    "            for i in range(len(spike_times)):\n",
    "                if spike_times[i] >= 5000:  # Start analysis at 5 seconds\n",
    "                    if spike_cells[i] < Nassemblies * Nexc:  # Check if spike is from assembly 1 pyramidal neurons\n",
    "                        if not in_burst:\n",
    "                            in_burst = True\n",
    "                            burst_count += 1\n",
    "                        duration += 1\n",
    "                else:\n",
    "                    if in_burst:\n",
    "                        in_burst = False\n",
    "\n",
    "            # Calculate burst frequency\n",
    "            burst_frequency = burst_count / (simConfig.duration - 5000)  # Duration in ms\n",
    "\n",
    "            # Save the data for each patient\n",
    "            output_dir = 'burstdata'\n",
    "            if not os.path.exists(output_dir):\n",
    "                os.makedirs(output_dir)\n",
    "            with open(os.path.join(output_dir, f'burstdata_patient_{patient}.txt'), 'w') as file:\n",
    "                file.write(f'Patient {patient}:\\n')\n",
    "                file.write(f'Burst Duration starting at 5 seconds: {duration} ms\\n')\n",
    "                file.write(f'Burst Frequency: {burst_frequency} Hz\\n')\n",
    "\n",
    "            # Print the burst duration and frequency\n",
    "            print(f'Patient {patient}:')\n",
    "            print(f'Burst Duration starting at 5 seconds: {duration} ms')\n",
    "            print(f'Burst Frequency: {burst_frequency} Hz')\n",
    "\n",
    "            # Return the burst duration and frequency\n",
    "            return duration, burst_frequency\n",
    "\n",
    "        \n",
    "        simConfig.analysis['analyzeRaster'] = {'func': analyzeRaster, 'kwargs': {'patient': patient}}\n",
    "        simConfig.filename = 'patient_{0}'.format(patient)\n",
    "        simConfig.savePickle = True  # Save data in pickle format\n",
    "\n",
    "        sim.createSimulateAnalyze(netParams=netParams, simConfig=simConfig)\n",
    "\n",
    "        \n",
    "        \n",
    "        # Create network and run simulation\n",
    "        sim.createSimulateAnalyze(netParams = netParams, simConfig = simConfig)\n",
    "        #sim.create(netParams = netParams, simConfig = simConfig)\n",
    "        \n",
    "        # Call analyzeRaster function to calculate burst duration and frequency\n",
    "        analyzeRaster(sim.simData, patient)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73c3f454-c3c0-436c-873f-8d2b82d0e465",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}