{ "cells": [ { "cell_type": "markdown", "id": "fed475d53adb2402", "metadata": {}, "source": [ "# Simulating a Hodgkin-Huxley Type Thalamic Network with `braincell`\n", "\n", "[![Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/chaobrain/brain-modeling-ecosystem/blob/develop/docs/braincell_HH_thalamus_oscillations.ipynb)\n", "[![Open in Kaggle](https://kaggle.com/static/images/open-in-kaggle.svg)](https://kaggle.com/kernels/welcome?src=https://github.com/chaobrain/brain-modeling-ecosystem/blob/develop/docs/braincell_HH_thalamus_oscillations.ipynb)" ] }, { "cell_type": "markdown", "id": "3ad151974837e0a4", "metadata": {}, "source": [ "This section introduces the simulation of a Hodgkin-Huxley type thalamic network using `braincell`.\n", "\n", "As a practical example, before learning how to model a Hodgkin-Huxley type thalamic network with `braincell`, you should already be familiar with modeling a Hodgkin-Huxley type [E–I network](braincell_HH_network-zh.ipynb)." ] }, { "cell_type": "markdown", "id": "1c15ce72366e2270", "metadata": {}, "source": [ "## Overview of the Thalamic Network\n", "\n", "The thalamus plays a critical role in the generation of thalamocortical oscillations, yet the underlying mechanisms remain unclear. To explore whether the thalamus alone can generate multiple distinct oscillations, we construct a computational model.\n", "\n", "Thalamic activity is modulated by various neurotransmitters, among which acetylcholine (ACh) and norepinephrine (NE) are key regulators that shape oscillatory dynamics.\n", "\n", "The thalamic network model we build consists of two nuclei: the lateral geniculate nucleus (LGN) and the thalamic reticular nucleus (TRN). The LGN contains three types of neurons: high-threshold bursting thalamocortical cells (HTC), relay-mode thalamocortical cells (RTC), and local interneurons (IN). The TRN contains a single type of neuron: reticular cells (RE).\n", "\n", "These different neuronal populations are interconnected through gap junctions and synapses. Under the modulation of ACh and NE, the network generates oscillations across four frequency bands: delta, spindle, alpha, and gamma." ] }, { "cell_type": "markdown", "id": "b250343a22c8f7ef", "metadata": {}, "source": [ "## Thalamic Neuron Modeling\n", "\n", "After understanding the connectivity structure of the thalamic network, we first model the thalamic neurons." ] }, { "cell_type": "code", "execution_count": 172, "id": "ecb18514bc24f8a7", "metadata": { "ExecuteTime": { "end_time": "2025-09-22T03:31:50.867830Z", "start_time": "2025-09-22T03:31:44.712649Z" } }, "outputs": [], "source": [ "from typing import Dict, Callable\n", "\n", "import time\n", "import braincell\n", "import brainstate\n", "import brainpy.state\n", "import braintools\n", "import brainunit as u\n", "import matplotlib.pyplot as plt\n", "import numba\n", "import numpy as np\n", "\n", "\n", "brainstate.environ.set(dt=0.1 * u.ms)" ] }, { "cell_type": "markdown", "id": "7d64f41194442c5", "metadata": {}, "source": [ "First, we define a base class for thalamic neurons, which inherits from `SingleCompartment`." ] }, { "cell_type": "code", "execution_count": 173, "id": "9931eb2b65f67602", "metadata": { "ExecuteTime": { "end_time": "2025-09-22T03:31:50.879907Z", "start_time": "2025-09-22T03:31:50.876617Z" } }, "outputs": [], "source": [ "class ThalamusNeuron(braincell.SingleCompartment):\n", " def compute_derivative(self, I_ext=0. * u.nA):\n", " I_ext = self.sum_current_inputs(I_ext, self.V.value) * self.area\n", " for key, ch in self.nodes(braincell.IonChannel, allowed_hierarchy=(1, 1)).items():\n", " I_ext = I_ext + ch.current(self.V.value)\n", " self.V.derivative = I_ext / self.C\n", " for key, node in self.nodes(braincell.IonChannel, allowed_hierarchy=(1, 1)).items():\n", " node.compute_derivative(self.V.value)\n", "\n", " def step_run(self, t, inp):\n", " with brainstate.environ.context(t=t):\n", " self.update(inp)\n", " return self.V.value" ] }, { "cell_type": "markdown", "id": "aff18f6c12238204", "metadata": {}, "source": [ "In this example, the thalamic neuron model differs from the [HH neuron in the E–I network](braincell_HH_neuron-zh.ipynb).\n", "Here, the interface `compute_derivative` is explicitly defined to update the membrane potential dynamics.\n", "\n", "This interface integrates external input currents with the ionic currents from all channel components, and then updates both the membrane potential and the derivatives of the ion channels. This corresponds to the standard differential equation update method for Hodgkin–Huxley type neurons.\n", "\n", "Additionally, the method `step_run` is implemented to facilitate later simulation and visualization.\n" ] }, { "cell_type": "markdown", "id": "9aa545188c69fa84", "metadata": {}, "source": [ "Next, we model the four different types of thalamic neurons.\n", "\n", "Here, we take the HTC neuron as an example.\n", "\n", "The HTC neuron inherits from the `ThalamusNeuron` class defined above and updates its membrane potential dynamics following the same procedure as `ThalamusNeuron`.\n", "\n", "According to the reference, its membrane area is $2.9 \\times 10^{-4} \\ \\text{cm}^2$.\n", "\n", "In addition, the HTC neuron contains multiple types of ionic currents:\n", "\n", "- a spike-generating fast sodium current ($I_\\text{Na}$)\n", "- a delayed rectifier potassium current ($I_\\text{DR}$)\n", "- a hyperpolarization-activated cation current ($I_\\text{H}$)\n", "- a high-threshold L-type Ca$^{2+}$ current ($I_{\\text{Ca/L}}$)\n", "- a Ca$^{2+}$-activated nonselective cation current ($I_\\text{CAN}$)\n", "- a regular low-threshold T-type Ca$^{2+}$ current ($I_{\\text{Ca/T}}$)\n", "- a high-threshold T-type Ca$^{2+}$ current ($I_{\\text{Ca/HT}}$)\n", "- a Ca$^{2+}$-dependent potassium current ($I_\\text{AHP}$)\n", "- a potassium leak current ($I_\\text{KL}$)\n", "- a leakage current ($I_\\text{L}$)\n", "\n", "We then add these ionic channels to the HTC model and specify their parameters according to the reference." ] }, { "cell_type": "code", "execution_count": 174, "id": "30e3057d2cf7f5ce", "metadata": { "ExecuteTime": { "end_time": "2025-09-22T03:31:50.888772Z", "start_time": "2025-09-22T03:31:50.884072Z" } }, "outputs": [], "source": [ "class HTC(ThalamusNeuron):\n", " def __init__(\n", " self,\n", " size,\n", " gKL=0.01 * (u.mS / u.cm ** 2),\n", " V_initializer=braintools.init.Constant(-65. * u.mV),\n", " solver: str = 'ind_exp_euler'\n", " ):\n", " super().__init__(size, V_initializer=V_initializer, V_th=20. * u.mV, solver=solver)\n", "\n", " self.area = 1e-3 / (2.9e-4 * u.cm ** 2)\n", "\n", " self.na = braincell.ion.SodiumFixed(size, E=50. * u.mV)\n", " self.na.add(INa=braincell.channel.INa_Ba2002(size, V_sh=-30 * u.mV))\n", "\n", " self.k = braincell.ion.PotassiumFixed(size, E=-90. * u.mV)\n", " self.k.add(IKL=braincell.channel.IK_Leak(size, g_max=gKL))\n", " self.k.add(IDR=braincell.channel.IKDR_Ba2002(size, V_sh=-30. * u.mV, phi=0.25))\n", "\n", " self.ca = braincell.ion.CalciumDetailed(size, C_rest=5e-5 * u.mM, tau=10. * u.ms, d=0.5 * u.um)\n", " self.ca.add(ICaL=braincell.channel.ICaL_IS2008(size, g_max=0.5 * (u.mS / u.cm ** 2)))\n", " self.ca.add(ICaN=braincell.channel.ICaN_IS2008(size, g_max=0.5 * (u.mS / u.cm ** 2)))\n", " self.ca.add(ICaT=braincell.channel.ICaT_HM1992(size, g_max=2.1 * (u.mS / u.cm ** 2)))\n", " self.ca.add(ICaHT=braincell.channel.ICaHT_HM1992(size, g_max=3.0 * (u.mS / u.cm ** 2)))\n", "\n", " self.kca = braincell.MixIons(self.k, self.ca)\n", " self.kca.add(IAHP=braincell.channel.IAHP_De1994(size, g_max=0.3 * (u.mS / u.cm ** 2)))\n", "\n", " self.Ih = braincell.channel.Ih_HM1992(size, g_max=0.01 * (u.mS / u.cm ** 2), E=-43 * u.mV)\n", " self.IL = braincell.channel.IL(size, g_max=0.0075 * (u.mS / u.cm ** 2), E=-70 * u.mV)" ] }, { "cell_type": "markdown", "id": "7731be4399e4a6cf", "metadata": {}, "source": [ "Similarly, we model the RTC, IN, and TRN neurons.\n" ] }, { "cell_type": "code", "execution_count": 175, "id": "f888237abfc9a094", "metadata": { "ExecuteTime": { "end_time": "2025-09-22T03:31:50.897739Z", "start_time": "2025-09-22T03:31:50.893554Z" } }, "outputs": [], "source": [ "class RTC(ThalamusNeuron):\n", " def __init__(\n", " self,\n", " size,\n", " gKL=0.01 * (u.mS / u.cm ** 2),\n", " V_initializer=braintools.init.Constant(-65. * u.mV),\n", " solver: str = 'ind_exp_euler'\n", " ):\n", " super().__init__(size, V_initializer=V_initializer, V_th=20 * u.mV, solver=solver)\n", "\n", " self.area = 1e-3 / (2.9e-4 * u.cm ** 2)\n", "\n", " self.na = braincell.ion.SodiumFixed(size)\n", " self.na.add(INa=braincell.channel.INa_Ba2002(size, V_sh=-40 * u.mV))\n", "\n", " self.k = braincell.ion.PotassiumFixed(size, E=-90. * u.mV)\n", " self.k.add(IDR=braincell.channel.IKDR_Ba2002(size, V_sh=-40 * u.mV, phi=0.25))\n", " self.k.add(IKL=braincell.channel.IK_Leak(size, g_max=gKL))\n", "\n", " self.ca = braincell.ion.CalciumDetailed(size, C_rest=5e-5 * u.mM, tau=10. * u.ms, d=0.5 * u.um)\n", " self.ca.add(ICaL=braincell.channel.ICaL_IS2008(size, g_max=0.3 * (u.mS / u.cm ** 2)))\n", " self.ca.add(ICaN=braincell.channel.ICaN_IS2008(size, g_max=0.6 * (u.mS / u.cm ** 2)))\n", " self.ca.add(ICaT=braincell.channel.ICaT_HM1992(size, g_max=2.1 * (u.mS / u.cm ** 2)))\n", " self.ca.add(ICaHT=braincell.channel.ICaHT_HM1992(size, g_max=0.6 * (u.mS / u.cm ** 2)))\n", "\n", " self.kca = braincell.MixIons(self.k, self.ca)\n", " self.kca.add(IAHP=braincell.channel.IAHP_De1994(size, g_max=0.1 * (u.mS / u.cm ** 2)))\n", "\n", " self.Ih = braincell.channel.Ih_HM1992(size, g_max=0.01 * (u.mS / u.cm ** 2), E=-43 * u.mV)\n", " self.IL = braincell.channel.IL(size, g_max=0.0075 * (u.mS / u.cm ** 2), E=-70 * u.mV)" ] }, { "cell_type": "code", "execution_count": 176, "id": "b25d90796c2bfaf9", "metadata": { "ExecuteTime": { "end_time": "2025-09-22T03:31:51.180937Z", "start_time": "2025-09-22T03:31:51.176721Z" } }, "outputs": [], "source": [ "class IN(ThalamusNeuron):\n", " def __init__(\n", " self,\n", " size,\n", " gKL=0.01 * (u.mS / u.cm ** 2),\n", " V_initializer=braintools.init.Constant(-70. * u.mV),\n", " solver: str = 'ind_exp_euler'\n", " ):\n", " super().__init__(size, V_initializer=V_initializer, V_th=20. * u.mV, solver=solver)\n", "\n", " self.area = 1e-3 / (1.7e-4 * u.cm ** 2)\n", "\n", " self.na = braincell.ion.SodiumFixed(size)\n", " self.na.add(INa=braincell.channel.INa_Ba2002(size, V_sh=-30 * u.mV))\n", "\n", " self.k = braincell.ion.PotassiumFixed(size, E=-90. * u.mV)\n", " self.k.add(IDR=braincell.channel.IKDR_Ba2002(size, V_sh=-30 * u.mV, phi=0.25))\n", " self.k.add(IKL=braincell.channel.IK_Leak(size, g_max=gKL))\n", "\n", " self.ca = braincell.ion.CalciumDetailed(size, C_rest=5e-5 * u.mM, tau=10. * u.ms, d=0.5 * u.um)\n", " self.ca.add(ICaN=braincell.channel.ICaN_IS2008(size, g_max=0.1 * (u.mS / u.cm ** 2)))\n", " self.ca.add(ICaHT=braincell.channel.ICaHT_HM1992(size, g_max=2.5 * (u.mS / u.cm ** 2)))\n", "\n", " self.kca = braincell.MixIons(self.k, self.ca)\n", " self.kca.add(IAHP=braincell.channel.IAHP_De1994(size, g_max=0.2 * (u.mS / u.cm ** 2)))\n", "\n", " self.IL = braincell.channel.IL(size, g_max=0.0075 * (u.mS / u.cm ** 2), E=-60 * u.mV)\n", " self.Ih = braincell.channel.Ih_HM1992(size, g_max=0.05 * (u.mS / u.cm ** 2), E=-43 * u.mV)" ] }, { "cell_type": "code", "execution_count": 177, "id": "31357e54da148d96", "metadata": { "ExecuteTime": { "end_time": "2025-09-22T03:31:51.259583Z", "start_time": "2025-09-22T03:31:51.255745Z" } }, "outputs": [], "source": [ "class TRN(ThalamusNeuron):\n", " def __init__(\n", " self,\n", " size,\n", " gKL=0.01 * (u.mS / u.cm ** 2),\n", " V_initializer=braintools.init.Constant(-70. * u.mV),\n", " gl=0.0075,\n", " solver: str = 'ind_exp_euler'\n", " ):\n", " super().__init__(size, V_initializer=V_initializer, V_th=20. * u.mV, solver=solver)\n", "\n", " self.area = 1e-3 / (1.43e-4 * u.cm ** 2)\n", "\n", " self.na = braincell.ion.SodiumFixed(size)\n", " self.na.add(INa=braincell.channel.INa_Ba2002(size, V_sh=-40 * u.mV))\n", "\n", " self.k = braincell.ion.PotassiumFixed(size, E=-90. * u.mV)\n", " self.k.add(IDR=braincell.channel.IKDR_Ba2002(size, V_sh=-40 * u.mV))\n", " self.k.add(IKL=braincell.channel.IK_Leak(size, g_max=gKL))\n", "\n", " self.ca = braincell.ion.CalciumDetailed(size, C_rest=5e-5 * u.mM, tau=100. * u.ms, d=0.5 * u.um)\n", " self.ca.add(ICaN=braincell.channel.ICaN_IS2008(size, g_max=0.2 * (u.mS / u.cm ** 2)))\n", " self.ca.add(ICaT=braincell.channel.ICaT_HP1992(size, g_max=1.3 * (u.mS / u.cm ** 2)))\n", "\n", " self.kca = braincell.MixIons(self.k, self.ca)\n", " self.kca.add(IAHP=braincell.channel.IAHP_De1994(size, g_max=0.2 * (u.mS / u.cm ** 2)))\n", "\n", " # self.IL = braincell.channel.IL(size, g_max=0.01 * (u.mS / u.cm ** 2), E=-60 * u.mV)\n", " self.IL = braincell.channel.IL(size, g_max=gl * (u.mS / u.cm ** 2), E=-60 * u.mV)" ] }, { "cell_type": "markdown", "id": "173d5286ddc1241d", "metadata": {}, "source": [ "After completing the neuron modeling, we simulate the constructed single neuron.\n", "Here, we use the HTC neuron as an example.\n" ] }, { "cell_type": "code", "execution_count": 178, "id": "2128597ec244da7f", "metadata": { "ExecuteTime": { "end_time": "2025-09-22T06:44:38.202667Z", "start_time": "2025-09-22T06:44:36.484687Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Elapsed time: 1.4957 s\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAGdCAYAAAA8F1jjAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAATClJREFUeJzt3QmYXGWZ9/+7931Jp5PuhOwJEJKwhi3sQoagqKPwuoy8SiCCYPAvgg5EZlicJYw6QeRCxBGCL+KgKKCIIBh2CVvYkkA2IKSzdHfW7k7vS/2v+6k6Vaeqq7rqVHXVOV39/XiV1TunujupX+7nfu4nx+fz+QQAAMCDct2+AAAAgFgIKgAAwLMIKgAAwLMIKgAAwLMIKgAAwLMIKgAAwLMIKgAAwLMIKgAAwLPyZYQbGBiQnTt3SkVFheTk5Lh9OQAAIAE6b7atrU0mTpwoubm52RtUNKRMnjzZ7csAAABJaGhokEmTJmVvUNFKivVAKysr3b4cAACQgNbWVlNosJ7HszaoWMs9GlIIKgAAjCzx2jZopgUAAJ5FUAEAAJ5FUAEAAJ5FUAEAAJ5FUAEAAJ5FUAEAAJ5FUAEAAJ5FUAEAAJ5FUAEAAJ5FUAEAAJ5FUAEAAJ5FUAEAAJ5FUAGS8OS6XfL4u7vcvgwAyHoj/vRkINP2tffIFb9+07x8zJSz5ZDqErcvCQCyFhUVwKEPdx8Mvrxtb4er1wIA2Y6gAjj0sS2cdPb2uXotAJDtCCqAQ1v3tgdf7uodcPVaACDbEVQAh7Y0h5Z+unr7JZts398hT65rlIEBn9uXAgAGQQVwqK2rLy0VlQMdPbLjQKe4RcPJ1+55Ta749Rp5+K0drl0HANgRVACHfOJLS0XlCz9fLafe+oxsbmoTN+xt75EP9/iXtdZ8vN+VawCASAQVwCGfbVWkq2/4gsrmwJLSX9Y2iht8tgfW2tXryjUAQCSCCpBKUBmmpR97SGjvcWcnkb0rpbMnu3pvAIxcBBUgBX39wxNU+m3Nqx0uBRW79m73rwEAFEEFSKFHxR4wUtHbbw8q/a5XijqzbDcTgJGLoAKk8ITeN0xBpcdWmXFry7M9gFFRAeAVBBXAIXs0iVVR0dH6y594X5pbuxwvIQ1XlSaligo9KgA8gkMJgRT0xuhRuXjla/LRnnZ5/aN98vA3T437dezRpM+2DJRJvhgVHgBwExUVIIVn9FjVDw0p6s1tBxx/eS+EBHvPDAC4iaACpNDLMVw9KvZll54+d4KKfYv0cO1mAoBUEVSAFEJFOvpJ3Kqo2B9XL2f9APAIggqQgmGrqNiqNG5VVOyoqADwCoIKkNKun2F6QrdXMzwQEjR/cYIyAC8gqAAp9XKkYenHtR6V9FSLACAVBBUgDXNUUvmabu24sS8/qb7hqhYBQAoIKkAK0lF16PZIRYUtygC8gKACeGDXj/1rDkQmhgyJ/K/SUAvACwgqgEO+NCyPhM1m8UhAoEcFgBcQVACnbBWPWBWVnJzkv7xb+cDeJKwIKgC8gKACpCDWk3l+bk4KJzK71KMS8bpXKjsARjeCCpDCE3qsokN+bvJ/tNzabEMzLQAvIqgAKTyhRy6XJF1Rsb3c71IzbWRNJdZjA4BMIqgAKTS+xupRyctzuvQT/jW9EBLcC0wAEEJQAVIQa+knN5VuWpcaaiNzSToOXAQAzwaVW2+9VXJycuTqq68Ovq2rq0uWLl0qY8eOlfLycrnwwgulqakpU5cEpD7zJMaTucOVnyjj6zPfqBL5SCioABg1QeX111+Xu+++W4466qiwt3/nO9+Rxx57TB566CF5/vnnZefOnXLBBRdk4pKAtA5n01CeCjcaaqmoABiVQeXgwYNy0UUXyf/8z//ImDFjgm9vaWmRe+65R1asWCFnn322zJ8/X1auXCkvv/yyvPLKK+m+LCBpiTS+Oq2oRPLCOTtuTcgFgIwGFV3aOf/882XhwoVhb1+zZo309vaGvX327NkyZcoUWb16dcyv193dLa2trWE3wC2+YepRifw6rlRUIhZ/CCoAvCA/nV/8wQcflDfffNMs/URqbGyUwsJCqa6uDnt7XV2deV8sy5cvl1tuuSUt1wskwr4jJ+bSzwjccTN46SfjlwAAmauoNDQ0yLe//W154IEHpLi4eNi+7rJly8yykXXT/w7gltgj9J3OUfG530wbWdWhogIgm4OKLu00NzfLcccdJ/n5+eamDbM//elPzctaOenp6ZEDBw6EfZ7u+qmvr4/5dYuKiqSysjLsBrg38C36x6TYS+uNpR+aaQFk89LPOeecI2vXrg172yWXXGL6UK677jqZPHmyFBQUyKpVq8y2ZLVx40bZtm2bLFiwIF2XBWRk4FuqPSreaKZ1+woAII1BpaKiQubNmxf2trKyMjMzxXr7kiVL5JprrpGamhpTGfnWt75lQsrJJ5+crssChlWs5RHHc1S8uD2ZpR8A2d5MG89tt90mubm5pqKiu3kWLVokP/vZz9y8JMDhHJX0TKb1REWFkgqA0RZUnnvuubDXtcn2zjvvNDdgZJ6eHKuZ1vbxPl/c5trIs33caGSlmRaAF3HWD5CO7cm2YJJMYcILxQwm0wLwAoIKkILYzbTxP2bIHhU3KiqDBr5l/BIAYBCCCuCQz+Fk2kRChxfO2WHpB4AXEVQAp3zOAkUyT/huZITI/yRLPwC8gKACpKGZNryi4vSrutVM6/41AEAkggqQgphLP7mpPeF7oZhBUAHgBQQVIIXKQ38iFZVEmmk90B/ihaFzABCJoAKkY45K2Mc4+5rmdQ/MUWEyLQAvIKgAKR5KGDVU2CoqyTSlurP0w6GEALyHoAKkKOrzuS28+Ebs9uSMXwIADEJQAVIejDb4Gd2X4hKKFxpZWfoB4AUEFcChRKofiRxcOFT48cIcFTf6ZAAgEkEFcCjy+Tve8/mI2fXjgeUnAIhEUAHSsERir5CMlDkqkRUUggoALyCoACmK2qPidOnHAxWVSB64BAAgqAApj5qPk0SSqUy4Mkcl4nWaaQF4AUEFSHWCqy/erJUEelQGzTCRjPNiVQcACCpAioZje7IXpsIODksEFQDuI6gAqVYeom5P9qVUHXFla/CgXT+ZvwQAiERQAVIe+Db0x4+UXT+Dr8EDFwFg1COoAA45XaZJLqi430xLUAHgBQQVIEXRl36c7frxwjk7NNMC8CKCCpDyqPmhPz65EfruN9PSowLACwgqQBqWfuxP+smEDi9UM7xwDQBAUAEcS+D05BSXftyoZiSymwkAMo2gAqSjR8X+/iSe773RTJvxSwCAQQgqgENOG18TCR2D+17cOD05fqUIADKNoAKkYRtv2MC3hCbTOpvNkg5eCEsAEImgAjgU+QQerQclbIR+EqnDC9UMDiUE4AUEFSBF8Z7PfSOlP8QDs1wAIBJBBXDIl0jlwZfa7hkvzFFh6QeAFxBUgDRMcHW662fw9mQ3mmnDX0/mMEUAGG4EFcChQY2v8U5PTqgy4X4zbSR6VAB4AUEFSFG8UJHMEoo725O919ALAAQVwCHr6Ts3J31LP14Y+EZOAeAFBBXAqcATeH5ubkKnJycTOtw5PTn+tmsAyDSCCpBsRSU30cm0iX9NL1VUWPoB4AUEFSBJeTk5w3Z6cuSHeCEjeOEaAICgAjhkBY/cQJNKOpZ+PLE9maQCwAMIKoBD1tN3vhVUop7142weiTcOBKRHBYD3EFQAh6wMkRcMKsM/j8SdZlr3rwEAIhFUAIes/pNcq0clzjN6Qj0qSXzOcPPCNQBAJIIKkKShl37sk2lHZn+IF64BAAgqgEPW83ewmTbOE/rImaMS/no/OQWABxBUgBSbaaMt/TieTDvorB9OTwaAtAeV5cuXywknnCAVFRUyfvx4+dznPicbN24M+5iuri5ZunSpjB07VsrLy+XCCy+UpqYmfjrwroiKSrTnc/vbknnCj7blOd28uPwEAGkNKs8//7wJIa+88oo8/fTT0tvbK+eee660t7cHP+Y73/mOPPbYY/LQQw+Zj9+5c6dccMEF6bwsYHgHvkWtqNh6VJIYTeuFHTeJbKsGgHTLT+cXf/LJJ8Nev++++0xlZc2aNXLGGWdIS0uL3HPPPfKb3/xGzj77bPMxK1eulCOOOMKEm5NPPjmdlwckxQoh1vbkeNuPR+oI/WS2VQPAiO5R0WCiampqzL0GFq2yLFy4MPgxs2fPlilTpsjq1aujfo3u7m5pbW0NuwFuzlGJtrST6mRaNzJC5OOgRwXAqAoqAwMDcvXVV8upp54q8+bNM29rbGyUwsJCqa6uDvvYuro6875YfS9VVVXB2+TJkzNy/UDsZtrYH2NeHrHbk92+AgDIYFDRXpV169bJgw8+mNLXWbZsmanMWLeGhoZhu0bAiUS3Jyc3mZZmWgBIe4+K5aqrrpI///nP8sILL8ikSZOCb6+vr5eenh45cOBAWFVFd/3o+6IpKioyN8At1pKI1Uwb96yfhCbTRm5PFte5sfMIADJaUdG/0DWkPPLII/LMM8/I9OnTw94/f/58KSgokFWrVgXfptuXt23bJgsWLEjnpQFJs56+g820UZ/QfSkt/bgzQt97YQkA8tO93KM7ev74xz+aWSpW34n2lpSUlJj7JUuWyDXXXGMabCsrK+Vb3/qWCSns+MFIPpQw/PRk58/4bpxczNIPgFEXVO666y5zf9ZZZ4W9XbcgL1682Lx82223SW5urhn0pjt6Fi1aJD/72c/SeVnAsAgGlYF0bE8W90foU1IBkO1BJZHydXFxsdx5553mBowkwdOTo/WopHzWj8/VAKYhhYIKAC/grB8gyfCd6OnJiQT2wTNMxL3emyGahAEg0wgqgAP25+70Lv240aPi/2/mBv5WIKgA8AKCCpCkIZtpU176ERcH2fn/WqBFBYAXEFQAB3xRBr5FPZQwbI6Kwy/sVjXDOhXa/7CoqADwBIIKkIYeFUm1ouJiOSNUKSKoAHAfQQVwwP7UPfRkWp+j0BE5bM2VOSqBa7B2Mw1EOcMIADKNoAIkKbT0E69HJYnJtOKFQXZUVAC4j6ACJLvrJ1B5iLr92OFZP15Y+ok8GoCgAsALCCqAA/Ylmry8oc76EYdzVMJfdzMkBJd+yCkAPICgAjhgzw9WM238ybTO/zv9Lo7Qt+aouHEwIgBEIqgAKVYeoq782JtpE6moDPH5ma4WWXNUOOsHgBcQVIAkWb0c8Z7QE2um9cCun0FzVDJ+CQAwCEEFSHGEfryln2SqI66M0A/c00wLwEsIKkCyzbSBJ/RoQSR8Mq3zpR83Z5gMtaQFAJlGUAGSZG1PjjpC3xY9os1ZicedEfo+R0taAJAJBBXAAXt+sAa+xXs+T2Z7crTlpHRj6QeAFxFUAAfsT93Bs37iHkro/FRCV05PDjbTsvQDwDsIKoAD9upIos20yYQOVybTRiz9UFEB4AUEFSBJoSf0KO902kzrocm0wd4bggoADyCoAKmenhynmTaZ53t3Tk+OnEzLdFoA7iOoACk30w5+MrfnjOQm04rrpye7dR0AYEdQAZyIdtZP1GZahyP0PbTrx2qmdes6AMCOoAIkaaiKSrLNtKHx9S72qNgqKjTUAnAbQQVIdjKt1aMS9VBC+8uJLP1E7Lhxc9ePraJCTgHgNoIK4IAvgaWfVA4YtJZd3BwKa1WK/NdBUgHgLoIK4IAvgaWfwVuNnQ9bc/P0ZCuAuXUdAGBHUAGSlBf40xMZVOK9Ho31EVZIcGNbsLX8FF5RyfhlAEAYggqQ9GTa3OhLP4M+J/Gvm5fn/rC18B4VkgoAdxFUgFQHvsWZKuukz2OoBt10Cy0/hd5GRQWA2wgqQLID36ytxIOaacVxULE+xtVdP4F7+9IPPSoA3EZQAZLo49DCh/WEHm+ZxkkzbfBEZjd6VAL/yRzJCYYwln4AuI2gAiQhx8HST0JzVKzx9VaPiouVDBPCPLBNGgAUQQVwIuysn0SXfhL/slb4ceWsH9uDCwUVkgoAdxFUAAesp+2cHF0eiTHwLeJzEjvrJ3xrcL+rSz/+qoq5DkoqAFxGUAGSfDIvDAxS6RsYGHrXz4DzOSpuVjI0pFhNvRRUALiNoAIk+WSeHwgqvf2pL/1YSSW45BKefTLCqur4m2ndD0wAoAgqQJJ9HPmBxtfe/ohUkcT2ZOvrWl/T1V0/OaGlH4IKALcRVIAkt/AWBLpp+yIqKoMHviX+dUM7idyboxK+64egAsBdBBXAgeDTdo5IQX70isrgEfqJ96gEB775Mj/DJPSfywm7DgBwE0EFSKqPQxtfrR6VOM20Tuao2KbCZrqYETyU0FRU/G+jogLAbQQVIMk+joJAP0nfsMxRCR+h78YW5fAeFfcHzwGAIqgASdAeFWvXT2SPir3hNpWKSqarGeG7fsKvCwDcQlABkjxh2Kqo9PQPhPWTRD65J/JkH+pRCf2RzPQWZXszrZtNvQBgR1ABHLCeuHVpxNr1E7lEEvncntDyiXV6cqig4kJFxX+vO36spR9WfgC4jaACJFN1sM08iexTSWrpJ0pFJeM9KvazfgKXQY8KALd5IqjceeedMm3aNCkuLpaTTjpJXnvtNbcvCRi6j8Ms/YT++Ojyj2UgmaWfYI+K7W2ZXvqxNdOGDkckqAAY5UHlt7/9rVxzzTVy0003yZtvvilHH320LFq0SJqbm92+NGAQK4SYpR9bqrA31EY+uTs5lNDa8uxGRSX42MJG6Gf0EgDAe0FlxYoVctlll8kll1wic+bMkZ///OdSWloq9957r9uXBkQRqqjoDh1r1HyfraIyeHty4ks/1tdL9PPSsfTDCH0AXuJqUOnp6ZE1a9bIwoULQxeUm2teX716ddTP6e7ultbW1rAb4MbpycpqqO0dopnWyQh9rdQEh61lupxh29EUOhyRoAJgFAeVPXv2SH9/v9TV1YW9XV9vbGyM+jnLly+Xqqqq4G3y5MkZulogVPmwnsitLcq9fQMxm2mdjNDXr+bW+PpQVYcR+gC8w/WlH6eWLVsmLS0twVtDQ4Pbl4RRuT3Z/3pw6Jtt6MlAUhUV3+CpsK4NfAtdA0s/ANyW7+Z/vLa2VvLy8qSpqSns7fp6fX191M8pKioyN8AN9oP7VHDoW19qzbT2rxoctpbhckbwP8dZPwA8xNWKSmFhocyfP19WrVoVfNvAwIB5fcGCBW5eGhB3C68qzBt8MGHkU3sigSNqj4qLA99Cu34IKgBGcUVF6dbkiy++WI4//ng58cQT5Sc/+Ym0t7ebXUCA19hPGFaF+bmD5qgMrqgk/nX1y+a61qMS5RoyPMsFADwXVL70pS/J7t275cYbbzQNtMccc4w8+eSTgxpsAW/t+vE/kRfl55n77t4UtyeHLbu4c3KxvVrE0g8Ar3A9qKirrrrK3ICRtvRTVGBVVPpDHyPJz1HJte24cWsqbPjAN4IKAHeNuF0/gDeWfnLCelTsFZXIJ/dECiPB3US2ZaXMT6YNLWtZDb22FS0AcAVBBXAgMnRYFZVu+xyVQWf9+JJrZHXprB+xVXUyHZYAIBJBBXDAPu/EXlHp6UuxohL4IB1069ayi72ZNhhU6KYF4DKCCuBA5Jk8wWbavv4oJyEnHjisMGPvUXG1mTZ4DRm9BAAYhKACOGBfoom19GMJNcXGX/6xlljM0k/gT2XGJ9OKLSy5dd4QAEQgqABJjpkPa6aNsvRjNaT6Py+xr6vhJj+QVDJfUYmy9EOPCgCXEVSAJA/ui9dMaz3ZJ7L8Y4US+wwTN5d+3Fp+AoBIBBUgqYFv4T0q9mZa66ndHlTiVSasPKBVmGBvS8bP+rEaet3rkwGASAQVwIFgr0nECH17M631hJ9vCyrxVlDsISE4mTbjc1QGb5EmqABwG0EFcKA/ov+kyDrrJ8rSj7VzJpGln+D25Jwcyc9zJySErsHZjiUASCeCCuBAZP9JqKIy+FDCAvvSz0Ciu37sU2Hdmkyru36oqADwBoIK4ECo6TXiUEJbUOkLfIwVYlS8uWn2ABSaYeLe0g+7fgB4BUEFSGbpJ/AnJ7T0E+pR6Y8SVOI94dsDkFXNyPSyS1hVxwoq/QQVAO4iqAAOWMs6wUMJoyz9WBUVnYdijVKJVx2xz14J7biRjLLPcglWdaioAHAZQQVwwAoPuUM001pNqdoUm2ivR9jJxYGQ0Jfhc3as/1xYVYceFQAuI6gADgQrH7nhQSVaRcVJZcIKCfYZJu4t/dCjAsA7CCpAklt4Yx1KaJ04rHNU8hPs9Qg76ydYhRGXln5CFSMOJQTgNoIK4IA9UMRa+rEqKmHbfONOpg2FBPcm04aWfkKzXEgqANxFUAEcGEhgjkq/vUclweFt9oFvbi27WNfoZlUHACIRVAAH7IEi1lk/1hN+Xm6ug2ZaGVSFsSozmRJe1Ql/GwC4haACJFN1yI19enJoe3Liw9uizTBx7VBCJtMC8BCCCuBAaN6J//XCQOmhu3fwwDdzbk+CO3iizjDJcEjoCzT8MkcFgJcQVIAkqw72ikpPf4yKSoLLOPbJtImGm+FmH1SX6G4lAEg3ggqQzMA3q5k2UFHp7fcFl2qse22kDU2ZTaxHxVQzXOpR6Qs8uII8KioAvIOgAjhgH3Wvigr8zbT2qoq9ohKsTCS86yfUyJrxpR/boDom0wLwCoIKkMzST8ShhKq7dyBs9oiTfpOwRlaXmmmtHpWCvFwm0wLwDIIK4IC9UVZpxcQ6eNCaThusTDjoN7FaQcJmmGQ4JPTaJuqG5qgQVAC4i6ACOGCfd2I1v0ae92M1oOrAt0T7TYJ9LQ6Wi4ZbaFBdbnAyLXNUALiNoAI4YA8UluAW5UBQsXpV9O2JLuNYn6PLLm5vTw7brcSuHwAuI6gAKWxPtjfUWtNprXsdr5/orh/7jptEzwcabr2BazCj/13aIg0AkQgqgAP2CbKW0NJPf1hlxR5U4i396PZmVWD7nEw301phyjTT0qMCwCMIKkCqSz/5sZZ+8hKuTFjVjIJcexUmLQ8h7jWET6bN7DUAQCSCCpBMM60tqEQeTBi29JNgM22oR8U+JC6zScW6Rg1Lwd1KVFQAuIygAiS1PTn0tshdP9F6VOI94QdnmOTnurY9OdhMa5tM25fhsAQAkQgqQDJbeK2Jb7aln1SaaXuj7BTK5NKPHopon6MSmkybuWsAgGgIKoAD9ifz2M20/vuiPOdBJWzHTQaXXbSZ1yrg6C6m4Bh/dv0AcBlBBUhqeSR3UFAJVlSs6oiDZRzrc+3j6zN5KGFnrz9cqZKCPCbTAvAMggrggH3eSWQzrdWj0tnjf9IvLsgLVl7ihQ7r/Wbpx1p2yWA1oysQVDQk6WNjMi0AryCoAEmeMBxr6acjEFRKC23bkx0s/bgxmdYKV1pN0WMBmEwLwCsIKkCSJwzHaqa1llHsQWWo0KHvswa+FefniVWsyWR/iHXNWgVSTKYF4BUEFWDYmmkHwioqJQkGFWvZxQoKeYEQZB1umMmgUlLo/28zmRaAVxBUgFSbaSPO+rGWUcoK8xNqprU3smroceOsny7b0o8KTaYlqABwF0EFSGqOSvTTk3UeSXtPX3Dpx/q4oSoToeZb/8nJVgbK5PbkYEUlEFSYTAvAKwgqgAP2ptdozbT+sCLBpZ9EGmOtJlyrP8SNybRWULGqQ6HJtAQVAFkaVLZu3SpLliyR6dOnS0lJicycOVNuuukm6enpCfu4d999V04//XQpLi6WyZMnyw9/+MN0XRIwrOfhRDuU0OpPUaWF+aHBaUNWVAbCqxl57u76UaHJtAQVAO7KT9cX3rBhgwwMDMjdd98ts2bNknXr1slll10m7e3t8uMf/9h8TGtrq5x77rmycOFC+fnPfy5r166VSy+9VKqrq+Xyyy9P16UBaaqoaFDpCxufb43aHzKoROy4cWPYWmuX/7orSwrMfbAJmB4VANkaVM477zxzs8yYMUM2btwod911VzCoPPDAA6bCcu+990phYaHMnTtX3n77bVmxYgVBBZ7UH22Oiq2Z1qpMaH9Koss4XTG2BmcyqLR0+CudVSX5EWEpY5cAAO73qLS0tEhNTU3w9dWrV8sZZ5xhQopl0aJFJtDs378/6tfo7u42lRj7DXB1jkrgZQ0c7bYdPyo44TWBikpJQfjW4EzOMGnp7DX3VYGKCpNpAYy6oLJlyxa544475Bvf+EbwbY2NjVJXVxf2cdbr+r5oli9fLlVVVcGb9rUArs5RCQSMHtvSjzbSquCE1wGPV1QCQaW6xP+PhtBkWkoqAEZYULn++uvNiO2hbtqfYrdjxw6zDPSFL3zB9KmkYtmyZaYyY90aGhpS+nrAcFVU9DDCtkCvR0Wxv6KSSDOtfeS+W0HlQERFJTSZNmOXAADD06Ny7bXXyuLFi4f8GO1HsezcuVM+8YlPyCmnnCK/+MUvwj6uvr5empqawt5mva7vi6aoqMjcAK+c9WMfod8aeMKvLLae8OM30x4MhpsC14atBZd+SgPXzWRaACM1qIwbN87cEqGVFA0p8+fPl5UrV0qubUunWrBggdxwww3S29srBQX+vyCffvppOfzww2XMmDFOLw1Iu76ou35CpycP2j2TQDNtW7f/c8qLAlUYFxpZWzrCKyrWH1V2/QDI2h4VDSlnnXWWTJkyxezy2b17t+k7sfeefOUrXzGNtDpvZf369fLb3/5Wbr/9drnmmmvSdVnA8MxRiXEoYVuXVVFJvJnWqqiUWUHFhamw+wK7fqqtigqTaQFk+/ZkrYxoA63eJk2aFPY+HTOutBn2qaeekqVLl5qqS21trdx4441sTYb356hEPZSwX1o7I5ZxEmimPdjdG9HXYn1OZkoqet0HAhWVuopi13YeAUBGg4r2scTrZVFHHXWUvPjii+m6DCA9hxLaljGL7D0qVkUlMI8kkXN7DkYu/WS4kXV3W3ewKdiqqGhTvHUN+g8L63UAyDTO+gEc0J099i3JkSP0YzbTDlGZONjdHxZUMj2ZtjkQVMZVFAUDib1ZmNUfAG4iqAAOaNXEviXZ3kzr71GJbKaV+Es/gSpMecTST8aCSqs/qIyvDO2ms5Z+MnkdABANQQVwwDrpOFpFRcPI/kBTqtVMm0hTqrX0UxGx6ydT/SHNbV3mfnxFKKjYN+jRpwLATQQVIEEaNnoDPSr2iooVVNSegz1hzbTW0s/QFZVAj0og3AS3Bme6ohJopDXXYKuoEFQAuImgAjjsT4kMJ1YzrdpzsDtim6/Erai0xWymda+iYu9RYekHgJsIKkCCtFk2si/F2qocuSmmtqwooWZa3VHTHnPgW2YCwq4Wf1Cpq4pVUcnIZQBAVAQVwGEjrSqwTabVnTLFEcElcntyrNCh1RTrXVYDbnCEfoYSwvb9neZ+8pjS4NvCdv2QVAC4iKACOG2kzc8dNFfEGtamasoKg++Pt9X4QLt/x09JQV7o9GTbDJN00xCywwoqNSXBt9tyCmP0AbiKoAI43Zps60mJFlTGlod6PazBcDGDSmf46PpMb0/WGSrae6P/zfrK0NKP/yR0/8tUVAC4iaACOOxRsTfPWqxlG1VbXhh8Od7Sz4GIwwAzfXpyw/4Ocz+xuljybTuZMl3ZAYBYCCqAw4qKvZHWYm1HVuNsFZXg0k+M0GHNXRlTWjg4IGQgIWwPBJVJ1aH+FDcCEwDEQlABHG5Pjrb0Yw14U5PGhHo9rNOTY1VUWgIj9+1LP8E5KhkICNv3De5PCV4HSz8APICgAiSouzf20o+9ojKpJlSdiNtMG1j6CetRCXyO5hTrpPF02bq3Y9COn8jrYI4KADcRVIAE9fT3x6yo2Kso02vLEm6mtZZ+qu1LPxkctvbhnoPmfsa48phLP0ymBeAmggowDBWVYyZXm/uywjw58pCq4NvjjcNvsSoqtmZae1AZavR+qrRa80GzP6jMHB8KV5HXQVAB4KbQwjqApHtUTpk5Vu666DiZOrYsOA8lbPkkTjNttO3J6Q4Je9t7pLWrz2xDnjZ2cFAJLVul7RIAIC6CCuCwomI/kNA+d+STR04Y9PZ4zbT7g9uTC6OOr0/n0o9VTdFlK3u4irwOelQAuImlH8DxZNrBT+qxxHuy39veHWX2in18vaTNh3vazf3MKP0p/usIXANLPwBcRFABEtTR4w8qpUWJB5Vgn0eMoLKnzb/0M85+crG9ouJLf0VlRm2MoBIc+EZQAeAeggrgNKgUOg8q0Zpi9dTkzl7/16y1D4nL0K6fDY1t5v6wuuhBxTqviKUfAG4iqAAJ6ujpM/elhYm3dg21c2bPwe7ggYRlRfkZ3XGjO37e29VqXp47MbRLKZPXAACJIKgADisqGiwSlT/EAYO72wL9KRWh/pRMDVtrbO2Sfe09JowcGqOiEgoqabkEAEgIQQVwGFTKipw300Zb+rEqKvazgRKdv5Kq93b6qymzxpVH3fFjriGwAsXSDwA3EVQAh0s/Jcks/USrqBzsGdSfkqlG1vWBoDJnYmXckMVZPwDcRFABnDbTFgxPM+2e4NJPtIpKTkYqKnMmVMa9dk5PBuAmggqQxqWfRJppo1ZU0tjIqo20bzXsNy/Ps437j1lRIacAcBFBBXDaTJvE0s9QzbTjbMPegp+XxvH1Ow50SlNrt7k264yiZGbAAEAmEFSABHUGtyc7qKjYqhJayYjceaPqq0InL2di6WfNx/5qytyJlVIyxGOhmRaAFxBUgAS1pzDwLdoT/q4Wf1CZUFWc0WbaN7b6g8r8qWOG/DgrLDFHBYCbCCpAgnSSrCqPGM6WcFCxPeH39A0Ee1SiBpU0VlTeCFRUjp9aM+THMUIfgBcQVIAE9PYPBHtUqkoKUq6oNLV2iT7/F+bnSk1ZYew5KsMcEnTI24ZG/46fE6YlVlFJR58MACSKoAIkoLWzN/iyk4qKtXMmMqjYl32sM3UyMZn2xc27TUCaXV8h4ysHV3Ki9qhQUQHgIoIKkIDWrtCyT35eruMR+oODSmfMZZ90Lv28uHmPuT/zsHFxP9a6hsgmYADIJIIK4KCi4mTZZ6iln1BFpSRjW4M1cGhFRZ1+6LiEq0Hs+gHgJoIKkIDWLn9QqShOfNlH6bJOtCWUXQeGrqgEQ8IwVjPW7mgx81P0UMXj4/SnpLuhFwASRVABEtDa6V/6qXRYUYn1hL9ziK3JsT4nVY+v3WXuzz5ifMyDCKNPpiWoAHAPQQVwUFGpLHYeVKItoTQmuvQzTCFBl33+Eggqn5o3IaHPYYQ+AC8gqAAJ2N/hP+m4utR5ULEaau1BZfv+DnM/sbokTriRYbFuR6s07OuU4oJc+cTs+P0pyuoZZukHgJsIKkAC9rT1xDxAMJ7IcfgHu/tkf4e/QjO5ZuiKynCFhN+vaTD358yuk9IEzypK58GIAJAoggqQgL3t1knHg4ezOX3Cb9jXEazOVMRYShrOqbCdPf3yyFs7zMtfOmFywp/Hrh8AXkBQARJgjbsfm0RQsZZ++gbCg8rkMaUxPyc4mXYYQoL2pugcmEljSuS0WbUJfx49KgC8gKACpHvpJ6Iy0bC/c8hln+FcdtEm2vtf+di8/OUTJgeXoRKRjlkuAOAUQQVwsPQztqwo6YrKQKAxNpGKSn6gpNLbn1pIePWjffJ2wwFzptAXHSz7pGuWCwA4RVAB4tBKiB7mp2ornC/9WFWMvkBSsXb8TKqJHVQK8gKfk+K2n7ue+8Dcf/H4STK+YuizfSKx6weAFxBUgASqKfpcrQWGmtLhaKYNLP2Mib30UxBICXpqc7LeaTggz2/abSbjXn76TMefb1VUOOsHQNYHle7ubjnmmGPMOPG333477H3vvvuunH766VJcXCyTJ0+WH/7wh5m4JCBhOwI9JXUVxY4OJIwMKn39PvOk3xCoqEwesqKS2tKP/nf+4y/vm5c/d+whMmVs7P9W/G3VSV0CAIycoPLP//zPMnHixEFvb21tlXPPPVemTp0qa9askR/96Edy8803yy9+8YtMXBaQkB2Bc3kOGaICMpQ8W6+HLiF19PSb6swhMYa9qfzA0k+yFZW/vd8sr320T4ryc+W75x6e8nUDgFucnbCWhCeeeEKeeuop+cMf/mBetnvggQekp6dH7r33XiksLJS5c+eaisuKFSvk8ssvT/elAY4qKkMFi8R2z4R2/Gh1ZqjzdgoDFRVrS7PTuSn/GaimXHra9JjTbxO/boIKgCytqDQ1Nclll10m999/v5SWDi49r169Ws444wwTUiyLFi2SjRs3yv79+2MuI2klxn4DMlFR0TkkyQgu/QwMhHb8DLE12V5R6elzXlH5ryc3yEd72qWuskiuPMt5b4olUFBhMi2A7Awquka+ePFiueKKK+T444+P+jGNjY1SV1cX9jbrdX1fNMuXL5eqqqrgTftagIxUVFIMKgP2/pQhtibbe1SsnUKJ+vuWPXLfy1vNyz/8P0cndYiihaUfACMyqFx//fWmKXao24YNG+SOO+6QtrY2WbZs2bBesH69lpaW4K2hwX+GCZD2HpUUl1C03cTa8TPU1uRkm2lbOnvluw+9Y17+vydPkTMPS+zwwVhY+gEwIntUrr32WlMpGcqMGTPkmWeeMUs7RUXhA7K0unLRRRfJr371K6mvrzfLQ3bW6/q+aPTrRX5NIF10hoguo6jptWWpVSYGBoIzVIbammyfo+KkmfaWx9bLrpYumTq2VL7/qSMkVdauH3IKgBEVVMaNG2du8fz0pz+Vf//3fw++vnPnTtN/8tvf/lZOOukk87YFCxbIDTfcIL29vVJQ4C9RP/3003L44YfLmDFjnF4aMOw0WHT3DZjJrpPiLNfEH/jms/WoDP21QpNpEwsqT65rlIff3GFmpqz44tEJn5A85HUHelQY+AYgK3f9TJkyJez18vJycz9z5kyZNGmSefkrX/mK3HLLLbJkyRK57rrrZN26dXL77bfLbbfdlq7LAhzZ3HTQ3M8cVx5cCnFKtwgr3ZZsLSPFCyoajKzZK4kcmHjDI2vNy984c6bMn1ojw2E4T3AGAM9uTx6KNsPq1uWlS5fK/Pnzpba2Vm688Ua2JsMzNjf7g8qh4/1BOxnlRf4/Zh/ubjc9J7qsU19ZnND5QD1xKiratL7s4bWyt71HZtdXyNULD5XhEhr4RlABMAqCyrRp06KO4j7qqKPkxRdfzNRlAI5samoz97NSCCplgaCysdG/lV7nmsSrzgR3/cSpqOhyz9PvNZnws+KLx0hRfuzZLE5RUQHgBZz1Awzhne0HzP2Rh1Ql/TXKCv3hYf1Of1CZEmfZJ9Fm2ubWLtNAq65eeJjMmVgpw8k6LiCR5ScASBeCChBDa1evWa5RR02qSrmi0tzWbe5nJLB7KN6hhFqd/Nc/rpPWrj4Tor5xxoykry/2NaQ2xh8AhgNBBYhh7faW4BTZseVFKQcVSyLbnPPjzFF5buNu+ev6JtPL8l8XHpXUYYkJhyV6VAC4iKACxPB2g3/Z5+hJ1Sl9HWvpxzJ9XHlK1Qxtbl3+ROgsn+Fe8hl0MGISY/wBYLgQVIAYXv5gj7k/fmpqM30iKyqJLP1YBxbqDJdIj7y1QzY1HZSqkgJZetYsSZeC3OQPRgSA4UJQAaLo6u2X17f6D8Y87dDURtGPqygKm6mSyGnGJYGgorNXIntT/ueFD83L3zhzhlSVJn+WTzwF+fSoAHAfQQWI4vWt+8zJxTrvZOa45EbnW+xnBM2dWJnQ4LiSwHKRBia7lz/YKxub2qS0ME8uOnGqpJPT6bgAkA4EFSCKFzbtNvenHVprDtpMxZSxpcHptPMTXEayKiqdERWV+1d/bO7/z/xJaa2mJHswIgBk1WRawIt0eeUvaxvNywuPGJ/y19MhbD/+wtFmp843E+wpsXpUOm0VFd0u/czGZvPyl08IP6IiHayG3j4qKgBcRFABIry7vcWcyaPLK2cdnnpQUZ85eqK5Jcpa+rEHlafWN5nlKJ2Se8SECkk3q6LSQ0UFgItY+gEiPPbOTnN/9uzxwcpGpllLPxpMrLN2nli7y9x/5qiJKS9HOdmeTEUFgJsIKoCNNq/+4c3t5uXPH3uIa9dhBRXrmrShdfWHe83rC+cMT5UnnkJrhD7bkwG4iKUfwObJdY2yv6NXJlYVD9uyTzKs5ltr+UfPCdKtymPLCuWI+vQMeItkTbvVqg4AuIWKCmBror3v5a3m5S+fOCWhbcTpkpubI8UF/j+eHd398tIW//C5U2bVmvdlQrCZdoCgAsA9BBXANqNEx+ZrNePLJ052+3LM5FnV0tkrL232b5c+fVZtxv77bE8G4AUEFSDgjmc2m/t/OnGKjK8odvtyZExpobnftq9D3gkckHjqoW4EFSoqANxDUAFE5NkNzfLKh/vMcsflZ8wQL6gp8weVv6zbZXb+6KnL9im36aYnMyuCCgA3EVQw6mmz6L89/p55+dJTpyd0Fk8mjAkElcff9W9LPi2Dyz72ikofSz8AXERQwah338sfyYe726W2vFCuOjt9pxE7VRNY+rGcnsFln/BmWp9pNAYANxBUMKptaT4oP35qk3n5nxfNlori9J6f48QM22GIGhp0x08mWduTFQ21ANxCUMGopRNXr33oHbP0c8Zh4+QLx08SL5kzITQvZcHMWikvyuzYI2vgm6JPBYBbCCoYtW772yZ5p+GAVBTny39deGRGxtI7ceyUMSasFObnyrfPyfySlDVCX9GnAsAtTKbFqPTX9Y1y57MfmJf/4/NHyoQqbzTQ2mlAeWTpKSYklGW4mmLf9aN6qKgAcAlBBaOyL+Xa371jXr7k1GnyWQenGmdaUX6euJBRDK0waW+M9qew9APALSz9YFRpbuuSS+57TQ5298mJ02rk+586wu1L8rTifP/hiN2c9wPAJQQVjBptXb1yycrXpWFfp0ypKZU7LzouOCsE0RUX+oNKZ0+/25cCYJTib2mMCrqz58pfv2lOIdYTiP/fpSfKuIoity/L80qtoNJLUAHgDoIKRkVIueo3b5oTiPWJd+UlJ8i02tCMEsRWUkBFBYC7CCoYFSHlqfeazC6au786X46aVO32ZY0YxVZQoaICwCUEFYyakPI/XzteTj90nNuXNaKw9APAbQQVZG1I+db/hoeUMw8jpCS/9NPn9qUAGKUIKsjakPLX9YSUVJWw6weAywgqyOqQ8ouvziekDENFpYOlHwAuYTItsiqkLP3Nm/L0e6GQctbh492+rKzoUenojh5Udrd1y9+37JFt+zpEB+7XlBfKjNpymXdIpadOogYwchFUkBW6+/rlm79+U1ZtaCakDKOq0kJzf6CzJ+zt7d19svyJ9+V/X2uQ/oHBBxbqMUFHHlIlJ88cK+fMrpP5U8dInu3sIABIFEEFI15Xb79c+es18uzG3VKUnyu/vJjdPcOlusRfFTnQ0Rt8mx4/8NV7XpW3th0wr2sgOWJCheTm5Ehja5dsbjooOw50yjvbW8zt7uc/lJqyQjl79nj5hzl1csah44K9LwAQD0EFIz6kfOP+NfL8pt1SXJAr91x8gpw6q9bty8oaY8oGB5XvP7zWhJSqkgK566Lj5JQo3++dBzrl1Y/2ygub9sgzG5plX3uP/H7NdnPTn5MGyXPn1Mk5R9SZEAMAsRBUMGLpTpTL739DXty8xzR93rP4eDllJiFlOFWX+EPE/g7/0o/2//zpnZ1mGefexSeYJZ1oJlaXyOePnWRuevLy61v3mc/V2/b9ncGXdTXoxOk1cu6celNtmVxTmtHHB8D7CCoYsSFlya9el5c/2Osfi7/4BDlpxli3LyvrVJf6Kyr723tkYMAnP/7rRvP610+fHjOkRNKDHzVA6u3GT8+R93e1yVPvNcpT65vkvV2t8sqH+8ztB39+T+ZMqJRz59aZ4KLLSTk59LUAox1BBSNOR0+fXHrf6+bJrawwT+679EQ5YVqN25eVlQ6pLjH32nuilZSNTW1SUZwv3zxzVlJfT4PHnImV5nb1wsOkYV+HGcr31PpGU3XR4KK3n/xts0waU2L6WnSZ6OQZNewiAkapHJ/PN7hlfwRpbW2VqqoqaWlpkcrKSrcvB2mmu00uWfm6vLZ1n5QX5cuvLtXlB0JKuuhfD/Nu+qu09/RLfm6O9A345NvnHCrf+YfDhv2/pX0sq95vMsHlhU27pbtvIPg+/W8fO6VaTps1Tk47tFaOnlQl+XmMgQJGw/M3QQUjhu42WXzva/LGx/ulQkPKkhPluCmJLT8geZ+54yVZu6PFvKzVlJeuO9s00qa7avbS5j2m/+jFzbtl696OsPdrJe24qWPM8pNW046ZXC1lRRSIgWx8/uZPNkaEtq5eWbzydVmjIaU4X+5fcpJ5ckL6nXPE+GBQ+cYZM9IeUlRpYb6cO7fe3JQuEb20xR9a/r5lr7R09gZCzB7zfm3u1f4WrbrMm1hllpYOq6swM3UAjGxUVDAiQsrF974mb247IJXF+fLrr58kR00ipGSyknXjH9dJbXmR/POiw11fctEBc5ua2kxl7Y2t++SNrfvN3JZIhXm5clh9uQkuh9dXyMxx5TJrfLlMqCqmSRfwAJZ+kBVaAyHFmtvxwNdPknmHVLl9WfCYXS2dJrCs29Ei63a2yLodrabqEo3uEtPQMnNcmUyrLZNJY0pN07A272qIcTuIAaNFqxeCyuOPPy4/+MEP5N1335Xi4mI588wz5dFHHw2+f9u2bXLllVfKs88+K+Xl5XLxxRfL8uXLJT8/8RUpgkp2h5Sv3fOavN1ASIEz+teazmtZHwgtm5vb5IPd7bJ1T7tpCI5F57pMqCqRidXFMq6iSMaVF5lKUm1F4L680Nzr8qM2c2dDqOnrHzCNy/5bvzkzy7zcG3p9QnWJTK8tc/tSkWVc71H5wx/+IJdddpn853/+p5x99tnS19cn69atC76/v79fzj//fKmvr5eXX35Zdu3aJV/72tekoKDAfA5GN/3X8NfufU3eaThgZnn8egkhBYnTpR0dHqe38+ZNCL5dh8/pAYpbmg+am/a+6LKRhhq91ydlvY+2lBSNNvXqtmkNLpUl/nsdPlhckGeOczA3fT1wr6/rXBkNRHqNeuyAvqz3uhplXg9kn/4BvQ3Y7n0mZOl9v88n/f3+e31M/lDhDxb2kGEPILE+ZojcFubQ8eVy3rx6WTS3XuZOrGT5DBmTloqKhpJp06bJLbfcIkuWLIn6MU888YR8+tOflp07d0pdXZ15289//nO57rrrZPfu3VJYmNhYbSoqWRpS7nnVnBOjIUUrKXMnElKQXjrQbs/Bbtl+oNMcAaAnQ+vre9p6/Pfm1iN727ulqze0dTqbFOTlSFG+P1BpI7IVrD6KqETpUpkGlkVz6+T4aTUcOImRV1F58803ZceOHZKbmyvHHnusNDY2yjHHHCM/+tGPZN68eeZjVq9eLUceeWQwpKhFixaZpaD169ebz4umu7vb3OwPFNkZUsaYkHKy2cEBpFtubo6Mryw2t3jb3rXyok3ebV195qbLlPp6a1efOX/KqlZooLEqGPp2fV0rIgM+n+g/EfVen//134vWy3qv9Mk/z1RYcswcGete32beF3ibNg0XFWioCK/iBF/Wt5v3hz6mOOzj/e/Xr6NfL9afy2c3NMuT6xrNuVpacbr37x+Z29iyQnP8gQaXU2aNNV8PGE5pCSoffvihub/55ptlxYoVprry3//933LWWWfJpk2bpKamxoQXe0hR1uv6vli0h0UrNcg++i9W7UnRyaR6UJ1WUo6YQEiB92i1YWx5kbmNBtoj9rljDzE3Pb5Ct4k/ub5RVr3fLHvbe+TB1xvMTft2PjF7vKm0nHX4ePM6kCpHnWDXX3+9WZcc6rZhwwYZGPCXRW+44Qa58MILZf78+bJy5Urz/oceeiilC162bJkpE1m3hoaGlL4evEFL7V+8e7UJKdqw+JvLCCmAF5UU5pn5Niu+eIy88S8LTf/YV0+eKnWVRWYr+2Pv7JSrfvOWHPdvT8uS+16X373RYKYOA8lyFHevvfZaWbx48ZAfM2PGDNMYq+bMmRN8e1FRkXmf7vRR2kT72muvhX1uU1NT8H2x6NfRG7KH7sS46JevmnLyxKpiMydlxrhyty8LQBzav6JHGujtls/Olbe3H5C/rm+Uv65rNNOEV21oNjfrlGxdHlp4BKdkI41BZdy4ceYWj1ZQNExs3LhRTjvtNPO23t5e2bp1q0ydOtW8vmDBAvmP//gPaW5ulvHjx5u3Pf3006ahxh5wkN02NLbKV+95zTQu6vZHDSnWQXgARg7tb9HeHr1df95s2dR00IQW7Wuxn5J9y2PvmcF758web5aJ9BgEDTxAxueoXH311fL73/9e7r33XhNOtJH2scceM0tDY8aMMduTtcF24sSJ8sMf/tD0pXz1q1+Vr3/96462J7PrZ+R69cO9cvn9a0yj3uz6CjMWX2dXAMguug1cQ4seOKnHYGhDsUW3dJ9x2DgTXM48bNyo6fuBuD/wTSso2k9y//33S2dnp5x00knyk5/8RObOnRv8mI8//tjs8nnuueekrKzMDHy79dZbGfg2CvzpnZ3y3d+9Iz39A3LclGpZufhEqSpN/xkyANzV0tErL2zebXYRPbdpd1j/io5mOeqQKjl1Vq25abVFZ9IgO7keVDKFoDKy6K/b3S98KLc+scG8rrsDbv/ysfxlBIxCWlnRydMaWp7Z0GyWiCJ3Vx0/dUwwuBx5SBUzW7IIQQWeoxM0b3lsvfz6FX9D9SWnTpN/OX8Of/EAMBpbuuTvW/b4bx/skabW0MwsVVGUL8dOHWPCi96OmVJtTtrGyERQgedmpCx94E159aN9pryrAWXJadPdviwAHqVPTXo+08sf7JGXNu+R1R/uNcP17PQfOUdMqJDjp9bIcVPHmGWjqWNLGe8/QhBU4Blrt7fIN+5/Q3a2dJkBUCu+eLSZwwAATpaJ3t/Vappx3/h4v6zZus/8nRJJm3PnTaySeYdUmvPBdLlo2tiymFN34R6CClynv1q/X7Ndbnh0nRk5PqO2TH7xtfkya3yF25cGIEsGRVqh5a2GA7Khsc38XRPt8MhZdRXmYEW9HVZXYbZI6ygEAox7CCpwlZ578q+PrpNH395pXteth7d9+RipLGZnD4D09cFtbjoo63a0yLqdLbJ2R4upwsQ6RLK0ME9mjiuXKWNLZWpNqUyxbmNLZUJVCf1zaUZQSZGeZ/H42l3y+WMP4ZfVoXcaDsj/9+Bb8vHeDvO9+87CQ+WbZ83iXy4AMq6vf0C27m03AWZz80HZ1NQmW5oPyoe72814hKFOktaKiwYWPR6grqpY6iqKpV7vK/VWZOY+ZfMhjD6fz4Q8PXRTd2bqmU/DiaCSotue3iS3r9ps1jdv/uwcmT+1Zti+djb/a+bu5z+Qn/xtszkSXv+Q//SfjuF7B8CTAebjfR3yQfNB2bavI+y2fV/nkCEmsiozprRQqksLwu719PeK4gIpLcqTssJ883FlReH3uv26IDdX8vNyzHRevTn9h7HP55Pefp+5Xl32Ct76/ad2dwde139861lM7Xrr6fffm5f13v++jp4+OdjdL22dvWYQpwYU/drqX84/Qr5++gxx4/mbfV0x1FYUmaYsLR1eeNdqueDYQ+T6T842R8BjsPU7W+R7D70bnIPwyXn1cusFRzHEDYAn5eflmmUfvUVr3G1q7TKhRe/11tjSHXq5tUuaW7tNOOjo6ZeOnk5zVtlw0A1LGl60oqNVaJ/PH0Y0LujLA4GX9f/0Jf1HYbrLDZqdNPC4hYpKnC21P3pyo/xuTYP5RdCGrEtOnW621Y4pKxzW/9ZI1dXbLz97dov87LkPzB8Y/dfEzZ+ZK/94zES2CALIWvrU2drVJwc6emR/R6/sb9d7/8v+t/WYSoVWLTTMaOWioztwH6hoaNAZ7mfgvNwcKczLlaKCXHOvVRu9lQWqOrrzUis6ZYFKj76sb9PKj3lfYb5UlhRIZUm+6SnUl/W5Lx1/n7P0M8w9Fzf9ab2ZoKj0h/Z/F0w1R5tPGjM6TwHVX5un32uSf3v8PWnY1xmsovzgH+dxXg8AJEirN7psrre+fp/0DgTu+wfM+3JzckyVJUf/F8gK5nV9u1aGcnNMn4wVSEZSTyVBZZgNDPjMgVo/XbU5uLyhvw+fOHy8fOWkKXL6oePML8lo8MHug+YE1Bc27Tav11cWy42fmSOfOnKC25cGABghCCppot+uv73fLCv//pG8/MHe4Nsri/PNELNPHVkvC2bUSklhXlaOt77jmc3y29cbzDKPlhW/fvp0WfqJWaZ8CABAoggqGaos/ObVbeYk4N1toTMp9AlcT/087dBaOWXmWJkzsXJEb2HT003vem6L/L/VHwcbqs6ePV5u/PQcmVZb5vblAQBGIIJKBuk64htb98lf1u4yfRuRY521e3vOhEo5enK12e58aF2FzBxXZraueX1o2y9f/Ejueekjs3VN6UFg3110uJw8Y6zblwcAGMEIKi7Rb+eHe9rN6Z8vbt5jzqXQikQ0OjBIxzjrJEQdKjSxukQmVhXLhOoSGVteaE4KTUentfbbHOjslX3t3bL3YI+5vj3tPbLrQKfsauky2+x2tXSapR5rD/3ciZUmoJx12Dh28wAAUkZQ8Qj99m7f32l2DOnuIW3E1amIzbaloli0e1snAVaXFJh5JLp1rCjQ2a1LSfqyDgjSvfQDgb32AwP+ffbaQ6LDe/x7/G1b5Lr7zLY5/fhEzBhXJtf+w+FmRw+TZQEAw4WBbx6h1YfJNaXm9pmjJwbfrlP/NLBon8uO/f4KhlXN0EqGhgpdUtJqR6yKTKq0AXhseZHUlBWam1XNmVBV7B8dHajwUEEBALiFoOISrZRow63eYg1Sa+3sNUs0B3SYUEePGYHc3RcYi9yro5H7zWhkDRK6114LHlr10FyRl5MjpWZ4T54Z8lNqBvv478eWFZqBdVqNAQDAywgqHqUHQOmNkf0AgNGMf1IDAADPIqgAAADPIqgAAADPIqgAAADPIqgAAADPIqgAAADPIqgAAADPIqgAAADPIqgAAADPIqgAAADPIqgAAADPIqgAAADPIqgAAADPGvGnJ/t8PnPf2trq9qUAAIAEWc/b1vN41gaVtrY2cz958mS3LwUAACTxPF5VVRXz/Tm+eFHG4wYGBmTnzp1SUVEhOTk5w572NAA1NDRIZWWljDY8/tH9+NVo/x6M9sevRvv3YLQ//nR+DzR+aEiZOHGi5ObmZm9FRR/cpEmT0vrf0B/MaP0FVTz+0f341Wj/Hoz2x69G+/dgtD/+dH0PhqqkWGimBQAAnkVQAQAAnkVQGUJRUZHcdNNN5n404vGP7sevRvv3YLQ/fjXavwej/fF74Xsw4ptpAQBA9qKiAgAAPIugAgAAPIugAgAAPIugAgAAPIugEsOdd94p06ZNk+LiYjnppJPktddek2x08803m4m+9tvs2bOD7+/q6pKlS5fK2LFjpby8XC688EJpamqSkeyFF16Qz3zmM2Yaoj7eRx99NOz92l9+4403yoQJE6SkpEQWLlwomzdvDvuYffv2yUUXXWSGH1VXV8uSJUvk4MGDkg2Pf/HixYN+J84777ysefzLly+XE044wUyzHj9+vHzuc5+TjRs3hn1MIr/327Ztk/PPP19KS0vN1/ne974nfX19ki3fg7POOmvQ78EVV1yRFd+Du+66S4466qjgALMFCxbIE088MWp+/nfFefye+9nrrh+Ee/DBB32FhYW+e++917d+/XrfZZdd5quurvY1NTX5ss1NN93kmzt3rm/Xrl3B2+7du4Pvv+KKK3yTJ0/2rVq1yvfGG2/4Tj75ZN8pp5ziG8n+8pe/+G644Qbfww8/rDvefI888kjY+2+99VZfVVWV79FHH/W98847vs9+9rO+6dOn+zo7O4Mfc9555/mOPvpo3yuvvOJ78cUXfbNmzfL90z/9ky8bHv/FF19sHp/9d2Lfvn1hHzOSH/+iRYt8K1eu9K1bt8739ttv+z71qU/5pkyZ4jt48GDCv/d9fX2+efPm+RYuXOh76623zPe0trbWt2zZMl+2fA/OPPNM83ef/fegpaUlK74Hf/rTn3yPP/64b9OmTb6NGzf6vv/97/sKCgrM92M0/Pz/FOfxe+1nT1CJ4sQTT/QtXbo0+Hp/f79v4sSJvuXLl/uyMajoE040Bw4cML+8Dz30UPBt77//vnlyW716tS8bRD5RDwwM+Orr630/+tGPwr4PRUVFvv/93/81r7/33nvm815//fXgxzzxxBO+nJwc344dO3wjSayg8o//+I8xPyebHr9qbm42j+f5559P+Pde/2LOzc31NTY2Bj/mrrvu8lVWVvq6u7t9I/17YD1Zffvb3475Odn2PRgzZozvl7/85aj8+dsfvxd/9iz9ROjp6ZE1a9aYcr/9PCF9ffXq1ZKNdFlDlwFmzJhhyvla0lP6fejt7Q37Xuiy0JQpU7L2e/HRRx9JY2Nj2GPWsyh0+c96zHqvyx3HH3988GP04/X35NVXX5Vs8Nxzz5ly7uGHHy5XXnml7N27N/i+bHv8LS0t5r6mpibh33u9P/LII6Wuri74MYsWLTKHt61fv15G+vfA8sADD0htba3MmzdPli1bJh0dHcH3Zcv3oL+/Xx588EFpb283SyCj7effH/H4vfizH/GHEg63PXv2mB+c/Qeg9PUNGzZIttEn4Pvuu888Ie3atUtuueUWOf3002XdunXmCbuwsNA8KUV+L/R92ch6XNF+/tb79F6fxO3y8/PNX/LZ8H3RfpQLLrhApk+fLh988IF8//vfl09+8pPmL6e8vLysevx6+vrVV18tp556qvkLWSXye6/30X5HrPeN9O+B+spXviJTp041/4h599135brrrjN9LA8//HBWfA/Wrl1rnpi1H0X7UB555BGZM2eOvP3226Pi5782xuP34s+eoDLK6ROQRZurNLjoL+jvfvc700iK0efLX/5y8GX9V5P+XsycOdNUWc455xzJJtowqaH8pZdektEq1vfg8ssvD/s90OZy/flreNXfh5FO/3GmoUSrSb///e/l4osvlueff15Gi8NjPH4NK1772bP0E0FLXfqvxsgOb329vr5esp3+K+Kwww6TLVu2mMerS2EHDhwYNd8L63EN9fPX++bm5rD3a7e77oTJxu+LLgnqnwv9ncimx3/VVVfJn//8Z3n22Wdl0qRJwbcn8nuv99F+R6z3jfTvQTT6jxhl/z0Yyd8DrZrMmjVL5s+fb3ZBHX300XL77bePmp9/YYzH78WfPUElyg9Pf3CrVq0KK43q6/b1u2ylW0w1NWuC1u9DQUFB2PdCy3/aw5Kt3wtd7tA/aPbHrOuu2nthPWa917/EdC3b8swzz5jfE+sPdDbZvn276VHR34lsePzaQ6xP0Frq1uvWn7ldIr/3eq+lc3tge/rpp81WT6t8PpK/B9Hov76V/fdgJH8PIunvb3d396j4+Q/1+D35sx/29tws2Z6suzzuu+8+s8Ph8ssvN9uT7R3O2eLaa6/1Pffcc76PPvrI9/e//91sN9NtZroLwNqmp9sWn3nmGbNNb8GCBeY2krW1tZktdXrTPwIrVqwwL3/88cfB7cn68/7jH//oe/fdd80OmGjbk4899ljfq6++6nvppZd8hx566IjZnjvU49f3ffe73zW7G/R34m9/+5vvuOOOM4+vq6srKx7/lVdeabaf6++9fftlR0dH8GPi/d5b2zPPPfdcs733ySef9I0bN27EbE+N9z3YsmWL7wc/+IF57Pp7oH8WZsyY4TvjjDOy4ntw/fXXmx1O+tj0z7i+rrvWnnrqqVHx879+iMfvxZ89QSWGO+64w/yi6jwV3a6s8yKy0Ze+9CXfhAkTzOM85JBDzOv6i2rRJ+dvfvObZutaaWmp7/Of/7z5C20ke/bZZ80TdORNt+VaW5T/9V//1VdXV2cC6znnnGNmDdjt3bvXPDGXl5ebLXmXXHKJeZIf6Y9fn6j0Lx/9S0e3aE6dOtXMU4gM6SP58Ud77HrTuSJOfu+3bt3q++QnP+krKSkx4V5Df29vry8bvgfbtm0zT0w1NTXmz4DOyfne974XNktjJH8PLr30UvO7rX/v6e+6/hm3Qspo+PlfOsTj9+LPPkf/b/jrNAAAAKmjRwUAAHgWQQUAAHgWQQUAAHgWQQUAAHgWQQUAAHgWQQUAAHgWQQUAAHgWQQUAAHgWQQUAAHgWQQUAAHgWQQUAAHgWQQUAAIhX/f9acsGY2X7eJAAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def try_neuron_simulation():\n", " brainstate.environ.set(dt=0.1 * u.ms)\n", "\n", " I = braintools.input.section(values=[0, 0.05, 0], durations=[50 * u.ms, 200 * u.ms, 100 * u.ms]) * u.uA\n", " times = u.math.arange(I.shape[0]) * brainstate.environ.get_dt()\n", "\n", " neu = HTC(1)\n", " neu.init_state()\n", "\n", " t0 = time.time()\n", " vs = brainstate.transform.for_loop(neu.step_run, times, I)\n", " t1 = time.time()\n", " print(f\"Elapsed time: {t1 - t0:.4f} s\")\n", "\n", " plt.plot(times.to_decimal(u.ms), u.math.squeeze(vs.to_decimal(u.mV)))\n", " plt.show()\n", "\n", "if __name__ == '__main__':\n", " try_neuron_simulation()" ] }, { "cell_type": "markdown", "id": "d3031c1f0d2730b3", "metadata": {}, "source": [ "## Thalamic Network Modeling\n", "\n", "After modeling the individual neurons, we proceed to model the complete thalamic network.\n", "\n", "Most common synaptic models in neuroscience are readily available in `braincell`.\n", "However, in this example, we encounter a relatively uncommon NMDA synapse with magnesium block, which needs to be modeled manually.\n", "\n", "The `MgBlock` class we define inherits from `SynOut`, indicating that it is a synaptic output model.\n", "\n", "According to its equation:\n", "\n", "$$I_{\\text{NMDA}}(V) = \\dfrac{g \\, (E - V)}{1 + \\exp \\!\\left( -\\dfrac{V + 25}{12.5} \\right)}$$\n", "\n", "it can be easily implemented in the model." ] }, { "cell_type": "code", "execution_count": 179, "id": "c1d2d018c1a967a6", "metadata": {}, "outputs": [], "source": [ "class MgBlock(brainpy.state.SynOut):\n", " def __init__(self, E=0. * u.mV):\n", " super(MgBlock, self).__init__()\n", " self.E = E\n", "\n", " def update(self, conductance, potential):\n", " return conductance * (self.E - potential) / (1 + u.math.exp(-(potential / u.mV + 25) / 12.5))" ] }, { "cell_type": "markdown", "id": "8cf931923fb72510", "metadata": {}, "source": [ "In addition, since the gap junction connections in this example are based on spatial distance and probability, we define a `ProbDist` class.\n", "\n", "The purpose of the `ProbDist` class is to generate a set of synaptic connections given the sizes of the presynaptic and postsynaptic neuron populations, following these rules:\n", "\n", "- Only connect neurons whose spatial distance in the 2D plane does not exceed `dist`.\n", "- Each neuron pair has a probability `prob` of being connected.\n", "- Each presynaptic neuron has a probability `pre_ratio` of participating in the connection.\n", "- The `include_self` parameter controls whether self-connections are allowed.\n", "\n", "Based on these rules, we can model gap junctions in a detailed and precise manner during network construction." ] }, { "cell_type": "code", "execution_count": 180, "id": "a39e1cb5dab0f04a", "metadata": { "ExecuteTime": { "end_time": "2025-09-22T03:31:52.209278Z", "start_time": "2025-09-22T03:31:52.204373Z" } }, "outputs": [], "source": [ "class ProbDist:\n", " def __init__(self, dist=2., prob=0.3, pre_ratio=1.0, include_self=False):\n", " self.dist = dist\n", "\n", " @numba.njit\n", " def _pos2ind(pos, size):\n", " idx = 0\n", " for i, p in enumerate(pos):\n", " idx += p * np.prod(size[i + 1:])\n", " return idx\n", "\n", " @numba.njit\n", " def _connect_2d(pre_pos, pre_size, post_size):\n", " all_post_ids = []\n", " all_pre_ids = []\n", " if np.random.random() < pre_ratio:\n", " normalized_pos = np.zeros(2)\n", " for i in range(2):\n", " pre_len = pre_size[i]\n", " post_len = post_size[i]\n", " normalized_pos[i] = pre_pos[i] * post_len / pre_len\n", " for i in range(post_size[0]):\n", " for j in range(post_size[1]):\n", " post_pos = np.asarray((i, j))\n", " d = np.sqrt(np.sum(np.square(pre_pos - post_pos)))\n", " if d <= dist:\n", " if d == 0. and not include_self:\n", " continue\n", " if np.random.random() <= prob:\n", " all_post_ids.append(_pos2ind(post_pos, post_size))\n", " all_pre_ids.append(_pos2ind(pre_pos, pre_size))\n", " return all_pre_ids, all_post_ids # Return filled part of the arrays\n", "\n", " self.connect = _connect_2d\n", "\n", " def __call__(self, pre_size, post_size):\n", " pre_size = np.asarray([pre_size[0] ** 0.5, pre_size[0] ** 0.5], dtype=int)\n", " post_size = np.asarray([post_size[0] ** 0.5, post_size[0] ** 0.5], dtype=int)\n", " connected_pres = []\n", " connected_posts = []\n", " pre_ids = np.meshgrid(*(np.arange(p) for p in pre_size), indexing='ij')\n", " pre_ids = tuple(\n", " [\n", " (np.moveaxis(p, 0, 1).flatten())\n", " if p.ndim > 1 else p.flatten()\n", " for p in pre_ids\n", " ]\n", " )\n", " size = np.prod(pre_size)\n", "\n", " for i in range(size):\n", " pre_pos = np.asarray([p[i] for p in pre_ids])\n", " pres, posts = self.connect(pre_pos, pre_size, post_size)\n", " connected_pres.extend(pres)\n", " connected_posts.extend(posts)\n", " return np.asarray(connected_pres), np.asarray(connected_posts)" ] }, { "cell_type": "markdown", "id": "b807bf03a5a834d9", "metadata": {}, "source": [ "In this thalamic network, neurons are arranged in proportionally distributed grids:\n", "\n", "- 49 HTC neurons arranged in a $7 \\times 7$ grid\n", "- 144 RTC neurons arranged in a $12 \\times 12$ grid\n", "- 64 IN neurons arranged in an $8 \\times 8$ grid\n", "- 100 RE neurons arranged in a $10 \\times 10$ grid\n", "\n", "All neuronal populations are placed on a two-dimensional grid.\n", "\n", "The network connectivity is illustrated as follows:\n", "\n", "![Network Connectivity Structure](../_static/image/braincell_HH_thalamus_oscillations_structure.png)\n", "\n", "Details:\n", "\n", "![Annotations](../_static/image/braincell_HH_thalamus_oscillations_structure_2.png)\n", "\n", "Noise in this example is implemented as Poisson input at 100 Hz, modeled as `noise2HTC`, `noise2RTC`, `noise2IN`, and `noise2RE`.\n", "\n", "HTC and RE neurons have intrinsic gap junction connections, modeled as `gj_HTC` and `gj_RE`, while RTC neurons are also connected to HTC neurons via gap junctions, modeled as `gj_RTC2HTC`.\n", "\n", "Synaptic currents mediated by AMPA receptors are modeled as `HTC2IN_ampa`, `HTC2RE_ampa`, and `RTC2RE_ampa`.\n", "\n", "Synaptic currents mediated by NMDA receptors are modeled as `HTC2IN_nmda`, `HTC2RE_nmda`, and `RTC2RE_nmda`.\n", "\n", "Synaptic currents mediated by GABA receptors are modeled as `IN2RTC`, `RE2RTC`, `RE2HTC`, `RE2IN`, and `RE2RE`.\n", "\n", "Combining the relevant parameters from the reference, the resulting network model is as follows:\n" ] }, { "cell_type": "code", "execution_count": 181, "id": "6d5e08a6da906173", "metadata": { "ExecuteTime": { "end_time": "2025-09-22T07:33:52.738185Z", "start_time": "2025-09-22T07:33:52.724770Z" } }, "outputs": [], "source": [ "class Thalamus(brainstate.nn.Module):\n", " def __init__(\n", " self,\n", " g_input: Dict[str, float],\n", " g_KL: Dict[str, float],\n", " HTC_V_init: Callable = braintools.init.Constant(-65. * u.mV),\n", " RTC_V_init: Callable = braintools.init.Constant(-65. * u.mV),\n", " IN_V_init: Callable = braintools.init.Constant(-70. * u.mV),\n", " RE_V_init: Callable = braintools.init.Constant(-70. * u.mV),\n", " ):\n", " super(Thalamus, self).__init__()\n", "\n", " # populations\n", " self.HTC = HTC(7 * 7, gKL=g_KL['TC'], V_initializer=HTC_V_init)\n", " self.RTC = RTC(12 * 12, gKL=g_KL['TC'], V_initializer=RTC_V_init)\n", " self.RE = TRN(10 * 10, gKL=g_KL['RE'], V_initializer=IN_V_init)\n", " self.IN = IN(8 * 8, gKL=g_KL['IN'], V_initializer=RE_V_init)\n", "\n", " # noises\n", " self.noise2HTC = brainpy.state.AlignPostProj(\n", " brainpy.state.PoissonSpike(self.HTC.varshape, freqs=100 * u.Hz),\n", " comm=brainstate.nn.OneToOne(self.HTC.varshape, g_input['TC'], ),\n", " syn=brainpy.state.Expon.desc(self.HTC.varshape, tau=5. * u.ms),\n", " out=brainpy.state.COBA.desc(E=0. * u.mV),\n", " post=self.HTC,\n", " )\n", " self.noise2RTC = brainpy.state.AlignPostProj(\n", " brainpy.state.PoissonSpike(self.RTC.varshape, freqs=100 * u.Hz),\n", " comm=brainstate.nn.OneToOne(self.RTC.varshape, g_input['TC']),\n", " syn=brainpy.state.Expon.desc(self.RTC.varshape, tau=5. * u.ms),\n", " out=brainpy.state.COBA.desc(E=0. * u.mV),\n", " post=self.RTC,\n", " )\n", " self.noise2IN = brainpy.state.AlignPostProj(\n", " brainpy.state.PoissonSpike(self.IN.varshape, freqs=100 * u.Hz),\n", " comm=brainstate.nn.OneToOne(self.IN.varshape, g_input['IN']),\n", " syn=brainpy.state.Expon.desc(self.IN.varshape, tau=5. * u.ms),\n", " out=brainpy.state.COBA.desc(E=0. * u.mV),\n", " post=self.IN,\n", " )\n", " self.noise2RE = brainpy.state.AlignPostProj(\n", " brainpy.state.PoissonSpike(self.RE.varshape, freqs=100 * u.Hz),\n", " comm=brainstate.nn.OneToOne(self.RE.varshape, g_input['RE']),\n", " syn=brainpy.state.Expon.desc(self.RE.varshape, tau=5. * u.ms),\n", " out=brainpy.state.COBA.desc(E=0. * u.mV),\n", " post=self.RE,\n", " )\n", "\n", " # HTC cells were connected with gap junctions\n", " self.gj_HTC = brainpy.state.SymmetryGapJunction(\n", " self.HTC, 'V', conn=ProbDist(dist=2., prob=0.3), weight=1e-2 * u.siemens\n", " )\n", "\n", " # HTC provides feedforward excitation to INs\n", " self.HTC2IN_ampa = brainpy.state.CurrentProj(\n", " self.HTC.align_pre(\n", " brainpy.state.STD.desc(tau=700 * u.ms, U=0.07)\n", " ).align_pre(\n", " brainpy.state.AMPA(self.HTC.varshape, alpha=0.94 / (u.ms * u.mM), beta=0.18 / u.ms)\n", " ).prefetch('g'),\n", " comm=brainstate.nn.FixedNumConn(self.HTC.varshape, self.IN.varshape, 0.3, 6e-3),\n", " out=brainpy.state.COBA(E=0. * u.mV),\n", " post=self.IN,\n", " )\n", " self.HTC2IN_nmda = brainpy.state.CurrentProj(\n", " self.HTC.align_pre(\n", " brainpy.state.STD.desc(tau=700 * u.ms, U=0.07),\n", " ).align_pre(\n", " brainpy.state.AMPA(self.HTC.varshape, alpha=1.0 / (u.ms * u.mM), beta=0.0067 / u.ms)\n", " ).prefetch('g'),\n", " comm=brainstate.nn.FixedNumConn(self.HTC.varshape, self.IN.varshape, 0.3, 3e-3),\n", " out=MgBlock(),\n", " post=self.IN,\n", " )\n", "\n", " # INs delivered feedforward inhibition to RTC cells\n", " self.IN2RTC = brainpy.state.CurrentProj(\n", " self.IN.align_pre(\n", " brainpy.state.STD.desc(tau=700 * u.ms, U=0.07)\n", " ).align_pre(\n", " brainpy.state.GABAa(self.IN.varshape, alpha=10.5 / (u.ms * u.mM), beta=0.166 / u.ms)\n", " ).prefetch_delay('g', 2 * u.ms),\n", " comm=brainstate.nn.FixedNumConn(self.IN.varshape, self.RTC.varshape, 0.3, 3e-3),\n", " out=brainpy.state.COBA(E=-80. * u.mV),\n", " post=self.RTC,\n", " )\n", "\n", " # 20% RTC cells electrically connected with HTC cells\n", " self.gj_RTC2HTC = brainpy.state.SymmetryGapJunction(\n", " (self.RTC, self.HTC), 'V', conn=ProbDist(dist=2., prob=0.3, pre_ratio=0.2), weight=1 / 300 * u.mS\n", " )\n", "\n", " # Both HTC and RTC cells sent glutamatergic synapses to RE neurons, while\n", " # receiving GABAergic feedback inhibition from the RE population\n", " self.HTC2RE_ampa = brainpy.state.CurrentProj(\n", " self.HTC.align_pre(\n", " brainpy.state.STD.desc(tau=700 * u.ms, U=0.07)\n", " ).align_pre(\n", " brainpy.state.AMPA(self.HTC.varshape, alpha=0.94 / (u.ms * u.mM), beta=0.18 / u.ms)\n", " ).prefetch_delay('g', 2 * u.ms),\n", " comm=brainstate.nn.FixedNumConn(self.HTC.varshape, self.RE.varshape, 0.2, 4e-3),\n", " out=brainpy.state.COBA(E=0. * u.mV),\n", " post=self.RE,\n", " )\n", " self.RTC2RE_ampa = brainpy.state.CurrentProj(\n", " self.RTC.align_pre(\n", " brainpy.state.STD.desc(tau=700 * u.ms, U=0.07)\n", " ).align_pre(\n", " brainpy.state.AMPA(self.RTC.varshape, alpha=0.94 / (u.ms * u.mM), beta=0.18 / u.ms)\n", " ).prefetch_delay('g', 2 * u.ms),\n", " comm=brainstate.nn.FixedNumConn(self.RTC.varshape, self.RE.varshape, 0.2, 4e-3),\n", " out=brainpy.state.COBA(E=0. * u.mV),\n", " post=self.RE,\n", " )\n", " self.HTC2RE_nmda = brainpy.state.CurrentProj(\n", " self.HTC.align_pre(\n", " brainpy.state.STD.desc(tau=700 * u.ms, U=0.07)\n", " ).align_pre(\n", " brainpy.state.AMPA(self.HTC.varshape, alpha=1. / (u.ms * u.mM), beta=0.0067 / u.ms)\n", " ).prefetch_delay('g', 2 * u.ms),\n", " comm=brainstate.nn.FixedNumConn(self.HTC.varshape, self.RE.varshape, 0.2, 2e-3),\n", " out=MgBlock(),\n", " post=self.RE,\n", " )\n", " self.RTC2RE_nmda = brainpy.state.CurrentProj(\n", " self.RTC.align_pre(\n", " brainpy.state.STD.desc(tau=700 * u.ms, U=0.07)\n", " ).align_pre(\n", " brainpy.state.AMPA(self.RTC.varshape, alpha=1. / (u.ms * u.mM), beta=0.0067 / u.ms)\n", " ).prefetch_delay('g', 2 * u.ms),\n", " comm=brainstate.nn.FixedNumConn(self.RTC.varshape, self.RE.varshape, 0.2, 2e-3),\n", " out=MgBlock(),\n", " post=self.RE,\n", " )\n", " self.RE2HTC = brainpy.state.CurrentProj(\n", " self.RE.align_pre(\n", " brainpy.state.STD.desc(tau=700 * u.ms, U=0.07)\n", " ).align_pre(\n", " brainpy.state.GABAa(self.RE.varshape, alpha=10.5 / (u.ms * u.mM), beta=0.166 / u.ms)\n", " ).prefetch_delay('g', 2 * u.ms),\n", " comm=brainstate.nn.FixedNumConn(self.RE.varshape, self.HTC.varshape, 0.2, 3e-3),\n", " out=brainpy.state.COBA(E=-80. * u.mV),\n", " post=self.HTC,\n", " )\n", " self.RE2RTC = brainpy.state.CurrentProj(\n", " self.RE.align_pre(\n", " brainpy.state.STD.desc(tau=700 * u.ms, U=0.07)\n", " ).align_pre(\n", " brainpy.state.GABAa(self.RE.varshape, alpha=10.5 / (u.ms * u.mM), beta=0.166 / u.ms)\n", " ).prefetch_delay('g', 2 * u.ms),\n", " comm=brainstate.nn.FixedNumConn(self.RE.varshape, self.RTC.varshape, 0.2, 3e-3),\n", " out=brainpy.state.COBA(E=-80. * u.mV),\n", " post=self.RTC,\n", " )\n", "\n", " # RE neurons were connected with both gap junctions and GABAergic synapses\n", " self.gj_RE = brainpy.state.SymmetryGapJunction(\n", " self.RE, 'V', conn=ProbDist(dist=2., prob=0.3, pre_ratio=0.2), weight=1 / 300 * u.mS\n", " )\n", " self.RE2RE = brainpy.state.CurrentProj(\n", " self.RE.align_pre(\n", " brainpy.state.STD.desc(tau=700 * u.ms, U=0.07)\n", " ).align_pre(\n", " brainpy.state.GABAa(self.RE.varshape, alpha=10.5 / (u.ms * u.mM), beta=0.166 / u.ms)\n", " ).prefetch_delay('g', 2 * u.ms),\n", " comm=brainstate.nn.FixedNumConn(self.RE.varshape, self.RE.varshape, 0.2, 1e-3),\n", " out=brainpy.state.COBA(E=-70. * u.mV),\n", " post=self.RE,\n", " )\n", "\n", " # 10% RE neurons project GABAergic synapses to local interneurons\n", " # probability (0.05) was used for the RE->IN synapses according to experimental data\n", " self.RE2IN = brainpy.state.CurrentProj(\n", " self.RE.align_pre(\n", " brainpy.state.STD.desc(tau=700 * u.ms, U=0.07)\n", " ).align_pre(\n", " brainpy.state.GABAa(self.RE.varshape, alpha=10.5 / (u.ms * u.mM), beta=0.166 / u.ms)\n", " ).prefetch_delay('g', 2 * u.ms),\n", " comm=brainstate.nn.FixedNumConn(self.RE.varshape, self.IN.varshape, 0.05, 1e-3, afferent_ratio=0.1),\n", " out=brainpy.state.COBA(E=-80. * u.mV),\n", " post=self.IN,\n", " )\n", "\n", " def update(self, i, t, current):\n", " with brainstate.environ.context(t=t, i=i):\n", " self.noise2HTC()\n", " self.noise2RTC()\n", " self.noise2IN()\n", " self.noise2RE()\n", "\n", " self.HTC2IN_ampa()\n", " self.HTC2IN_nmda()\n", "\n", " self.IN2RTC()\n", "\n", " self.HTC2RE_ampa()\n", " self.RTC2RE_ampa()\n", " self.HTC2RE_nmda()\n", " self.RTC2RE_nmda()\n", "\n", " self.RE2HTC()\n", " self.RE2RTC()\n", " self.RE2RE()\n", " self.RE2IN()\n", "\n", " self.gj_HTC()\n", " self.gj_RTC2HTC()\n", " self.gj_RE()\n", "\n", " htc_spike = self.HTC(current)\n", " rtc_spike = self.RTC(current)\n", " re_spike = self.RE(current)\n", " in_spike = self.IN(current)\n", "\n", " return {\n", " 'HTC.V': self.HTC.V.value,\n", " 'RTC.V': self.RTC.V.value,\n", " 'IN.V': self.IN.V.value,\n", " 'RE.V': self.RE.V.value,\n", " 'HTC.spike': htc_spike,\n", " 'RTC.spike': rtc_spike,\n", " 'RE.spike': re_spike,\n", " 'IN.spike': in_spike,\n", " }" ] }, { "cell_type": "markdown", "id": "3e9210e97320f3e2", "metadata": {}, "source": [ "For this thalamic network, different parameters can generate oscillations at different frequencies.\n", "\n", "The correspondence between the relevant parameters and oscillation frequencies is shown below:\n", "\n", "![Parameter–Frequency Mapping](../_static/image/braincell_HH_thalamus_oscillations_parameters.png)\n", "\n", "To facilitate simulation, we directly package these parameters as follows:" ] }, { "cell_type": "code", "execution_count": 182, "id": "d29da2c5fd6d7064", "metadata": { "ExecuteTime": { "end_time": "2025-09-22T03:31:52.249622Z", "start_time": "2025-09-22T03:31:52.245670Z" } }, "outputs": [], "source": [ "states = {\n", " 'delta': dict(\n", " g_input={'IN': 1e-4 * u.mS, 'RE': 1e-4 * u.mS, 'TC': 1e-4 * u.mS},\n", " g_KL={'TC': 0.035 * u.mS / u.cm ** 2, 'RE': 0.03 * u.mS / u.cm ** 2, 'IN': 0.01 * u.mS / u.cm ** 2}\n", " ),\n", " 'spindle': dict(\n", " g_input={'IN': 3e-4 * u.mS, 'RE': 3e-4 * u.mS, 'TC': 3e-4 * u.mS},\n", " g_KL={'TC': 0.01 * u.mS / u.cm ** 2, 'RE': 0.02 * u.mS / u.cm ** 2, 'IN': 0.015 * u.mS / u.cm ** 2}\n", " ),\n", " 'alpha': dict(\n", " g_input={'IN': 1.5e-3 * u.mS, 'RE': 1.5e-3 * u.mS, 'TC': 1.5e-3 * u.mS},\n", " g_KL={'TC': 0. * u.mS / u.cm ** 2, 'RE': 0.01 * u.mS / u.cm ** 2, 'IN': 0.02 * u.mS / u.cm ** 2}\n", " ),\n", " 'gamma': dict(\n", " g_input={'IN': 1.5e-3 * u.mS, 'RE': 1.5e-3 * u.mS, 'TC': 1.7e-2 * u.mS},\n", " g_KL={'TC': 0. * u.mS / u.cm ** 2, 'RE': 0.01 * u.mS / u.cm ** 2, 'IN': 0.02 * u.mS / u.cm ** 2}\n", " ),\n", "}" ] }, { "cell_type": "markdown", "id": "8f479cfd02ca79e6", "metadata": {}, "source": [ "To simulate periodic stimulation of neurons in the network, we define a function `rhythm_const_input`, which generates rhythmic square-wave current inputs." ] }, { "cell_type": "code", "execution_count": 183, "id": "78d933d6785efcf3", "metadata": { "ExecuteTime": { "end_time": "2025-09-22T03:31:52.265785Z", "start_time": "2025-09-22T03:31:52.262633Z" } }, "outputs": [], "source": [ "def rhythm_const_input(amp, freq, length, duration, t_start=0., t_end=None, dt=None):\n", " if t_end is None:\n", " t_end = duration\n", " if length > duration:\n", " raise ValueError(f'Expected length <= duration, while we got {length} > {duration}')\n", " sec_length = 1 / freq\n", " values, durations = [0. * u.mA], [t_start]\n", " for t in u.math.arange(t_start, t_end, sec_length):\n", " values.append(amp)\n", " if t + length <= t_end:\n", " durations.append(length)\n", " values.append(0. * u.mA)\n", " if t + sec_length <= t_end:\n", " durations.append(sec_length - length)\n", " else:\n", " durations.append(t_end - t - length)\n", " else:\n", " durations.append(t_end - t)\n", " values.append(0. * u.mA)\n", " durations.append(duration - t_end)\n", " return braintools.input.section(values=values, durations=durations)" ] }, { "cell_type": "markdown", "id": "319cb0a5da72c8a3", "metadata": {}, "source": [ "For convenient plotting, we encapsulate a small utility function to draw time-series curves." ] }, { "cell_type": "code", "execution_count": 184, "id": "2bf9c9d0d1baf331", "metadata": { "ExecuteTime": { "end_time": "2025-09-22T03:31:52.285123Z", "start_time": "2025-09-22T03:31:52.282836Z" } }, "outputs": [], "source": [ "def line_plot(ax, xs, ys, ylabel=None, xlim=None):\n", " ax.plot(xs, ys)\n", " ax.set_xlabel('Time (ms)')\n", " if ylabel is not None:\n", " ax.set_ylabel(ylabel)\n", " if xlim is not None:\n", " ax.set_xlim(xlim)" ] }, { "cell_type": "markdown", "id": "10a0ce4641d6ccf6", "metadata": {}, "source": [ "Similarly, we encapsulate a function to plot spike raster diagrams of neurons." ] }, { "cell_type": "code", "execution_count": 185, "id": "9ee7e65430da50", "metadata": { "ExecuteTime": { "end_time": "2025-09-22T03:31:52.301567Z", "start_time": "2025-09-22T03:31:52.298979Z" } }, "outputs": [], "source": [ "def raster_plot(\n", " ax, times, spikes, ylabel='Neuron Index', xlabel='Time [ms]',\n", " xlim=None, marker='.', markersize=2, color='k'\n", "):\n", " elements = np.where(spikes > 0.)\n", " index = elements[1]\n", " time = times[elements[0]]\n", " ax.plot(time, index, marker + color, markersize=markersize)\n", " if xlabel:\n", " ax.set_xlabel(xlabel)\n", " if ylabel:\n", " ax.set_ylabel(ylabel)\n", " if xlim:\n", " ax.set_xlim(xlim)" ] }, { "cell_type": "markdown", "id": "9da31a29beb40d5", "metadata": {}, "source": [ "Simulate the network." ] }, { "cell_type": "code", "execution_count": 186, "id": "8e11bdda9e775e25", "metadata": {}, "outputs": [], "source": [ "def try_network(state='delta'):\n", " net = Thalamus(\n", " IN_V_init=braintools.init.Constant(-70. * u.mV),\n", " RE_V_init=braintools.init.Constant(-70. * u.mV),\n", " HTC_V_init=braintools.init.Constant(-80. * u.mV),\n", " RTC_V_init=braintools.init.Constant(-80. * u.mV),\n", " **states[state],\n", " )\n", " brainstate.nn.init_all_states(net)\n", "\n", " duration = 3e3 * u.ms # 3 seconds\n", " currents = rhythm_const_input(\n", " 2e-4 * u.mA,\n", " freq=4. * u.Hz,\n", " length=10. * u.ms,\n", " duration=duration,\n", " t_end=2e3 * u.ms,\n", " t_start=1e3 * u.ms\n", " )\n", " indices = np.arange(currents.shape[0])\n", " times = indices * brainstate.environ.get_dt()\n", " mon = brainstate.transform.for_loop(net.update, indices, times, currents, pbar=200)\n", "\n", " fig, gs = braintools.visualize.get_figure(5, 2, 2, 5)\n", " line_plot(fig.add_subplot(gs[0, :]), times, currents, ylabel='Current', xlim=(0, duration / u.ms))\n", " line_plot(fig.add_subplot(gs[1, 0]), times, mon.get('HTC.V'), ylabel='HTC', xlim=(0, duration / u.ms))\n", " line_plot(fig.add_subplot(gs[2, 0]), times, mon.get('RTC.V'), ylabel='RTC', xlim=(0, duration / u.ms))\n", " line_plot(fig.add_subplot(gs[3, 0]), times, mon.get('IN.V'), ylabel='IN', xlim=(0, duration / u.ms))\n", " line_plot(fig.add_subplot(gs[4, 0]), times, mon.get('RE.V'), ylabel='RE', xlim=(0, duration / u.ms))\n", " raster_plot(fig.add_subplot(gs[1, 1]), times, mon.get('HTC.spike'), xlim=(0, duration / u.ms))\n", " raster_plot(fig.add_subplot(gs[2, 1]), times, mon.get('RTC.spike'), xlim=(0, duration / u.ms))\n", " raster_plot(fig.add_subplot(gs[3, 1]), times, mon.get('IN.spike'), xlim=(0, duration / u.ms))\n", " raster_plot(fig.add_subplot(gs[4, 1]), times, mon.get('RE.spike'), xlim=(0, duration / u.ms))\n", " plt.suptitle(f'Thalamus Network State: {state}')\n", " plt.show()" ] }, { "cell_type": "markdown", "id": "ba41962ce9c94deb", "metadata": {}, "source": [ "Here, we run the simulation using parameters that generate `alpha` and `gamma` oscillations as examples." ] }, { "cell_type": "code", "execution_count": 187, "id": "b6cc2fef8d85cac2", "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "e3f65cdb8e4f4682bd55b40b893f87c2", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/30000 [00:00" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "433f851cd1934d368be21b0b981b4c93", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/30000 [00:00" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "if __name__ == '__main__':\n", " try_network('alpha')\n", " try_network('gamma')" ] }, { "cell_type": "markdown", "id": "60dafad6bdfd1b69", "metadata": {}, "source": [ "Observing the results, the constructed thalamic network model successfully exhibits the corresponding oscillatory rhythms." ] }, { "cell_type": "markdown", "id": "4561b83096fc6a2b", "metadata": {}, "source": [ "The example models in this article are adapted from:\n", "\n", "- G, Li., Henriquez, C. S., & Fröhlich, F. (2017), Unified thalamic model generates multiple distinct oscillations with state-dependent entrainment by stimulation., PLoS Comput. Biol., 13, 10, e1005797" ] } ], "metadata": { "kernelspec": { "display_name": "Ecosystem-py", "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.11.13" } }, "nbformat": 4, "nbformat_minor": 5 }