Construct_rho_alpha_extraction_all_morphs.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-07-04T00:20:39.136134Z",
     "start_time": "2022-07-04T00:20:39.107938Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>.container { width:100% !important; }</style>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using matplotlib backend: Qt5Agg\n",
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "# 11.9 22:17\n",
    "from neuron import gui,h\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import pandas as pd\n",
    "from tqdm.auto import tqdm\n",
    "from IPython.core.display import display, HTML\n",
    "import os\n",
    "h.load_file(\"nrngui.hoc\")\n",
    "h.load_file(\"import3d.hoc\")\n",
    "display(HTML(\"<style>.container { width:100% !important; }</style>\"))\n",
    "\n",
    "%pylab\n",
    "%matplotlib inline\n",
    "cells = {'JM072303':{'Rm':1000},'Hay_2011_passive':{'Rm':15000}, 'NMO_130658':{'Rm':3000}, 'NMO_161366':{'Rm':3500}}\n",
    "\n",
    "cell_name = 'NMO_161366' # decide which cell to process\n",
    "#NMO_161366  - pyramidal - http://neuromorpho.org/neuron_info.jsp?neuron_name=20180102_sample_1 - https://www.nature.com/articles/s41586-020-2907-3\n",
    "#NMO_130658 - basket - http://neuromorpho.org/neuron_info.jsp?neuron_name=2017-16-05-slice-1-cell-6-rotated -  https://pubmed.ncbi.nlm.nih.gov/31519874/\n",
    "#JM072303 - spiny - http://neuromorpho.org/neuron_info.jsp?neuron_name=JM072303 - https://pubmed.ncbi.nlm.nih.gov/16337918/\n",
    "\n",
    "\n",
    "syn_conds = [0.001] # synaptic conductance to use\n",
    "nstim_start = 40\n",
    "cur_amp = 0.1 # nA\n",
    "\n",
    "\n",
    "def create_cell():\n",
    "    \n",
    "    h.load_file(\"models/Hay_2011_passive/L5PCtemplate.hoc\")\n",
    "    h.load_file(\"models/Hay_2011_passive/L5PCtemplate_swc.hoc\")\n",
    "    \n",
    "    if cell_name == 'Hay_2011_passive':\n",
    "        cell = h.L5PCtemplate('models/Hay_2011_passive/cell1.asc')\n",
    "    elif cell_name == 'JM072303':\n",
    "        cell = h.L5PCtemplate('models/JM072303/JM072303.asc')\n",
    "    elif cell_name =='NMO_130658':\n",
    "        cell = h.L5PCtemplateswc('models/NMO_130658/2017-13-02-slice-2-cell-3-rotated.SWC')\n",
    "    elif cell_name =='NMO_161366':\n",
    "        cell = h.L5PCtemplateswc('models/NMO_161366/20180102_sample_1.SWC')\n",
    "    h.area(0.5,sec=cell.soma[0]) # make sure diam reflects 3d points\n",
    "\n",
    "    \n",
    "    soma = cell.soma[0]\n",
    "    for sec in h.allsec():\n",
    "        sec.insert(\"pas\")\n",
    "        sec.g_pas = 1/cells[cell_name]['Rm']\n",
    "        sec.cm =1 \n",
    "        sec.Ra = 100\n",
    "    \n",
    "    h(\"freq = 100\")#      // Hz, frequency at which AC length constant will be computedf\n",
    "    h(\"d_lambda = 0.1\")\n",
    "    h(\"forall { nseg = int((L/(d_lambda*lambda_f(freq))+0.9)/2)*2 + 1  }\")\n",
    "\n",
    "    vvec = h.Vector().record(soma(0.5)._ref_v)\n",
    "    tvec = h.Vector().record(h._ref_t)\n",
    "\n",
    "    iclamp = h.IClamp(soma(0.5))\n",
    "    seclamp = h.SEClamp(soma(0.5))\n",
    "    seclamp_current = h.Vector().record(seclamp._ref_i)\n",
    "\n",
    "    h.v_init = soma.e_pas\n",
    "\n",
    "\n",
    "    #check cell input resistance\n",
    "\n",
    "    imp = h.Impedance()\n",
    "    imp.loc(0.5,sec=soma)\n",
    "    imp.compute(0,1)\n",
    "    Rin = imp.input(0.5,sec=soma)\n",
    "    return cell, vvec, tvec, iclamp, seclamp, Rin\n",
    "\n",
    "\n",
    "def create_synapse(seg, e):\n",
    "    nstim = h.NetStim()\n",
    "    # set other variables\n",
    "    nstim.number = 1\n",
    "    nstim.noise = 0\n",
    "    nstim.start = nstim_start\n",
    "\n",
    "    \n",
    "    dend_syn = h.Exp2Syn(seg)\n",
    "    dend_syn.tau1=0.0001#0.3\n",
    "    dend_syn.tau2=9e9#10\n",
    "    dend_syn.e = e\n",
    "    \n",
    "    dend_ncon = h.NetCon(nstim,dend_syn)\n",
    "    dend_ncon.weight[0] = 0\n",
    "    \n",
    "    return nstim, dend_syn, dend_ncon\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "def get_cell_specs():\n",
    "    cell, vvec, tvec, iclamp, seclamp, Rin = create_cell()\n",
    "    all_segs = []\n",
    "    for sec in h.allsec():\n",
    "        for seg in sec:\n",
    "            all_segs.append(seg)\n",
    "    return Rin, len(all_segs), cell.soma[0].e_pas, cell, vvec, tvec, iclamp, seclamp, Rin\n",
    "\n",
    "    \n",
    "def calc_syn_effects(x):\n",
    "    '''\n",
    "    x - list of inds for segs to calculate\n",
    "    '''\n",
    "    st = x[0] # do seg from st to en\n",
    "    en = x[1]\n",
    "    cell, vvec, tvec, iclamp, seclamp, Rin = create_cell()\n",
    "\n",
    "\n",
    "    h.tstop = nstim_start+100 + 300\n",
    "    iclamp.delay  = nstim_start+200\n",
    "    iclamp.dur = h.tstop + 10\n",
    "    \n",
    "    iclamp.amp = cur_amp #nA\n",
    "\n",
    "\n",
    "    h.dt=0.2\n",
    "    h.steps_per_ms = 1/h.dt\n",
    "\n",
    "\n",
    "    to_save = {'syn_dv':[], 'iclamp_dv':[], 'seg_dist':[],'seg_length':[],'voltage':[],'synapse_reversal':[],'seg_name':[],'syn_g':[]}\n",
    "    #syn_dv - dV of synapse\n",
    "    #iclamp_dv - dV of iclamp, start after synapses\n",
    "    #syn_g - nanoS\n",
    "\n",
    "    all_segs = []\n",
    "    for sec in h.allsec():\n",
    "        for seg in sec:\n",
    "            all_segs.append(seg)\n",
    "    \n",
    "    if st>len(all_segs):\n",
    "        return []\n",
    "    \n",
    "    all_segs = all_segs[st:en]\n",
    "    \n",
    "    for seg in all_segs:\n",
    "        for syn_cond in syn_conds:\n",
    "            for syn_e in [0,-80]:\n",
    "                nstim, dend_syn, dend_ncon  = create_synapse(seg, syn_e)\n",
    "                dend_ncon.weight[0] =syn_cond\n",
    "\n",
    "                syn_g = dend_ncon.weight[0]*10**3 #from micro to nano\n",
    "\n",
    "                h.run()\n",
    "                #plot(tvec,vvec)\n",
    "                iclamp_dv = vvec[-1]-vvec[int((iclamp.delay-10)/h.dt)]\n",
    "                syn_dv  = vvec[int((iclamp.delay-10)/h.dt)]-vvec[10]\n",
    "\n",
    "                to_save['syn_dv'].append(syn_dv)\n",
    "                to_save['iclamp_dv'].append(iclamp_dv)\n",
    "                to_save['seg_dist'].append(h.distance(cell.soma[0](0.5),seg))\n",
    "                to_save['seg_length'].append(seg.sec.L/seg.sec.nseg)\n",
    "                to_save['synapse_reversal'].append(syn_e)\n",
    "                to_save['voltage'].append(list(vvec))\n",
    "                to_save['syn_g'].append(syn_g)\n",
    "                to_save['seg_name'].append(str(seg))\n",
    "                nstim.number = 0\n",
    "                dend_ncon.weight[0] = 0\n",
    "                del nstim,dend_ncon,dend_syn\n",
    "\n",
    "    return to_save"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calc voltage deflaction per synapse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-07-04T00:21:30.759067Z",
     "start_time": "2022-07-04T00:20:41.143212Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "if number of segments is larger than 3000, this will be a probalem\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7bd2d36cbbb24aa784874d5ea14fbc47",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/250 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## test\n",
    "# data = calc_syn_effects([10,20])\n",
    "# pd.DataFrame(data)\n",
    " \n",
    "    \n",
    "chunk_size = 10\n",
    "max_send = 2500\n",
    "to_send =[[i*10,(i+1)*10] for i in range(300) if i*10<max_send]\n",
    "print(\"if number of segments is larger than 3000, this will be a probalem\")\n",
    "\n",
    "from multiprocessing import Pool\n",
    "\n",
    "outputs = []\n",
    "pbar= tqdm(total=len(to_send))\n",
    "with Pool(processes=5) as pool:\n",
    "    for out in pool.imap_unordered(calc_syn_effects,to_send):\n",
    "        pbar.update(1)\n",
    "        if len(out)>0:\n",
    "            outputs.append(out)\n",
    "\n",
    "#merge output\n",
    "full_dict = {k:[] for k in outputs[0]}\n",
    "for out in outputs:\n",
    "    for k in out:\n",
    "        full_dict[k].extend(out[k])\n",
    "df = pd.DataFrame(full_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-07-04T00:28:17.749755Z",
     "start_time": "2022-07-04T00:28:17.337278Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "62.23537626767067\n"
     ]
    }
   ],
   "source": [
    "Rin, segs_num, e_pas, cell, vvec, tvec, iclamp, seclamp, Rin = get_cell_specs()\n",
    "print(Rin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calc Rho and Alpha"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-07-04T00:39:38.629998Z",
     "start_time": "2022-07-04T00:39:37.570779Z"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "df['g_eff'] = None\n",
    "df['alpha'] = None\n",
    "df['rho'] = None\n",
    "df['alpharho'] = None\n",
    "\n",
    "plotting = False\n",
    "_,axes = subplots(ncols=4, sharex=True,figsize=(10,3)) if plotting else _,None\n",
    "for jj,row in df.iterrows():\n",
    "    \n",
    "    g_eff = 1/(row['iclamp_dv']/cur_amp)*10**3 #mv/nA ?=uS --*10**3--> nS ?\n",
    "    df.loc[row.name,'g_eff'] = g_eff\n",
    "    \n",
    "    if plotting:\n",
    "        \n",
    "        sca(axes[0])\n",
    "        plot(row['syn_g'],g_eff,'-o')\n",
    "        xlabel('gs (nS)')\n",
    "        ylabel('g_eff (nS)')\n",
    "       \n",
    "\n",
    "    ## calc ar ###\n",
    "    ##############\n",
    "    gl_S = 1/Rin*10**-6 # uS -> S\n",
    "    cur_amp_A = cur_amp*10**-9 # nA\n",
    "    dvs_V = row['iclamp_dv']*10**-3 # mV ->V\n",
    "    g_syns_S = row['syn_g']*10**-9 #  nS ->S    \n",
    "    ar = ((cur_amp_A/dvs_V)-gl_S)/g_syns_S\n",
    "    df.loc[row.name,'alpharho'] = ar\n",
    "    \n",
    "    if plotting:\n",
    "        sca(axes[1])\n",
    "        plot(row['syn_g'],ar ,'-o')\n",
    "        xlabel('gs (nS)')\n",
    "        ylabel(r'$\\alpha \\rho$')\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "    ## calc alpha ##\n",
    "    ################\n",
    "    dvs_fromsyn_V = row['syn_dv']*10**-3 # mV ->V\n",
    "    g_eff_S = g_eff*10**-9\n",
    "    V_diff = (e_pas-row['synapse_reversal'])*10**-3 #mV ->V\n",
    "    g_syns_S\n",
    "    a = -1*(dvs_fromsyn_V*g_eff_S)/(V_diff*g_syns_S)\n",
    "    df.loc[row.name,'alpha'] = a\n",
    "    \n",
    "    if plotting:\n",
    "        sca(axes[2])\n",
    "        tm = ar/a\n",
    "        plot(row['syn_g'],a ,'-o')\n",
    "        xlabel('gs (nS)')\n",
    "        ylabel(r'$\\alpha$')\n",
    "        tight_layout()\n",
    "        #ylim(min(list(tm) + list(a))-0.01,max(list(tm) + list(a))+0.01) \n",
    "    ## calc rho   ##\n",
    "    ################\n",
    "    df.loc[row.name,'rho'] = ar/a\n",
    "    if plotting:\n",
    "        sca(axes[3])\n",
    "        plot(row['syn_g'],ar/a ,'-o')\n",
    "        xlabel('gs (nS)')\n",
    "        ylabel(r'$\\rho$')\n",
    "        #ylim(min(list(tm) + list(a))-0.01,max(list(tm) + list(a))+0.01) \n",
    "        tight_layout()\n",
    "        #suptitle(f'dist from soma = {row[\"seg_x\"].iloc[0]:9.5} um | loc = {row[\"seg_obj\"].iloc[0]}')\n",
    "        tight_layout(rect=[0, 0.03, 1, 0.95])\n",
    "        \n",
    "    if jj>4:\n",
    "        plotting = False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-07-04T00:28:02.392310Z",
     "start_time": "2022-07-04T00:28:02.389015Z"
    }
   },
   "source": [
    "# plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-07-04T00:39:40.514908Z",
     "start_time": "2022-07-04T00:39:39.916853Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0.98, 'NMO_161366')"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA70AAAEmCAYAAACu4/cFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABjQUlEQVR4nO3df5xcZX33/9cnmbALSBZXfiWgYTHaiGkSvk3Z1NpWQ2uMkBDuuz9iWm1tlf5SzN16l0QpTRUh3Hd/UFCrVKm1NabVlhDAFO+b1Fp7m1U0BIPGElkjMQTUwAaCm7DJ5/vHnDM5e/bMzJmZMzNnZt7Px2Mfu/PrzDVn55zPua7rc12XuTsiIiIiIiIi3WhauwsgIiIiIiIi0iyq9IqIiIiIiEjXUqVXREREREREupYqvSIiIiIiItK1VOkVERERERGRrqVKr4iIiIiIiHQtVXpFRERERESka6nSKyIiIiIiIl1LlV4REel5ZvYdM3vCzE6P3PdWM/t88LcHjxcijxfM7Ekz89i2rjCzL5vZETP7oZl90swuSFGGWWa21cwOBO93YcJzft7MvhZs+zEz++XIY7eb2bfM7ISZ/UbsdauDx8aCMv+dmc1MeM43g21/28x+JvLYaWb2ITP7QbCNL1T7PCIiInmhSq+IiEhRAXhnhcefBpZHbr8BeCr6BDP7RWAT8FfAWcArgaPAF83shVXe/wTwr8B/T3rQzC4Otv0eYABYBHw18pRdwO8BX0t4+X8CP+3uA8BFFD/rDZFt/wJwM/AW4AzgZ4FHI6+/HRgEXhH8/h9VPouIiEhuqNIrIiJS9L+Bd5nZmWUe/3vgzZHbbwY+Ed4wMwP+HLjB3T/p7j9y94PAW4FnqVJRdPcn3P1DwFfKPOU64CPuvs3dJ9z9h+7+7cjrP+ju9wPjCdt+zN1/ELnrODA3cvtPgfe6+w53P+Hu33P37wWf68eAlcDV7v59dz/u7tHKtoiISK6p0isiIlL0APB54F1lHt8C/KyZnRlUjH8GuCvy+I8BLwE+HX2Ru58A/hn4hQbLtwTAzL5uZo+b2T+Y2WDaF5vZq81sDHiGYm/yLcH904HFwNlmttfM9pvZB8zs1OClw8A+4E+D9Oavm1lib7SIiEgeqdIrIiJy0vXAO8zs7ITHxoG7gV8BVgNbmdyrelbw+/GE1z4eebxeFwBvolhhfRlwKnBb2he7+xeD9OYLKPZqfyd46FxgBvCLFCvyi4BLKPYsh+87HxgDZgNvB/7OzF7R0KcRERFpEVV6RUREAu6+G7gHWFfmKZ+gmNY8KbU5EKYPz0p43azI4/X6EfC37v5f7v4scCPFccU1CdKW/xXYHNkuwG3u/niQBv0XkW3/CHieYtr2MXf/d+DfgNfV/1FERERaR5VeERGRyf4EeBtwfsJj/0GxAnsu8MXYY98C9gO/FL3TzKZR7J29v8FyPQR41WelUwBeCuDuT1Esd7ltP5TRe4qIiLSFKr0iIiIR7r4X+EfgmoTHHFgBrAz+jj/2LuA6M1tjZqea2XnAR4GZwF9We28z6wf6gpt9we3Q3wJvMbOLzOw04FqKvdLha08Jnm/ADDPrDyrcmNmvmtlLrGgO8H4mV8L/lmJa9znBLNNrI9v+AvBdYH2wTNNPA68B7qv2eURERPJAlV4REZGp3gucnvSAuz/s7g+XeewfKY67/R8U05m/QXHs7U+7+w9TvO+PKM70DLCHk6nHuPsdFFOqRyhOLHWUyRXzzwXPfxXFJYZ+RHHpIYCLgf8XbPs/KfZKvy3y2vdRnDX6v4BvAjspVoxx9+eBKymmO48BfwO82d33pPg8IiIibWexhmoRERERERGRrqGeXhEREREREelaqvSKiIi0iJl92MyeTfj5cLvLJiIi0q2U3iwiIiIiIiJdSz29IiIiIiIi0rVU6RUREREREZGupUqviIiIiIiIdC1VekVERERERKRrqdIrIiIiIiIiXUuVXhEREREREelaqvSKiIiIiIhI11KlV0RERERERLqWKr0iIiIiIiLStVTpFRERERERka6lSq+IiIiIiIh0LVV6RUREREREpGup0isiIiIiIiJdq9DuArTKWWed5RdeeGG7iyEiIl3iq1/96g/c/ex2l6OTKTaLiEiWysXmnqn0XnjhhTzwwAPtLoaIiHQJM9vX7jJ0OsVmERHJUrnYrPRmERERERER6Vqq9IqIiIiIiEjXUqVXREREREREupYqvSIiIiIiItK1VOkVERERERGRrtWRszeb2enAh4BjwOfd/ZNtLpKIiIiIiIjkUG56es3sDjN70sx2x+5/vZl9y8z2mtm64O7/BnzG3d8GrGx5YUVERERERKQj5KbSC3wceH30DjObDnwQWA5cDLzRzC4GLgAeC552PM3G9xx8hk0jyUsqbhrZx5Kb7i/7uIiIiGRPsVlERFohN5Ved/8CcCh296XAXnd/1N2PAZuBK4H9FCu+UOEzmNnVZvaAmT0w/sxT3LZ9b+mxaDC9dfteDo6NT3pcREREsqfYLCIirZabSm8Z53OyRxeKld3zgX8B/ruZ/TVwd7kXu/vt7r7Y3Rf3n/FC3rF0bumxaDBdMjTINIPhocHE7ai1WUREJBtpY/M1S+cya6B/0uNRis0iIpJW3iu9lnCfu/sRd3+Lu/9u2kms5p13BmuG55RuR4PpjtFDnHAYGZ3c0bx2806G1t3Lu+/czcGxcd57zzcUYEVERDJSKTYDeIXXqidYRETSyvvszfuBF0duXwAcqGUDZrYCWDF37uSW4jXDcyYF2tu2753Smrx114FJAXf8+ROTAuytQUt0dDsiIiJSWZrYvOSm+0sxNynOXrN0bmLsFhERict7T+9XgJeZ2ZCZnQKsBrbWsgF3v9vdrx4YGCj7nDXDc/jS+sumBNWVC2djQMFgoL/AqkWzGegvcOToBBu37eHg2DjvvnM38/54W2LvbzT1qtzfIiIivSZNbC6X3hzGUIAvrb8MQDFVREQqMvdKyUOtY2afAl4DnAU8AfyJu3/MzN4A3AJMB+5w9/fXuN2wNfltjzzySMPlDFueZ/YXODw+Ubp/1kB/KfhCMShft2U3J7z4mAMHx8aZNdDPkaMTHB6fYGZ/gYc2LJv0GvUei4h0BjP7qrsvbnc5OlEjsXnBhvtK8XdmfzFh7fD4xJQ4LCIivadcbM5NT6+7v9HdZ7n7DHe/wN0/Ftz/WXd/ubu/tNYKb/D6qq3J5ST1yIYtz+uWz2PVotml+6OTYG0a2cd77ixWeA14x9K5iS3W0QHLYSU5aXySeoZFRKSbNBKbo8LKbxhfFS9FRCRJ3sf0tlV0koyw5zU63mjN8Bx2jB7i4Nj4pEmwbt2+tzQWuK9gk54fio9DunX7Xk44TDOmpHJt3LaHw+PFlOroNtQzLCIivWbd8nls3LaHYxPH6StM59rl86aMA964bY/io4iIlOSmp7dZzGyFmd0+NjZW82urLZdQ7jnnndFX+ruvMB042Vv76o338+47d3PuGX1TZqyc2V/gBX2FSc+PtlbHp7JOmrky/jq1eouISN40EpvXDM/hoQ3L2HPDG9i1YRlrhueUYt2SoUH6Z0zj8PhE2ZmdFRdFRHpPbsb0NtvixYv9gQceaMl7XbT+Xk4Eu3Vmf4Gl887hrgcPTFl64TsbLy/9XWkM8DsiM1TGe3rj94et3OHYpujt4aFBtu46wMqFs7ll9SWlbag1XESkdhrT27h6YnNS3IrGuicOj5di8I1XzS89Z9PIvlLmFJxMiQ63BVqVQUSk0+V+TG83WblwNtMM+gvG4fGJKUsfASy6oDiOadPIPhZsuI93VxgDHM4uDZNnqEyadTre8xz2IB85OsGWBw9wwotLMYXKrXOolnAREcmjpLgVjX1hDF61aPak+Hjr9r2lCm84lCi6rXBVho3b9rT8M4mISHN1faW3kRSqet2y+hIevelyrl/xSmYN9LNy4exJj39n4+VsefurgclBGOCM/gL/9JXHSinQUKzort28s+xEV1HxivCa4Tmc1lecabp/xjSmGZPKUy6FO03qtIiISD2yHnoUjX1hDA4zmqKvG+gvMLO/wA2rij3AaSaZXLDhPhZsuE+xT0Skg3V9pTerGSLrEQbhW1ZfUurZPfPUAkPr7mXedZ8tViKHBotrAU8rpkKvWz6PB/cXLwIe3D9Wquhu3XVg0kRXtVRAw6B+/RUXT7oQCFO9jhydKPua+GRbqgiLiEijGonNSVlOcUmxac3wHF477xyePTrBlyOTT4aZWOuWz2PWQD/XLp9XeixsmD48PlGxwVlERPKt6yu9Wau3knfwmaMAPP2jCRwYn3Dec+dutgRjfc8+o5+Hggk5wgpywShVdGcP9AOw4PwB1gzPKaVhvfvO3VVboMtdIFQK5mlSp8NtKD1aRETypFxsChuQw2E+4fNuDmZ7js+dEe0djmdErd28k4vW38vazTub/4FERKQhqvTWqFwgrSasMC66YKCUOhUd53vk6ESpgrjl7a/mxqvmc1pfgYEgDWv/0+NAsfd308i+SSnR9bZAVwrmSZIqwkuGBplmk9cpBs0sLSIi7ZPUSLt2887SBFfhMJ/weUcnjpcqv3AyPgFcu3wep/VNXeExXoEOKbaJiORP11d6sx7Tm2YZoyRhhXHL21/N6MbLufGq+aUKZzjh1bvv3F0KkmEv7Gl9BdYMz6F/RvFf1T9jGrcmVHCHhwZLY39rCban9hVTquudqXLH6CFOOJPWKYZ0vcLR27pIEBHpHc2ebyOpkTasnE4zJo33dWDieLE2fHTiODA5PpWLVeGEWfF5O+ptHBcRkebp+kpv1mN604wlSrudXRuW8dCGZZwSrOULcN2WYsU3Xmm8/oqLmTXQz+tfeR7PHZ1goH9yq/OO0UMcHBtny4MHKgbbaMDOIjCXawRIO7N0dMkIrakoItIb2jHfRlhJXXD+wJQ4GNR5S/E4jE/DQ4OlmDs8NDhpQslKE2bF46LimIhIe2md3jYKA+6SoUHu2XWAieBf0V8wjh33SevphqJrET4+Nl66/8ar5vPuO3eXbofrA+8YPVR2LcNy6/+2Q9Kaw6H42sMiInmgdXob147YHMaUmf0Fjk0c5+iEs/CCAZ545uiUGBSNPw4cHBtnmlGa/RnSrXcfj2PR14DWBxYRyYrW6W2SRlpvwxbmkdFDnIjcf3TCOeGw5cEDUybJiLYgr1p0ci1CmLzMwuHxiVKv73vvfjjx9fHe2Ha2RFfqQa+UUq7WcxERqWdFAyhOKunAoz84MiUGbRrZx9NHjmIUhxCFr7th1XyAmrKm4nFM6wOLiLSWenob1EgvZLR388ujh9i660BpbFA4QUYoXM6oWityEgNGN15ees9yLcrhNuKt2HmmXmARaRf19DYuq9hcTyzYNLKP99y5GydYLQEmZVhF42p8u+Fj/QXj6ITTVzCuX/HK1HEzXDIQ4NjEccYnnIH+Ars2LCs9rt5fEZHa9WxPb7Mny6h3YqswoIU9rtGxQZcODfKCvgKFyH8nnKF508g+Fmy4b8oyRdfE3n/VotmlHmAHXnHdZ7lo/b289+6Hyy4x9NzRCYziMkn1jvNtdc9rtf2vnmARkfzJwySTa4bncGWQMTXhlDKsFmy4j7Wbd/Lc0Qn6Z0wrrXAQLlG06gNfLI3zPRr0FB877qwZnsPazTsZWncv8677bNWlBE/rK3B4fIK+wvTE9YG1HKCISHbU09smlVqlF2y4r7Qk0UB/oRhQJ46XJtgIH+ufMY0971s+5XXR1uIL19075b3DXmM4OY4oDLAz+wuc3leoe5xvHnpeoy3k4eeKlyc6nnr7nicBGprFWkR6j3p6G9fu2BztsR2fOHk9NM2KleBo7Lho/b2TMrDCia7CLK1bVl8y5Tk3XlU+a6rcXBbRXuB4XMpDjBURybOe7enNqzSt0kZxfcCHNizjzNP7Jq3NCzD+/IlJt9ctnzeltfiCM/unbPeZ8Qnee/fDvOfOk7NQhuVZt3xeQ7NT19vznaVoC3m58oRjqLY8eIDD4xMcHp/gui27Wbt5Z2JPuoiIdJ8wRly/4pXceNX84jKCM6ZxynSbsoZ9OPvzogsGSnHl0qFBzpnZz6XBWvXx5YsqZU2Vm8siXLLw2aMTU16jmaFFROqjnt6MNTIOJ9r7ODJ6qNT6G2/1jY7/vWX1JazdvHPS7bjoTJXPHp2Y1AoNxVTopNd1qkozQYfCXvHixc00nhmfwDnZug9M6jEXEYlTT2/j2t3TmySaNXXt8nkVY3r43DCD6tbte/nB4fHSagxhT28t1wabRvZx3ZbdU3qay0nq/dWYYBHpVeVisyq9GWsk9ajca6NB9aFIJSwMak+MjRP+F/tnTOP1rzyPHaOHWDI0WPodVqKB0sQdoTCF+rmjExx3uDJWCe7G4BmvGIe3h4cGuevBAzhM2d/x13fbPhGR2qjS27g8V3oLwVhfOBkP4uf+aAX51L4CB8fGMYpzaYQNp2s372TLgwcmbaeaNI23lZ6rNGgR6VVKb26RRtJ7q73WYrfDNN6+wslHxp8/wdZdxaWK7gqWLPq3PU+WUqjWDM/h/VfNL02SNWO64RTHCU94MVCHwTn+PvVMbpXXtKt4Wll4+5bVl/D+q+aXUr3LaWSfiIhIfoVDhY5HWoePTRwHpp77o8OKwhh+5aLZpfs2jeybFFOfGZ+oGg/jE11Wk5QmvWRokGlWXGopvu08xmQRkWZTpTdjldabrfe1SWN1YfJYpHCm5hnTrTQWaXrw3x0bn5i01u+a4TmcdUZxrO9ZL+hLrNyt3byTedd9lqF191IwEoNnGp1YOUzzP9Ss0SIi+dXIOTiMAVcuOjk+ty+YSDJ+7o/Gi/DvS4cGS9lUt8Zin1N9dYQs4uaO0UOccBgZPZT5tkVEOlHXV3qbvWRRGo0E30otvtFge8vqS/jOxst50Qv6GJ9wjh0/wURknqstDx5gyU33s+oDX+Si9fdy3hl9pcC9ZnjOlF7krbsOMB4sxbD/6fHE4JlGHia2gur/g1r/R9UqxrqwEBEpr9mxOYtz8C2rL+HGq+aXhgCt3byTW4NhMLcGSwiGojEknCjx5m17uGbpXGb2FxjoL7Aq6AE+94w+Llp/76TG6Ojr43GznmuIcrFXE2GJSK/SmN4WaMY433LCsT2Pj41Pur9/xrRJsz1PM3j0pstLt9du3sldDx5gusHpfQVeO+8c7tl1gAkvzgB93ElcVqFTxrVW249Zj3+qNB6r0nIUItI5NKa3cc2KzbWMia0muoxgVDReRGPIkaMTifNwhKJLCX5n4+VTXh+PQc0enxvffifFdhGROI3pbaNmjvONC3sgVy0qLq2watFsbrxqPqdMn/yvPmW6sXbzTpbcdD9rN+9kx+ghzugvMOFwWl+B7/zgSGkCjwNj4wwPDbJx2x4WbLiv9LqwNbvWlvR2tCpX249Z90hX6gkOl6M4PD6RuO/U6i4i0phGhhqlMc3g3DP6EntnwyFJ4dChVR/4Iheuu5dVH/giUGyEDn+v3byTi9ZXHkbU7Iyp+PbDXvKbt+1RLBKRrqGe3jbIYlmjWl4btuKWEy7TM9Bf4LS+4rqE775zd+JzoqLPD8uSpny9vrxC2NMbrsMc/7zRJaZO6yv0xD4R6UTq6W1cnmJzOWGv8Y+OTfD0jyZKszOHcbFSD+ymkX2T4umNV80H4OZte3AoLZcXim+rHbEx/Lxhj7VmgBaRTqKe3hwpN9YozbjT67bsrrl3NWzFDccTrVo0e9IY3mkUK7BDZ53OE4fH+fLoIS44s7/0uFHsGe4vWOl106xYYXvH0rmTxjalGUeV1GqdpzGwze5pXTM8h4c2LGPXhmWJFzHh/gHK7hP1BouIZC/p3Br2GkdTnAf6C0yjGB/D3tmk18Ynsrpt+17WDM/h1L4Ch8cnSqsvzAgmoIz35oaxcWMLe13Dzxv2WLd7Tg4RkSyo0tsG5VKVqlX8bt2+lxNerHDWEoSiy/GEv8/oL5QeD1OaH/reGCe8OIlVmNo8zeCM/gLjE84phen0BRXflQtns2Z4TqnM123ZnTgBRyh6MRBPO9s0so/njk4kBvxmKldxbHcFPM0FR3y/J1HFWESkNpXO/ysXFocNXbloNqf2FUrL/G3ddWDKBFaheKNzeD4/74w+AOadN5PzBvp5/rhzel+BL48emjTJVZpG0GaplCKu+CIinUaV3jYoF0iqzaoYPn7DqvkNpzmtWz6Pmf0F+mdMY2Z/geGhQaZPK7Y4zx7o57mjE/TPmMYL+gosnXcOBSuu5RvO6HzXgwdKZTKKKV43b9vDmuE5U3p/ofKFRDjG9fS+QkvTeMuVKS8zTle64Lhm6dxSal2lRpK89J6LiLRT2kpapfP/Lasv4dGbLueW1ZeUZmUO49/GbXt4JugJjqYrxxudw/P5Q98bK/2OvufWXQc44cUVF6KNxI32umZdSVV8EZFOo0pvDoTBCJhSyYkGlqRKUL2BLEyx3fO+5Ty0YRk7Rg/x/PFiqN7/9DiHxyc4NnGCw+MTjIweKvX8hsK1BtcMzymlZx2dOD6pzNFeyEoXEu2qZJZ732ZPgJKFNcNzuGHV/Ir7bcnQYKr1ldViLyLdLm0lLe35f83wHNYtn8cZ/QVmBplTYZhcOu+c0vPi59fw9uyg93bB+QOlxuKN2/ZMuigLe4wrLV2YVtaV1Lw0DouIpKVKbw5UCkbl1utrdAblpPcJx+uGvb+nFKaVUo7PPLUY1E8tGANBkB8eGixV1gFOKUwvbSupF7LclGlJ6c6tqISlvbhpdnkaabioVP4do4dSra+c5mJIFWMR6WTNqKRFs5TWLZ9HkCw16ZwbP7+Gtw8Ek0s+9L2xUqX28PjEpAbm8M+ktOlatXKFAsULEcmjrq/0mtkKM7t9bGys3UUpq1IwigeWMGBu3XWgNCNzFoFszfAc3n9Vsefw+isu5rS+AuPPn+DZo8V0rXACj/EJ59S+AheddTpbHiyW4eiEM7O/UFqeYc3wHBacPwAUl3SIlrvSGNRQu9Km2jXGt1nbT3ORk3Y8tVLZRCRLrY7NzcjgiZ5jo9k3YYPw2s07efq5Y8DJjJvwNeH44LBx+JrI+be/YMwa6GfpvHMmrRHcyFobrcxgUrwQkTzq+kqvu9/t7lcPDAy0uyhlxYNRuG5fOJFFVDRghusAJo2hbbQc8d7asBI73YqTaTy4/+SFikNpPG5YcQwfD8ctpRmDGv+MrU6bSgrUYaVwoImTbDXr86a5yEk7njptGdXCLyJpdEJsriQp5Tg85+4YPcTBsXHuevAA48+fAE72/kbH+EaHqKwZnsONV81noL/AKYXpvGPpXHaMHpo0Y3SYNp3386xSn0Ukj7RObw5dtP5eTnhxKYRzB/qrrs8XtgSHva1ZrekXrtUXVqoPjo2X1ub9/uHxSTM8r1w4mx2jh3guWNevf8Y0jk2cYOXC2dyy+pLS9sK1CdclrE9bqRytWKcw+nnD9+n2NYWTPnP88Vo+a9L+anSbInmldXob10mxOarSuS48rz515CjjQaBctagYC6ud/6LbfcfSuVy/ZXcp1oa9yFsePAAUl006tUPXclccEJFm0Tq9HSRMe+orWE0pQka2aUXxnt9ZA/1cu3weX1p/Ge9dNZ+Z/QUG+gvcsGo+2/c8ycGxcY5NHC+lSD960+VcGknzunX7XpxiqnQt5WtVqlRSz2jSZFDdlLpVrTe42meN9zikaeHXGGIR6XRphiVdv+KVpXG+2/c8yYIN9/HuO3dzcGyc99y5mwUb7ivNzxGe66IxZ83wHE5EthvO7hxy6l/GqB3n2Oh7dlMcFZHOoEpvDoXLIly/4pWTxgeVC07hUgbXLp9XCsTnntFXNkW6HkmVo9P6ClwbjOMNl2roK0yvOAbZmDoGuVrwbWeqVNJkUL2UulXts8YvXNKkVGdVMRYRaUQjFb8057roOF9gUqpy2AB8VzA3Rniui8eclQtnl17z5dFDpUbxVYtmN7SMUTvOsdH37KU4KiL5oEpvjoXLGIQVxnLBKR58HXhw/xgnnEmtwlkKZ5PcuG1PqQfXgvfeNLKvdDGxZGhw0hjkaxPGIFcLvu2a3bnceN5oebq9R7LahV3aybKi+yirinHStkVE0mpFxS+6zu5Af4H+gk16vC+YtCo818XPfZdGsoy27jrALasvYeXC2Wx58ADvvecbdS9j1I5KZ3zir6Q4oHO6iDSLxvTmXDi+Z5rBDavmA6QaD9Q/Yxrjz5+gf8Y0rr/i4szHzITjiAf6i729t23fy5FgPO+sgf5S2lVY7kkpw0EZw/FIS4YGGRk9NCl4Vxrvk2bcaBbSvE+rytLJmrmPNIZY2kljehvXzthcbU6DWrazMVhOqNp8FeE5yyhWeE8pTGfpvHPYMXqIJUOD7Bg9NOlcFT4fTo4LDuf9gM4e15skvLboD/YN1DYHiIiIxvR2qLBlNKw4Vlv6J3z+9VdczHkD/Yw/f6IprdjRlOqwRxooLX0Tn6052nobljGsGIcV3rS9v0kt1M1oHU7TEl5uHWW1Up9UbT82ss+UKi0i9cpqGZ9wFvw081WE56z3XzWfM0/v4/D4RCmbK1wGMBrfr1k6l5n9BWb2F0q9vtHVFMbGJ0pZV90Ue45OeGmfXrdlN2s372TBhvtYsOG+rvmMItJa6untMJtG9nHdlt2ccKr2nEVnS1467xy273kSaE6raXQG6Yc2LAOKSy9t3XWgNLNzvOe33OzQee79rSYv5egk1fZZoz21aXtz1CMstVJPb+O6ITaHPb0GpYbgtK+7bftezj2jb9IygFDswd0VxNL4OTKaAXYidgkXfV0WWn1eDPfJ8NAg2/c8yTPjEzhM+qxZf0YR6S7q6e0CYfAJx8dWm+BqzfAcTu0rlFqSw1bTMA0ri/LE3z86Wik6IUdSz290jcNrls5loL9QarXevudJ4s0xtfb+tkO0HOr1TafWybJqlbY3Rz3CIlKPNcNzeGjDMnZtWFZTxTA8Nx185mjpvnDI79j4RGkiyvgqAuE5M7wWiMbd8HX1xJ6kmNXq82J0HeOHNizj/VfNL33W8HP2RleNiGRNld4OEgafkdFDfGn9ZaXe05srpDVFg2Po2PHiIgiNVsqiwTCa7hxu+7mjE6V05+gsltGe3TCNK6ygA6UlHuKBNmn5oFC7JruqVA5VotJpxmRZ9WjV+4iIRIXnnhuvms97g7k7ALY8eIAFG+7jrgcPlBqQ127eyXVbdjM8NMilQ4M4sPCCgUnbC9Olyw2DKicpZjW7QbnaOTVaCQ4rwOuC6wwRkVp0ZKXXzC4ys4+Z2WfaXZZWigef+NjYpMpVNGDM7C9WKk+ZXvy3R2dgrkd8PcH4UkWHxyc4va8wqTITttAm9fyGMyXfsGp+qRId7c1OWj4I8tE6nUTjfbORpqc2i/93Vu+j/7NIZ8ri2K1nG9Fzz5rhOaxadLJX83CQ3gtw7hl9bN11oLQyQ3g+euKZo3xn4+XcGOkVDeNrWPFNU66kCm5W457LqeXc3eyyiEh3a3ml18zuMLMnzWx37P7Xm9m3zGyvma2rtA13f9Tdf6u5Jc2f+Ak/vL103jlTekCTAlxYkVy3fB5rN+8srRkYTY2qJWCXq4RCcq9sNLgl9fweHp/gmaMTkz5b2Jtdbl2/cIxzq1un00hqCGh3RbxbNXOyrFreB/R/FulUWRy74TYamVjqltWXcG6wtm/Ug/vHShXalQtnJ56PnOIyR2F2V9iwnOaztaNSmfW5W42OIlJOO3p6Pw68PnqHmU0HPggsBy4G3mhmF5vZj5vZPbGfc1pf5HzZNLJv0iyGSZXPpAAXDWjR9XuvjaQK1RL0KwWraJnia/aGz4+WJ9rzG03Xjq/rlzTL8wkvpkS3snW6HnmoiHerav/vrCqiWmNYpHuVO3ZrOVbDbUD5DKxaynLjVfPpn1G8VJsx3fjX3Y9zwuFfdz8OMCkmRs9zO4LrgTA2hrNAP3Xk6KQZkCt9tlaco7I+d6vRUUTKaXml192/AMS7Bi8F9gY9uMeAzcCV7v51d78i9vNkq8ucN/HlEZICdbUL77C1eNWi2ZPGwUbH4YbKBb5KwSr6/vGxyEnjbqM9v9F07Wq9pfElncqVNQ8VjDxWxHtFK8frZjlxVh6+tyK9otyxW08K7rrl8xjoL3Dk6ERdx2+0LNdfcTGzBvo5dcZ0xieKyc7jE87GbXsmlS08zw0PDfLc0Qn6Z0zjBcFcGWuG53BaX4HxYCmg8LNU+mx5qEDW2lisxmURKScvY3rPBx6L3N4f3JfIzF5kZh8GLjGz9RWed7WZPWBmD3z/+9/PrrRtFs50HJ0kKh6oq03sdMvqS3j0psu5ZfUlpcev27I7cRxuPYEv3os7MyH4x7cbvVgoNxY2HtDSphDnIXhHqTLTWlmM1836f6ZUaelVnRab66lIRVdPyCrDZN3yefQXTg5IMk4OJTr3jL7Sigg7Rg9xeHyCYxMnJr1/GIsHIg3bS4IhSE89d2zKuS0PFchaG4vVuCwi5eSl0msJ95Wdld7df+juv+PuL3X3myo873Z3X+zui88+++xMCpoHa4bnsGvDstJ6uNEL8XIX5tUunsulCdfa+1uuvKdFgn+5dOfo88tVZJMCWqVKcSiL4J1lpUeVmfxp9tJJcUqVll7VabG53opU1pXGNcNz2HPDG0oTVl27fF5pZucH94+VVnOIL2kUbSSOL60UpkGPP39iyrmtXRXIZp7PdK4U6V15qfTuB14cuX0BcKDMc2tiZivM7PaxsbHqT+5A8aV/yl2Yl5sEKl5ZjKcJb9y2J1Xv79rNO7lo/b2ldQXjktKdt+46UOqprqSWyki5SnF0TeD4Z08ry0pPHlrQZbJGl05qxsWUUqWlW3VCbG70mGlWpTG63XjvgEcev2X1Jaka1uK9v+W06hxS7nyWxfurwVmkd+Wl0vsV4GVmNmRmpwCrga1ZbNjd73b3qwcGBqo/uQNFJ4DauG1PYq8sJAffSpXF8DGDxFTjaC/tppF9bHnw5DIKSYGp3KRV0Z7fauv0pa2MxLeXFOTqCXxZVlTLfSZVTPKrVZNl1UOp0tJpOiE2Vzpmaj1XV2sYrtcFZxYnzTrz1MKUNWzTlDGp97ecVp1Dyp3Psnh/NTiL9C5zL5tF3Jw3NPsU8BrgLOAJ4E/c/WNm9gbgFmA6cIe7vz/L9128eLE/8MADWW4yNzaN7OO27Xs5crQ4udWsgX6+tP6y1K8LT/63BhNhhD29t8V6R6GYSn1wbHzSe4T3QXFirHCZoZn9BU7rK5S2We69w0CWttzVxMuY9FnKfb5q+yu6j5ohaf9KZ6j0ndo0sq+0Hva65fPaMt6sXPnyULZOZWZfdffF7S5HJ8tzbK50TNd6rr5o/b2lIUSP3nR5s4rMqg98kQf3j7HoggEOPnO0aiyuRT1xM0vNfn+dC0W6Q7nYnLrSa2bnAq8DFgJnAk8Du4D/4+4HMytpxsxsBbBi7ty5b3vkkUfaXZymCgPC8NAgO0YP1RTgogE8rIhGK8DhbaBqBbLWSnjWgaxcJbdahbXac1pRIW33RYU0R7RhKG8NGtGyhRfHS+o4h/QiVXrr1+mxudZz9drNO9m66wArF84uTSDZDBeuu7f096pFsxkZPTQpFg8PDbakHFlqRYMz5Ps8LSLplYvNVdObzewVZvYZ4BvAm4AZwMHg95uAh83sM2Z2ccZlzkQnpFBlJUy/3L7nSQ6OjZdaLNOIzrC8cdueSSlE1cbMwuRZx8rNwgzJ6VbVZpquVbVU7nKqPacVaVGaebI7xWdcT9Ku1PZo2YDSeHulQkszdXpsrjVuJa2Y0OzjPVwmMBqLt+46ORSpUd0wxjeq2vhmDT8S6WxpxvR+nOK6ubPd/XXu/k53vy74/TpgFvCPwMeaWE6pw7HjJ+qaYTk6jjecvTkaBNKMmQ23mabymXZ7jag27jfpOXGqkEq9ojOu521McLRs4cXxgvMHmGYwHCxnoos96XXVjoFaj99mHe+rFs3GgP6CMTw0yIIN97Fx255Sj/TKhbOZZrBy4eyatpv0+bthjG9UtfHNmhdBpLNVrfS6+7C7f8bdj5rZejM7Pfb4MXf/tLv/VPOKWb9OmCEya+GF6ynTp3FwbJx337k71eQZYWC5dvm8UuXu1u17OTw+wWmR2ZvjJ/7wdcNDgyy56X7Wbt5Z9uIgKXiV216WPapp1vNVpVZaLc1yW0nPbZbwGDj4zFFOeLGnCHSxJ9nrtNgcXykhrpa4VW4pwCzcsvoSRjdezp4b3lBar/fw+ATvvfthltx0P5cODXLDqvnsGD006XPUU6lv1aRQSbF57eadPDE2Tn/BWjYpVfxaR42AIp2l1tmbJ4ARM7s2rPya2Qwz+0L2RctGp6dQ1SOaXhxKk8qUFFiiQa3S+roOpbTqux48UDa9utp7hM8JxxU3K6hUW8JJpBXSDB1Iem6SLL+/8eMjTaaESC06LTbHVx2Iq6XRNGxMDpcCbNbxdE0kvo1PeOn8Uc+KBkkxM+uhSbXYuusADhw77i1rqA4/bzhZZxazeotI69Ra6d0fvOYm4DEz+3pw34+yLpg0bs3wHFYtqi+VKbqNaK/vwbHx0jihcssb9RUMAAu2Uc+SRM3uWWrkPRXUJCu1NL7Usl51o+LHR7VMCR0T0u3WDM/hhlXza+rpS3ssZx3vwveFYpozQMEolT06ZKlSY3ZUmkp9KzNCktK0W3UeqnQuVlaMSH7VWum9FfifwLnAjwN3AkeAP8q4XJnptBSqrGU5eUa51J5oWvQ7ls7llMJ0ZvYXuDboac5iTdxWBLO076mgJlmppfGl1vWqo5ox4UszL9ql+3VibE7T0xeV9lhOyqhq5FiNvu8phekAnN5XKJX98PgEHjwvnLgy3pgNtZ83WrkGbvzaBlp3Hqp0Lo5OCqpGQJF8qbXS+yzw7+7+fXf/nrtfD7wR+Ofsi5aNTkuharZGgkK5gB/vDT48PsGzRydKr1syNFjzpDhpxuBmLe17anF7aaZ6v1+VLsSaMeFLMy/apft1cmxOe4ymfV5SRlUjx2o05oZzfISN0EuCOHx4fKKUoVWuklZrWZLOQc04H9SbDZPV+1QSnRS0FUNRRCS9Wiu9dwD/YmbRq6p9wIuyK5I0U9qL00qPVQos8XFPm0b2lZZIaGRSnHaMwS33nrcGayFnOeY4/lkUFHtXucprFlkatU6SlfY9s75ol/wzs3PN7E1m9mdm9tHg95vM7Lx2l63Zwu87UPPQnSS1TGiXxo7RQ6WYGy/DjiAOA0wzeO28c3g2WMc3TQNvreehZpwP0vSgZ91jXotWDkURkfRqrfTeCHwN2GNmD5rZZ4Pbd2VeMmmKtBenlR6rFMjj456u27KbE14MrmEAiPf81lruNGXMQqX3zHod0/hnUVCUuCyyNGqdJCurBirpHmb2CjP7DPAN4E3ADOBg8PtNwMNm9hkzu7iNxWy6TSP7uG7L7kzO07VMaJdGpWNwydBgaUmjcBbnMEbHV2AAGo67zTgfpNlmFjG0GVk3jWxXRBpTU6XX3Y+7+zrgHOBPgLuB3wbe1oSyZaITxw1lrZ5UoEZOytE06DCY3rBqPlBsFd++58lJPb/1akfgCN9z5cLZmS5dUG22XJFmZjuU+75ds3QuA/0FnjpylAUb7ku95reW/upqHwc2A7Pd/XXu/k53vy74/TpgFvCPwMcqbaTTY/Ot2/dOadBNY+3mnVy0/t5Jywhmfb6vdAzuGD2EAy88vY81w3NKY1Bf0FdItQJDrWVtxgzPac4xWezTZp7LPMVzlPElki1zT3Podb7Fixf7Aw880O5itMWSm+7n4Ng4swb6SylZScLU3WuCReyr3V/N2s072brrACsXzuaW1ZeUyjHQX+C0vuLMkZ1+Ybxgw30cHi/OhLlrw7J2F0d6TNpjO4v3AKa8T/wYh/rPF53IzL7q7ovbXY5WM7P3uPv7s9hWp8bmTSP7uG373lIcS/u9v2j9vaXK8qM3Xd7CEhdFyw3FyvtzQXrzzP4CBoyNF+fkGAgmpKz2uWo55ltxzkqSp/NS2n3Qrn0l0unKxeaqPb1mdo2Z9VV5Tp+ZXdNIAaV50rZ4ZpniCJPHFUXLce3yeanGQ3WSaNORWmelVVqRERD2BIVLnESF4/Wj64CHs8Em9RJJ17i23QVot3ITH163ZXfFc3/SUjtRzY4fSUOcngkqvOuWz+PUvkLpua+dd06q+F/LNUK7spjyNGQo64nQRCSdNOnN5wF7zewjZrbGzH7CzF4e/H6jmX0EeIRiyrPkUNoUnXKTXMXvX7DhvlSpjvHXxSeAylMQqlc4M+a6YGZMyFdwle7WjEmvkt7joQ3L2LVh2ZT3qXQBb5G/1RDUdaz6U3pLfBLHcpKW2olqZfwIy+wUlzQK052nBf/dkdFDqebgqKVy1q7hD3mqQKbZB+E1UrWMOJ1bRdJLld5sZmcBvwEsp7g+75nAU8BDwGeBT7j7D5tWygx0agpVu5RLq4mmOoYXuztGD1VNGQrTgENh8KknPSxv4uWuN+1NJCvtSIsLv+dLhgYZGT006WKtW9P0eji9+bC7z8xiW90Qm+Pf/eGhwVRxsdy2bktR2clKUrpz9BgOK+HRuTnSxrO8xr68litKKdAi9as7vRnA3X/g7n/m7pe5+znufoq7n+vuv+Duf5nnCm+nT5bRLvFW0bA1ccnQYGnczwmHLRUmvEjSP2NaabvtWIu3GeLlrva5am2ZVUuu1KodS3yF3/OR0UNTejHKnU/0ne5YBTN7rZktLfdTbQPdFJvj3/34WvaVxI+FLHtC02RmJaU7b911oBSjrwmO2RMO77374ZpmrE4T09txLuiEaw2lQItkr9YlizqOu9/t7lcPDAy0uyi5UM+6mzA5qD+0YRnvv2p+KQUKque6hWnA119xcdmZHNOuIZw3YbnLzeQcD0q1BtxOCNCSL+1Y4iv+PY/OUpumgauTjnnhSeAOijM0J/18tNoGuik2NzL7ftqxwPW4dfteDo+fXIO32jGWlKK9ZngOM/uL43yPTnjijNX1rBARLWOr41snNMKlbfyoZy1okV6l2Zt7TL2pMEkpV5tG9nHztj04sHTeOanSuaJpRUBpHd9yadSdlLJTyyzZSfuyXLpVq9PdpDvV+r1rVKVZapPK0onHfK+mN2epm2NzmuMrXO83KQ6m3Ual9w9j9LpgFuZqx1ildOfh2NCFsGzh7M/1HLt5iG+deO6p9XvRiZ9RpF7lYrMqvT0miwCTdLKtZ/yJw6SxQvFlktodCGvVSJkVkKQdkr53WVWEo8sZXRpMYFft4r+VlfIsqNLbuG6OzeFcFjP7CzxUYUm7SrEjy9hQa4yq9N7RyvrM/gKnx5YhzPuxG9WJ1xu1fi/SfMZO+p+JVKJKbxcH1lYrd6GcJmjEW5E7LdBUU2/Q0ORX0g6t6nGtd5t5bwzqxUpvsDzhR9z9aIXn9AG/7e63VtteN8Xm+Hm73nXc4xlR7YqTleJSfIKrrCrrtca+Xo2Vzaio5/18K5JWQxNZiUQljYdJM7U+TB1/0m1NLvWOT0oz5jE6RlKkUeWO26RxeI1+9+odQ6dJWnJJyxiWER+fG85lcW1kSbs0wrWub962J5OJreodz1kpLoXHZlKFF+o/dts950WnjH2t9r2o53Ok+Z91yv4RSVJTT6+ZnUJx6aJFwAuij7n7m7MsWFbMbAWwYu7cuW975JFH2l2cXMmqhTRsHRzoL3BqX6Gnx5jEx1A1kkIeb8WtNEZSpFa1HH/x716j5474e9eyvTz17PRiTy9ks4xhN8bmauNzw+dU+/5WSouu5/ufVaxtpHcxbblrfY+sezy75bqkWZ+jW/aPdLesenr/DlgLPAN8O/aTS900Q2TWsmohDVsHwzG6aWeh7NTZmitZMzyHU/sKpdkyG9lOvBV35cLZpbWRRRpVS09M/LsX9kRVWqqs0jHdyIzmms28/bJYxrAbY/Oa4TncsGp+4mz+4fEQHjuVvr9hD/G6hB7ier7/WWVMVOpdrBbD05Y76T0qbXvN8JzSesJZXD90S3ZJsz5Ht+wf6U219vQ+BQy5+9NNK1GTdNO4oaxk3UKappW7nHnXfZbxCae/YOy54Q0Nl6Wdwv06PDSYakZrkU4Q7anZuG1PaazitcGssPHveS09ArWci/I06Uyv9vRmqVtjc/z7H82IOi026VM1m0b2lRqYls47Z9IMyrVoZpZEteO9kYnqqm1bvY/1aeb3IU8ZOdJ7surp/S7Ql02RpN2yGCsU317Yyl3r2L2jEz7pdycL9+uO0UNNW39RpNWiPTXRsYrlenCqrV8dVcu5KOvzlkgzlFu799rl82r+/kbX2x0ZPVR6fa0ZUs3MkqjWA9jI2uHVtq3ex/o08/ugjBzJo1orvZ8A7gomq1ga/WlG4aR1skovTjMhU5IrFxXTJ69c1D2pu9csncs0gxNO5if+bkkHl84RvbCMHuflLjjjjT9pjwF9t6UbRI+RSpM9pvm+X7N0LgP9BWb2FyYdZ7VWLJYMDTLNYHhosL4PVUE9jVFpK6vVtq2GsPo0s7FADRGSR7WmN4+Wecjd/aJsitQc3ZpClZWs1wJMWm4BSLVWZzelxDQrFVPpXJIHScdr/L4s1wbNG6U3N66bY3N4LDx3tNhLm/SdbuT7Hl0L+5bVl1R9fiMTyNWj2+J5r9H3QzpVJunN7j5U5ifXFV6pLstWuWjrc7QFNt4qndTCHV/yodOFnx/ItPdKraiSB0k9TfH7au2F0Xe7c5jZX5jZonaXI6/CYwEo+52u9/u+aWQfW3cd4ITDyOihVK9pZAK5enRTimsvZqDo+yHdpuZ1es3sZWZ2fbBG3/Vm9rJmFExaK8v0oHJBPE3AbWZKcDtlXZlvVmVapJroxV/Ssd5opVWpih1lBnCfme02s2vN7IJ2FyhPwmNhXTCOF6aer2udrTh06/a9pWXE0s6hEX+vZjcwdVMDVi9W0PT9kG5Ta3rzCuCTwD3APuAlwBXAm9x9a1NKmJFuTqHKUivTTaJpj8CUlOhumv24kZmtK+mkVFDpDuW+c+VSncNZZxtZtzqvlN4MZjad4nq9v0rxemCE4vwf/+Luz1Z7fS/F5rTn6zTPSxo2EL5umsENq+aXPd5qTYuWfM0anxetul5UGrTUKqvZm28ErnT3Ne6+3t1/FbgyuD+XzGyFmd0+NjbW7qJ0hKxbM6utr5eU+lzvBDh5ljSzdRbUUiqtVu47Vy7VOZx1tt4sh15MK+wk7n7c3e9x9zcCS4CzgY8DB83so2Z2ftLrejE2R2c0X7DhPhZsuC/VOtZJknqI02ZKhWnRW3cdaOwD1ajZx3Izt68MlKla1fvdi73s0hy1VnovAP4jdt8Xg/tzyd3vdverBwYG2l2UjpB1JareJQk2jezj6SNHMZoz02Q7NCMlWYFYWm3N8BzesXQut27fO+k7XC7VeaC/ANQ/ZEEXPPlmZjPN7LfM7N+AL1Ds6f0Z4BXAs8C2pNf1amx2YPueJ0uNQUnf6/jMz2niRdgbtnLh7NJxWO61KxcWV0tYubC1qyVojGh3aVWje5r3UeOopFFrevO/Af/q7jdH7vsj4A3u/prsi5edXkqhykoWKSX1pgSFaVqQbTpwHjQrJVkpQNIq1b7Djc7gHN9WXtMKez292cw+AyyjWNn9BLDF3Y9GHp8GjLn7GeW20UuxOTxuBvoLOGDAtVXS/qPHWtjYlHSOTzom06Y7t0qzj+VuH04h5Wmol0Rlld78e8BbzeyAmY2Y2ePA24DfzaKQki9ZtJrW2xMZXZdweGiwq1rwmtU6qlZuaZXw+DxydCLxuEyawTmpdzgNZTPk2g7gZe5+ubv/Y7TCC+DuJ4Bz21O0/AnP/dcun8dDG5axa8Oyqt/raLyodI4vl2lRLd057CFbu3ln0+Nss4/lNcNzOK2vULYHvR7qQWxMq/afhnpJGrUuWfRNiilLvwz8OfBLwMXB/dJl2nkSWTM8h10blvHQhmVdNbYX6ktdS0MnfalXrd/DNcNzOLXCxWXSd3Hjtj0cHBvnPXc2NoO5LkLzxd0Pxu8zsz+IPP5ca0uUX/VU+qINRkuGBsumLse3HU93jjYeR18bVqS37jrQljib9fHcriFakqxV+y/NsaXYIVUrvWb2s5G/lwI/C5wC/CD4/TPB/dJlmlU5q1W8V6mbTlxZBgQtYyT1qud7WOnistJ30WlsOTJdhObK9WXuv66lpehy4Xd+ZPRQ2XXvq70m2ngcfe2SoUGmGSw4f6AtjaZZH89Z9yarMbkxedp/ih2Spqf3Q5G/P1bm56PZF03ypJ0ni3ivUrQsnV4BjgaErD6LTuxSq3ouTNJcXEa/i+uWz2Nmf4GBBocs5OkiqleZ2dKgsXu6mb02vB38vBV4pt1lzLtazvf1rIcdfzzaeBztMd4xeogTDqM/OEL6GV6yk/fjWcMrGlNt/7XyGi7td63TryulvJomsupkvTRZRjO0ezKZ+Jq+4d/hRXU3TF6Q1UQM7f5fiYTKfRe7ZdKRXp3IysxGgz9fAnw38pADTwA3ufvWNNvq1dic5TGQdhLDpEmxlgwNMjJ6iCNHi7NJd/oxKZ0lj7Egj2WS2mQykZWZvavM/X+QdL90j3anzsZbC8OmmqSljjq1hS6rFu+8pKVLb1i7eScXrb+XtZt3TnmsXCt/3nt3pDJ3H3L3IeCT4d/Bz0Xu/lNpK7y9LMtjoJ6lAePpz+uWz8vFMZlVzFLs6wx5jAV5LJNko9bZmzV+p8flIXU2Wob4RXUeylevZqRRdfL+kM6wddcBTnjxdy2iOUa6QO0c0Xk+gI/HUpuXRlKfpYJGZjSPCy/Sqw0biMaYpBnY43l/7Tgus4pZ4Xau29LYxHlxOldlK48TUOWxTJKNVJVejd+RUB5awCqVIQ/ly5NmjBkWiVq5cDbTrPg7rfiFrRpnOorm+chIVvNThBfptax0EJ0rY+O2PVy3ZfeU1yYdl82OI1nF8DTLNdVD56rWy+M+z2OZpLq0Pb1hIOsH7mByYPst4B1NKZ3kTh5SZyulWmvSicmi+0MnaWmGW1ZfwqM3Xc4tqy+p+Lzo+SJpkh01VnUGd58f+XuozM9F7Sxjp0i7Bm8926vl+VCsHE4zqk6W1ew4khTD67nWWDM8hxtWzc/8vKJzVevlcZ9rUqzOVNNEVmb2CXd/cxPLk5qZrQIuB84BPujun6v0/F6dLKOZ2j3Yv93v32k0wZVUk3ZCnCRrN+9k664DrFw4O7ECnPZ4baQMrdaLE1mlTV129+1pnqfYXNTO83P43sNDg+wYPcSS4HfSMdiOcirWS6fSd7c9MpnIyt3fbGbnmtkKM3uLmf1m+FNjYe4wsyfNbHfs/teb2bfMbK+ZratSli3u/jbgN4BfqeX9JRvtTp3NY+tfnrV7MjLJv0Z6caqN7U17vCojIffKpTQrvbkBabOUmhFr4+nRW3cdKHsMxsvZitivWC+V5Lk3Vd/dfKl19uZVwLeB9wIfoZjW/BHgTTW+78eB18e2PR34ILAcuBh4o5ldbGY/bmb3xH7Oibz0uuB10mLtTp3NQ6p1J1KlQsppJEBXG9ub9qJeFwn5ViGlWenNGakUz5p5/l4yNMg0g9kD/UwzGB4abGt5Qo1WtHV90N3yfE2jSbHypdbZm28A3uLulwBHgt9XA1+tZSPu/gXgUOzuS4G97v6oux8DNgNXuvvX3f2K2M+TVnQzsM3dv5b0PmZ2tZk9YGYPfP/736/xo0ot2n2hmueTXt60u4de8quRMfHh2N5Lq8wgm2UZ9P3tTIrNlVWKZ82MtTtGD3HC4cDYOCccRkYP1VWeZh+Xtcb7Zs3kLPnQ7uvPRpX7Piu+Za/WSu9L3P3Tsfv+DshinO/5wGOR2/uD+8p5B/DzwC+a2e8kPcHdb3f3xe6++Oyzz86giFJOu3tdO/2k10rt7qGX7tbK75S+v60XXbKo3HJF1cb9KjZXVimeNTPWhu+7cuHsshPNxd8zqZGq2cdlPRN2NWMmZ8mHag2lea88lvs+q7Eme7VOZLUX+Gl3f8LMdgK/B/wA2OHuL6rpjc0uBO4JZ4I0s18Clrn7W4PbbwIudfdMZobWZBmto4H7nUOTW0nWNo3sY+O2PQCsWz6vqd+rdn9/e3Qiq92RuD1a5mmeNsVZsbl+rYi18Ynl0rxnu4/LJHksk7RGp16TbhrZx3VbdnPCYWZ/gdP6Ch0xwWMeZDKRFfA3wKuDv/8S+DdgF5PX7avXfuDFkdsXAMkzotQgmHTr9rGxsUY3JSmFrVbDDaY5SvNpcispp97W8TXDczgtWP8zqVcly1Z3LVHWelktWaTYPFml46LcY5V6PLM6zuK9tml6WdtxXFb7vM0oU957EKWoUzMBo8tuARwcG+fmbXv0nWtATT29U15s9hLgdHf/Zh2vvZDJPb0F4L+Ay4DvAV8B1rj7w3UXMEKtya3Xqa1rvUj/K4lr5DtRqVellu3mffmiXuzpjTKzUyhOJvlGYDbFhurNwPvdfTzNNhSbiyodF2mOmfiSYVmd0zulhzRtD3SW5xPFTWmV8Dh86shRxiecmf0FHtqwrN3Fyq2senoncffv1lnh/RTwJeDHzGy/mf2Wu08AbwfuA74J/FMWFV61JrdPp7au9SJNbiVxjRy/lTIIatmuxuzm3oeBpcA1wE8Gv3+OFNlfis2TVTou0hwz4ZJhdz14gCU33c+SocFM4m+nZAOl2Ufh+WRjRr1lusbpfJtG9rFgw30s2HBfbr/bcPI4PKUwHQAL7tf1Wm1qHdP73nKPufv1mZSoSdSa3F5577GRk9R6LVlp9LuU914m9fTaD4GXuvvTkfsGKa7EMJhmG4rN2Qh7ek+ZboxPOAac0V/IbFx9p8WFpGuO8Hxy5OgEh8cnOuazSPOE32ugI74P8ZjYacdlq2TV0/vi2M9PAu8CXtpwCaWrqcemc6j1WrLS6HdJY3Zz7yBwWuy+U4HH21CWnlCuZydcMuz6Fa9kmoFD2XH1ad4j3vvVaXEh6ZojPJ+sWz6voz6LNM81S+cy0F9gZn+hI74P8ZgYPy7V81tZQ2N6Aczs9cAb3f3XsylStsxsBbBi7ty5b3vkkUfaXZyeFW2dAtTrK9IhOilLo9Vl7cWe3thyRJcCa4DbODkZ5e8Dm9z95irbUWyuQ9qxqzdv24MDS+edw47RQzUdE9V6vzrhnBDvEeuEMks+ddJ3J35+6KSyZ6lcbM6i0jsNeMrdBxraUJMphSo/lI4h0jmaebx2+sQyPVrpLbdMUZSWLGpApeOi1pT/eo6JcNkxA65NSI/uxBherjKwZGiw5kYB6S2d9H1X+nNRJunNZnZR7Gc+cAPwWFYFle6nJY1EOkcz0xqzHvbQaSmYnajCMkU1LVkk5VU6LmpN+Q/TN48cnUgda9cMz+GhDcvYtWFZ4vt04nEWL3O4j7fuOtDwOUgppd0n+j/N4/e93HeuUvqzvqe1T2R1InbXc8CDwDvd/asZliszSqHKr15tgRKRolp7rfKWqtWLPb1ZUWwur5bjIs0xoVg7Vdibfez4CU6ZPq2hCb+0f7tP3pfWq+c710vf00x6et19WuznBe7+6rxWeAHc/W53v3pgINfZ1z1JLVAinaMZx2itvVaaEC9/zOzcYPmht5jZb4Y/1V6n2FxeLcdFmmMijz1V7bZmeA6n9RUYf/4Ep/cVGqqsaP92n7wvrVfPd07f09p7ek8BfgNYBLwg+pi7vznLgmVN44byrZdaoEQ6UR6O0bwtYdTrPb1mtgr4B+AR4JXAw8B84Ivu/to021Bsbkyzj4m8ZVdkKdx3wxrXKw1oV1zq5mOzUVktWfQJYC3wDPDt2I9I3dQCJZJvzT5G0/Qk19ozrAySprsBeIu7XwIcCX5fDeQ2+6vbNHtZr43b9nBwbJybt+1pyvbbKdx3O0YPcXBsnI3b9uh8IWWlHUeb9nWNUuZT7Wqt9C4DXuXu17r7n0Z/mlG4LARpV7ePjY21uyhSQfSkoQtVkfxp9sV1MwK4Lgqa7iXu/unYfX8HVM38UmzOXjNjZ2PrfORb2KAH6HwhZdUbT5oVh2ptiNa1de2V3u8Cfc0oSLNo3FDn0YWqSL41I3g2oydZGSRN96SZnRv8/R0z+yngpcD0ai9UbE6nlmOtWuys57hdt3weswb6Wbd8XurXdJqwQS/8rNHzRa37TBWL7lVvPGlWHNKcGLWrJ735LjN7o5ktjf40o3DSmzTBlUi+NSN4hgEcSHXMNyMdWmr2N8Crg7//Evg3YBfwobaVqMvUcqxVu7iu57jttWMo3qNd6z5TxaJ71XssVHpdK69x1Qhce6X37cC5wI3AxyI/H824XNLDoicIBRCR/MnD2r06N7Sfu9/s7v8c/P0J4OXAT7j7H7e3ZN2jlmMtfnEdv6DWRW9lSeeU+D6rVknRPu49jVRcWxnHeq0BK0mtSxZpIXppKQUQkfyptVe2FmmPeY1naj8ze62ZDQV/nwe8D/ifwd+SgUYuVOMX1LrorSzpnBLfZ9UqKdrHvaeRiquucVurpiWLOpGZrQBWzJ07922PPPJIu4sjItI18rCMUVrNKKuWLLJvAsvc/btmtim4+0fA2e6+ssprFZubrNJSKvUsd9JLS6SU+6xJ+7SX9otM1cwli5rx3eqF72tWSxZhZr9gZneY2d3B7cV5HtOryTJERJojD8sYpX2eWtSb4vygwluguLrD1cDvAq+q9kLF5trUk6lQqdexnt6pXhpSUO6zJu3TpOcqs6R3NLN3X6saZKumSq+ZvQP4a+C/gJ8N7v4RxbX6RESkhzR7qbEsx/cq7bApDgezN/8c8A13fza4f0Yby9SVGr1QzWJ8by81HNXyWZOe28sVC8mu0UOrGmSrpvRmM/s2cJm7f8fMnnL3F5rZdOBJd39R00qZgcWLF/sDDzzQ7mKIiHSlZqQPp00ba2Z6WSVKb7Zrgd8HTgHWuvtmM3stsNHdh9NsQ7E5nUa/4500FCFv6k0Fb8c5SfKhVcfb2s072brrACsXzuaW1Zc07X06TbnYXGul90lglrsfN7ND7j5oZv3AqLvPyrC8mVNgFRFpnl68yOv1Si+Amb0cOO7u347c7nP3r6d5vWJza/Ti8ZkVNRhIrVp1vF20/l5OOEwzePSmy5v2Pp0mqzG9XwDWxe67huLafCIi0qPykj6ssXTNZWY/G/l7KXABMMfMlkZun92u8kmyassZSXmNpoNqX/eeVsXDlQtnM82Kv6PWbt7JRevvZe3mnU19/05Ta0/vLOBu4CzgfOBR4DCwwt0PNqWEGVFrsohI92tlr0wv9vSa2W53nx/8PVrmaZ52KUPF5vZo5Djphdlfs6SeYgm16thJ6gEO33vJ0CA7Rg919fGbSU+vuz8O/CTwy8Aa4NeB4TxXeM1shZndPjY21u6iiIhIk/XyJB2t4O7zzWxa8PdQmZ+qFV7F5vpl0XPYyHGiSZpqo3OShJpx7CSdD5J6gMP33rrrQM8ev7X29L7L3f8s4f4/cPe/yLRkGVNrsoiIZKkXe3oBzOxx4O+BT7j77ka2pdhcu6x7DmvtfdL4YJH6NOPYqXY+iPbwjoweYjjye/ueJwFYt3xeVx3LWY3pvb7M/dfVXiQRERHpQL8DDAFfMbOvmdk7zUzjeFsk657DWnuf8jJ+X6TTJB07jWZuVDsfhMf3yOghvrT+Mm5ZfQlfWn8ZO0YPcXh8gsPjEz3T65uq0huZoKJgZq8Nbwc/bwWeaW4xRUREJA/c/S53/yVgFvAR4JeAx8xsq5n9dzPTOr1NlHWls5FKtCZpSk/7SpI0mvJc7XxQ7vi+ZulcZvYXGOgvlB7r9u9oqvTmyGQVLwYeizzkwEHgJne/O/viZUcpVCIikqVeTW9OYmZDwJuAtwKnuftZaV6n2NyYZkyMU8s2NUlTetpXkqTVwwUqHd9J39FOnLiuofTmcHIK4Fh8sgp3fxXwsqwLLCIiIvlnZn0UJ7kcBs4FUq3RK41rxsQ4tWxTkzSlp30lSbJKeU77mkrHd9J3tJsmrqt1TO9Emfs1pldERKSHmNmrzex24AngBmAH8HJ3f217S9Y7mlGRqmWbGt+bnvaVpFVPRTPta8Lje3hokAUb7mPBhvtKFeWk72jS+aBT06BrHdM7XWN6RUREepeZbTCzbwPhsKbL3f3l7v4+d++sq6AO14yKVHSbtVzcduqFcDtoX0kl9VQ00zZWhcd3dCKrm7ftqfr86DmmU3t/0/b0fiz46QfuiNz+KPBbwDuaUjoRERHJmyXAe4BZ7n61u/9nuwskzVHLxW2nXgi3g/aVlBOOoY2P8a32nam1AeyaSOU4/eK1J1/bian6tY7p/WTCmN6fAu5pbjHrZ2YrzOz2sbGxdhdFRESk47n76919s7uP17sNxeZsNavnMH5xW+l9OvVCuB20r6SccpXbNL2/tZwH1gzP4car5jNroJ91y+fV9PpOTdVPNXtz2Reb/TjwZuBX3X12ZqVqAs0QKSIiWdLszY1TbM5Gq2YGTvM+nTjbaztpf0lULbM5x4/HRs8Djbw+T9/jhmZvjm3o7GAh+q8BDwKXAu9svIgiIiIiUqtW9RymeR+l7tZG+0uiaulFjR+PS4YGMeCpI0frmvl5ydBg6syOuE74HqedyGpGsOD83cD3gN8G7gSeBn7J3T/dvCKKiIiISDmtSjdMM8mVUndro/0l9Yof99v3PIkD4xOeuvK5aWQf123ZzcGxcUZGD03aXrctX5a2p/cJ4CPAt4Al7n6xu78PONa0komIiIhITVo1M3C5C+J6Z3/uVeH+ArSvJBMGqSuft27fywmHaTb1Nd22fFnaSu9DwJkUF57/STN7YdNKJCIiIiJ1aVWaYZpJrjoh5TEvtK+kUeuWz2PWQD/vv2o+kK4RJTyOb1g1f0qFtRMqsrVIO3vza4CXAp8D3gUcDFKdTwdmNK10IiIiIpJaq9IM4xfESZW2Tkh5zAvtK2lU9JhM24iStmLbDVkbqSeycvd9wcLzLwMuAx4HTgC7zOx/NauAIiIiIpJOu3pnopW28AIZ6KqeomZSWrhkZdPIPp4+chQDhocGM9lmN2Qi1Dx7M4C7f9HdrwbOA94B/HimpRIRERGRjlFPL5Mk0/6TRmzctofxCceBkdFDQPWe2mqPd0MmQl2V3pC7j7v7p9x9eVYFEhEREZHGtavHMM14XymvGyoY0n7RCa2qNaRUe7wbxvc2VOkVERERkXxqV49hmvG+Up5mdJZGRCe0Co/Bag0pvdDQ0pGVXjN7hZl92Mw+Y2a/2+7yiIiIiORNXi5k81KOTqPGAqlH0vhwqDy+vtae3E7M3mh5pdfM7jCzJ81sd+z+15vZt8xsr5mtq7QNd/+mu/8O8MvA4maWV0RERKQT5WVypG5IjWyHpMnBOqmSIe23cdseDo6Ns3HbntJ9WXyXOrFBph09vR8HXh+9w8ymAx8ElgMXA280s4vN7MfN7J7YzznBa1YCXwTub23xRURERDpLJ16k9jpNDiaNOnb8xKTfkM25oBOzN1pe6XX3LwCHYndfCux190fd/RiwGbjS3b/u7lfEfp4MtrPV3V8F/GprP4GIiIhIZ+nEi1Q5Sf+/3pNFj+wp06dN+g3ZfJc6MXuj0O4CBM4HHovc3g8Ml3uymb0G+G9AH/DZCs+7Grga4CUveUkGxRQREZFGKDa3x5rhOR11gSqT6f/Xe6I9svX+79ctn8dt2/cyPDTIgg33le4LJ0rrJXmZyMoS7vNyT3b3z7v7Ne7+2+7+wQrPu93dF7v74rPPPjuTgoqIiEj9FJvbS2NDRTpDlj2yO0YPcXh8gsPjE01Jke+E80peKr37gRdHbl8AHMhiw2a2wsxuHxsby2JzIiIi0iDF5vbR2FCRzpBlCvE1S+fSXzAMGB4aLN2fVWW1E84rean0fgV4mZkNmdkpwGpgaxYbdve73f3qgYGBLDYnIiIiDVJsbh+NDRXpPWuG53Dm6X04sH3Pk6WKblaV1U44r7R8TK+ZfQp4DXCWme0H/sTdP2ZmbwfuA6YDd7j7w60um4iIiEg309hQkd50zdK53LZ9L0eOTpQquuF9jVZWO+G8Yu5lh852BTNbAayYO3fu2x555JF2F0dERLqEmX3V3bVWfB0Um/Mh7Om5Zunc3F+wikg2No3sK1V0u/G4Lxeb85Le3DRKoRIREckXxeZ86IRxeCKSvaQuz06YjKoRXV/pFREREZGpOmEcnohkq1xjV7c3gnV9pVczRIqIiOSLYnM+RGeH7fZeHpFOlfWxWa6xq9sbwbp+TG9o8eLF/sADD7S7GCIi0iU0prdxis35seSm+zk4Ns6sgX6+tP6ydhdHRALNPDYrjevv1DH/PTumV0REREQq6/ZeHpFO1cxjs1JKc7elO3d9pVcpVCIiIvmi2Jw/YaozoDRnkRyJDkPIWqUKdbc1hCm9WUREpA5Kb26cYnP+KM1ZROI6KdVZ6c0iIiIiUlG39e6ISONqSXXO66R4qvSKiIiICKAZnUVkqloaw/I6FrjrK70aNyQiIpIvis2dIa8XryLSWrWMK85rtkjXV3rd/W53v3pgYKDdRREREREUmztFXi9eRXpVJ2RfNHPirUZ0faVXRERERGqnVGeRfFH2Rf1U6RURERGRinSxLdJ+zcq+6IVGLVV6RURERKSi6MV2L1wgi+RRs1KHe6FRq+srvZosQ0REJF8UmztP9GK7Fy6QRXpJLT3Indro1fWVXk2WISIiki+KzZ1NE1yJtFfWFc9aepA7tdGr6yu9IiIiIpKdvM7OKtIr2lnxrLXRKy89w6r0ioiIiIiIdIh2ZlvU2uiVl55hVXpFREREREQ6RDOzLbLumc3LcIhCW99dREREREREciHaM5tFpXrN8JxcDIXo+p5ezRApIiKSL4rNIiKNa8Z42bz0zGbN3L3dZWiJxYsX+wMPPNDuYoiISJcws6+6++J2l6OTKTaLiNRvyU33c3BsnFkD/Xxp/WXtLk4ulIvNXd/TKyIiIiIi0m26tVe2GVTpFRERERER6TDNXj4sTfp0rSnW7VrCSJVeERERERGRDtaMymSa5YZqXZKoXUsYqdIrIiIiIiLSwZpRmUyTPl1rinW7UrI1kZWIiEgdNJFV4xSbRUSysWlkH7dt38s7ls7NxRJB7VIuNmudXhERERERkQ6Wl/Vw80rpzSIiIiIiIl2gGWN72zX5VJa6vtJrZivM7PaxsbF2F0VERERQbBYRaZZmjO1Nu808z+Tc9ZVed7/b3a8eGBhod1FEREQExWYRkWZpxkRRabeZ55mcu77SKyIiIiIi0guia/dm1ZOadpt5nslZE1mJiIiIiIh0mWhPalaTXFXaZq2TabVy8i319IqIiIiIiHSZdqY65416ekVERERERLpMM3pSO3VpJPX0ioiIiIiISGbytsyRKr0iIiIiIiKSmVbOzJyGKr0iIiIiIiKSmbyN/dWYXhEREREREclM3sb+qqdXREREREREmiIP43tV6RUREREREZGmyMP43o6t9JrZ6Wb2VTO7ot1lERERERERkanyML635WN6zewO4ArgSXefH7n/9cBfAdOBj7r7xiqbuhb4p6YVVERERERERBqSh/G97ZjI6uPAB4BPhHeY2XTgg8AvAPuBr5jZVooV4Jtir/9NYAHwDaC/BeUVERERERGRDtXySq+7f8HMLozdfSmw190fBTCzzcCV7n4TxV7hSczstcDpwMXAj8zss+5+IuF5VwNXA7zkJS/J9HOIiIhI7RSbRUSk1fIypvd84LHI7f3BfYnc/T3uvhbYBPxNUoU3eN7t7r7Y3RefffbZWZZXRERE6qDYLCIirZaXdXot4T6v9iJ3/3j2RREREREREZFukZee3v3AiyO3LwAOZLFhM1thZrePjY1lsTkRERFpkGKziIi0Ul4qvV8BXmZmQ2Z2CrAa2JrFht39bne/emBgIIvNiYiISIMUm0VEpJVaXuk1s08BXwJ+zMz2m9lvufsE8HbgPuCbwD+5+8MZvZ9ak0VERHJEsVlERFrJ3KsOne0KZvZ9YF+7y5FDZwE/aHchckj7JZn2y1TaJ8l6Yb/McXfNxNQAxeayeuH4qYf2SzLtl6m0T5L1wn5JjM09U+mVZGb2gLsvbnc58kb7JZn2y1TaJ8m0X0Tqp+MnmfZLMu2XqbRPkvXyfsnLmF4RERERERGRzKnSKyIiIiIiIl1LlV65vd0FyCntl2TaL1NpnyTTfhGpn46fZNovybRfptI+Sdaz+0VjekVERERERKRrqadXREREREREupYqvV3MzF5sZv9mZt80s4fN7J3B/YNm9n/M7JHg9wsjr1lvZnvN7Ftmtqx9pW8+M5tuZjvN7J7gds/vFzM708w+Y2Z7gu/NT2m/gJn9j+AY2m1mnzKz/l7cL2Z2h5k9aWa7I/fVvB/M7CfM7OvBY7eambX6s4i0i2JzZYrNUyk2J1NsLlJsTkeV3u42Afyhu78CWAL8vpldDKwD7nf3lwH3B7cJHlsNvBJ4PfAhM5velpK3xjuBb0Zua7/AXwH/6u7zgIUU909P7xczOx+4Bljs7vOB6RQ/dy/ul49T/ExR9eyHvwauBl4W/MS3KdLNFJsrU2yeSrE5RrF5ko+j2FyVKr1dzN0fd/evBX8/Q/EkeT5wJfB3wdP+DlgV/H0lsNndj7r7KLAXuLSlhW4RM7sAuBz4aOTunt4vZjYT+FngYwDufszdn6bH90ugAJxqZgXgNOAAPbhf3P0LwKHY3TXtBzObBcx09y95cVKJT0ReI9L1FJvLU2yeSrG5IsVmFJvTUqW3R5jZhcAlwAhwrrs/DsXgC5wTPO184LHIy/YH93WjW4A/Ak5E7uv1/XIR8H3gb4PUso+a2en0+H5x9+8BfwZ8F3gcGHP3z9Hj+yWi1v1wfvB3/H6RnqPYPMUtKDbHKTYnUGyuSrE5RpXeHmBmLwD+GVjr7ocrPTXhvq6b3tvMrgCedPevpn1Jwn1dt18otpj+f8Bfu/slwBGCdJgyemK/BONgrgSGgNnA6Wb2a5VeknBf1+2XFMrtB+0fERSb4xSby1JsTqDYXLeejc2q9HY5M5tBMah+0t3/Jbj7iSCNgeD3k8H9+4EXR15+AcVUkW7z08BKM/sOsBlYamb/gPbLfmC/u48Etz9DMdD2+n75eWDU3b/v7s8D/wK8Cu2XUK37YX/wd/x+kZ6h2JxIsTmZYnMyxebKFJtjVOntYsGsax8DvunufxF5aCvw68Hfvw7cFbl/tZn1mdkQxUHsX25VeVvF3de7+wXufiHFwfzb3f3X0H45CDxmZj8W3HUZ8A16fL9QTJ1aYmanBcfUZRTH4PX6fgnVtB+CNKtnzGxJsD/fHHmNSNdTbE6m2JxMsbksxebKFJvj3F0/XfoDvJpiasJDwIPBzxuAF1Gcye2R4Pdg5DXvAb4NfAtY3u7P0IJ99BrgnuDvnt8vwCLggeA7swV4ofaLA/wpsAfYDfw90NeL+wX4FMWxU89TbBX+rXr2A7A42JffBj4AWLs/m37006ofxeZU+0ixefL+UGxO3i+Kza7YnPbHgg8pIiIiIiIi0nWU3iwiIiIiIiJdS5VeERERERER6Vqq9IqIiIiIiEjXUqVXREREREREupYqvSIiIiIiItK1VOkVAczsw2b2x+0uR73M7HfN7Akze9bMXtTu8uSBmb3OzLY0adt/YWa/04xti4hIkWJz91FslnbRkkXS9czsO8C5wARwnOKi7p8Abnf3E3Vs663u/n8zLmbdzGwGcBhY4u672l2evDCzB4C3u/uOJmx7FsVF7V/q7sey3r6ISLdTbO5Nis3SLurplV6xwt3PAOYAG4FrgY+1t0iZORfoBx5OetDMCq0tTvuZ2U8CA80IqgDu/jiwB1jZjO2LiPQIxeYeotgs7aRKr/QUdx9z963ArwC/bmbzAczs42Z2Q/D3WWZ2j5k9bWaHzOw/zGyamf098BLg7iBV6Y+C53/azA6a2ZiZfcHMXhm+X7DdD5rZvWb2jJmNmNlLI4+/0sz+T/A+T5jZu4P7p5nZOjP7tpn90Mz+ycwG45/HzF4OfCu4+bSZbQ/udzP7fTN7BHgkuO9tZrY3eK+tZjY7sh03s98zs0eCcr7PzF5qZl8ys8PB+5+StE/NbK6Z/Xvw+X9gZv8YeexVZvaV4LGvmNmrIo993sxuMLP/F+zPu83sRWb2yeA9v2JmF0ae/1dm9ljw2FfN7Gcq/KuXA/8eee2FwWcsRO77vJm9Nfj7N8zsP83sL4P/+6NB2X8jeM8nzezXY+/xeeDyCmUQEZEUFJsVm2Pvr9gsmVOlV3qSu38Z2A8knZz/MHjsbIotte8uvsTfBHyXYsv0C9z9fwXP3wa8DDgH+Brwydj23gj8KfBCYC/wfgAzOwP4v8C/ArOBucD9wWuuAVYBPxc89hTwwYTP8V9AGMjPdPelkYdXAcPAxWa2FLgJ+GVgFrAP2Bzb3OuBnwCWAH8E3A78KvBiYH7wOZK8D/hc8PkuAG4LPt8gcC9wK/Ai4C+Ae23yuKbVwJuA84GXAl8C/hYYBL4J/EnkuV8BFgWPbQI+bWb9Zcr045y84EhrGHgoKOsmivvnJyn+X34N+ICZvSDy/G8CC2t8DxERKUOxWbE5gWKzZEKVXullByiepOOepxh85rj78+7+H15h8Lu73+Huz7j7UWADsNDMBiJP+Rd3/7K7T1AMuouC+68ADrr7n7v7eLCNkeCx3wbe4+77I9v9RastHeomdz/k7j+iGCDvcPevBdtbD/xUtLUWuNndD7v7w8Bu4HPu/qi7j1G8eLikzPs8TzE1bXbwOb4Y3H858Ii7/727T7j7pyimHa2IvPZv3f3bkff4trv/32BffTr6nu7+D+7+w2Bbfw70AT9WpkxnAs+k2ksnjbr737r7ceAfKV5QvNfdj7r754BjFINs6JngfUREJDuKzYrNUYrNkglVeqWXnQ8cSrj/f1Ns9f1ckEqzrtwGzGy6mW0MUp0OA98JHjor8rSDkb+fA8IWyRcD3y6z6TnAnUE6z9MUWy6PU2zdTuuxyN+zKbYgA+DuzwI/pLgPQk9E/v5Rwu1oS2rUHwEGfNnMHjaz30x6z8C+et/TzP7QzL4ZpGM9DQwweT9HPQWcUeaxcuLvjbtX2gdnAE/X+B4iIlKZYrNic5Ris2RClV7pSVacTOF84Ivxx4JW3T9094sotnz+gZldFj4ce/oa4Erg5yme6C8M3yJFMR6jmDZU7rHl7n5m5Kff3b+XYruljxL5+wDFYF0snNnpFFOFatle8pu4H3T3t7n7bIqt4B8ys7nx9wy8pJ73DMYIXUsxBeyF7n4mMEb5/fwQ8PLI7SPB79Mi951XazliXgFoRk4RkYwoNis2o9gsTaJKr/QUM5tpZldQHBPyD+7+9YTnXBFMAGEUlxs4HvxAscXxosjTzwCOUmyZPQ24sYbi3AOcZ2ZrzazPzM4ws+HgsQ8D7zezOUGZzjazK2vYdtwm4C1mtsjM+oJyjrj7dxrYJkHZfsnMLghuPkUxoB8HPgu83MzWmFnBzH4FuJji567VGRSXtfg+UDCz64GZFZ7/WYpjrgBw9+9TDOi/FvQA/CblL2rS+jmKaV8iItIAxWbFZsVmaTZVeqVX3G1mz1BspX0PxYkb3lLmuS+jOInFsxQnb/iQu38+eOwm4LogteldFNcU3EfxpP0NIPU0/O7+DPALFFusD1KcyfG1wcN/BWylmMb1TLDd4aTtpHyv+4E/Bv4ZeJxiUFld7/ZifhIYMbNnKZb5ne4+6u4/pDg26g8pXnj8EXCFu/+gjve4j2IQ+y+K+3ucySlik7j714CxyIUKwNuA/xmU5ZXA/6ujHEBpLcCLgS31bkNERBSbUWxWbJaWsApzAIiIdCwzex3we+6+qgnb/nOKE3t8KOtti4iIdCvFZmkXVXpFRERERESkaym9WURERERERLqWKr0iIiIiIiLStVTpFRERERERka6lSq+IiIiIiIh0LVV6RUREREREpGup0isiIiIiIiJdS5VeERERERER6Vqq9IqIiIiIiEjX+v8BbR1c528+wMkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1152x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df_ex = df.query('synapse_reversal==0')\n",
    "fig, axs = subplots(1,2,figsize=(16,4),sharex=True,sharey=True)\n",
    "plt.sca(axs[0])\n",
    "plt.scatter(df_ex.seg_dist,df_ex.alpha,s=3)\n",
    "plt.ylabel(r'Attenuation ($\\alpha$)',fontsize=12)\n",
    "plt.xlabel('Distance from soma (um)',fontsize=12)\n",
    "plt.sca(axs[1])\n",
    "plt.scatter(df_ex.seg_dist,df_ex.alpharho,s=3)\n",
    "plt.ylabel(r'Visability ($\\Gamma$)',fontsize=12)\n",
    "plt.xlabel('Distance from soma (um)',fontsize=12)\n",
    "plt.yscale('log')\n",
    "plt.xlim(0.00001,plt.xlim()[1])\n",
    "plt.suptitle(cell_name)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.7"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}