{ "cells": [ { "cell_type": "markdown", "id": "fed475d53adb2402", "metadata": {}, "source": [ "# 使用 `braincell` 对 Hodgkin-Huxley 型丘脑网络进行仿真\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": [ "本节介绍使用 `braincell` 对 Hodgkin-Huxley 型丘脑网络进行仿真。\n", "\n", "这是一个实战仿真示例，在学习使用 `braincell` 对 HH 型的丘脑网络进行建模之前，你应该已经学会了如何对 HH 型的[E - I 网络](braincell_HH_network-zh.ipynb)进行建模。" ] }, { "cell_type": "markdown", "id": "1c15ce72366e2270", "metadata": {}, "source": [ "## 丘脑网络概述\n", "\n", "丘脑在丘脑皮质振荡的产生中发挥着很关键的作用，但是其背后的机制尚不明确。为了探究独立的丘脑能否生成多种不同的振荡，我们需要对其进行建模。\n", "\n", "多种神经递质能对丘脑进行调控，其中乙酰胆碱（ACh）和去甲肾上腺素（NE）是塑造振荡的关键。\n", "\n", "我们准备构建的丘脑网络模型由两个核团构成，分别是外侧膝状体核（LGN）和丘脑网状核（TRN）。其中 LGN 包含三种神经元，分别是高阈值爆发性丘脑皮质细胞（HTC）、中继模式丘脑皮质细胞（RTC）、局部中间神经元（IN）；而 TRN 则包含一种神经元，为网状细胞（RE）。\n", "\n", "这些不同神经元通过 gap junction 和突触形成网络，在 ACh 和 NE 的调控下，生成 delta 、spindle 、 alpha 、 gamma 这四种不同频段的网络。" ] }, { "cell_type": "markdown", "id": "b250343a22c8f7ef", "metadata": {}, "source": [ "## 丘脑神经元建模\n", "\n", "在了解丘脑网络连接结构后，我们首先对丘脑神经元进行建模。" ] }, { "cell_type": "code", "execution_count": 137, "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 brainpy.state\n", "import braincell\n", "import brainstate\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": [ "首先，我们先建立一个丘脑神经元的基类，其继承自 `SingleCompartment` 。" ] }, { "cell_type": "code", "execution_count": 138, "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": [ "本例中的神经元与[E - I 网络中的 HH 神经元](braincell_HH_neuron-zh.ipynb)建模有一些区别，相较之下本例显式的写出了更新膜电位动力学的接口 `compute_derivative` 。\n", "\n", "此接口整合了外部输入电流并加上所有例子通道的电流，然后更新膜电位与离子通道的的导数，是典型的 HH 型神经元的微分方程更新方法。\n", "\n", "同时设置 `step_run` ，便于后续仿真与可视化。" ] }, { "cell_type": "markdown", "id": "9aa545188c69fa84", "metadata": {}, "source": [ "接下来，我们来对不同的四种神经元进行建模。\n", "\n", "此处我们以 HTC 神经元为例。\n", "\n", "HTC 神经元继承自上面建立的 ThalamusNeuron 类，按照 ThalamusNeuron 的方式来更新膜电位动力学。\n", "\n", "依照参考文献，其膜面积为 $2.9 \\times 10^{-4} \\ \\text{cm}^2$ 。\n", "\n", "同时，HTC 神经元拥许多种不同的离子电流：\n", "\n", "- 快速钠电流 ($I_\\text{Na}$)\n", "- 延迟整流钾电流 ($I_\\text{DR}$)\n", "- 超极化激活阳离子电流 ($I_\\text{H}$)\n", "- 高阈值 L 型钙电流 ($I_{\\text{Ca/L}}$)\n", "- 钙激活非选择性阳离子电流 ($I_\\text{CAN}$)\n", "- 低阈值 T 型钙电流 ($I_{\\text{Ca/T}}$)\n", "- 高阈值 T 型钙电流 ($I_{\\text{Ca/HT}}$)\n", "- 钙依赖性钾电流 ($I_\\text{AHP}$)\n", "- 钾漏电流 ($I_\\text{KL}$)\n", "- 漏电流 ($I_\\text{L}$)\n", "\n", "\n", "给 HTC 模型加入这些离子通道，并依据参考文献确定相关参数。" ] }, { "cell_type": "code", "execution_count": 139, "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": [ "同理，对 RTC 、 IN 、 TRN 神经元进行建模。" ] }, { "cell_type": "code", "execution_count": 140, "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": 141, "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": 142, "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": [ "完成建模后，对建立好的单神经元进行仿真。此处以 HTC 神经元为例。" ] }, { "cell_type": "code", "execution_count": 143, "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.0922 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": [ "## 丘脑网络建模\n", "\n", "在完成对神经元的建模后，我们开始对完整的丘脑网络进行建模。\n", "\n", "目前神经科学领域的大部分常见突触连接模型都可以在 `braincell` 中找到，但在此例中，出现了相对比较少见的 NMDA 突触的镁离子阻塞效应，我们需要手动对其建模。\n", "\n", "我们要建立的 `MgBlock` 继承自 `SynOut` ，表示这是一个突触输出模型。\n", "\n", "根据其公式：\n", "\n", "$$I_{\\text{NMDA}}(V) = \\dfrac{g \\, (E - V)}{1 + \\exp \\!\\left( -\\dfrac{V + 25}{12.5} \\right)}$$\n", "\n", "我们可以轻松对其建模。" ] }, { "cell_type": "code", "execution_count": 144, "id": "c1d2d018c1a967a6", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "C:\\Users\\Administrator\\AppData\\Local\\Temp\\ipykernel_33864\\2860468575.py:1: DeprecationWarning: 'brainstate.nn.SynOut' is deprecated and will be removed in a future version. Please use 'brainpy.state.SynOut' instead.\n", " class MgBlock(brainstate.nn.SynOut):\n" ] } ], "source": [ "class MgBlock(brainstate.nn.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": [ "同时，由于本例中 gap junction 的连接是基于空间距离和概率的，因此我们在此建立一个 `ProbDist` 。\n", "\n", "此 `ProbDist` 类的作用是，给定前神经元群体和后神经元群体的规模，会根据以下规则生成一组突触连接：\n", "\n", "- 只连接在二维平面中空间距离不超过 `dist` 的神经元\n", "- 每对神经元有 `prob` 的概率被连接\n", "- 每个前神经元有 `pre_ratio` 的概率参与连接\n", "- `include_self` 控制是否允许自连接\n", "\n", "基于以上规则，我们之后在进行网络建模时可以很具体精细的对 gap junction 进行建模。" ] }, { "cell_type": "code", "execution_count": 145, "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": [ "在这个丘脑网络中，神经元按照比例排列分布：\n", "- 49 个 HTC 神经元，排列成 $7 \\times 7$ 的网格\n", "- 144 个 RTC 神经元，排列成 $12 \\times 12$ 的网格\n", "- 64 个 IN 神经元，排列成 $8 \\times 8$ 的网格\n", "- 100 个 RE 神经元，排列成 $10 \\times 10$ 的网格\n", "\n", "所有神经元群体都被放置在一个二维网格中。\n", "\n", "其网络连接如下：\n", "\n", "![网络连接结构](../_static/image/braincell_HH_thalamus_oscillations_structure.png)\n", "\n", "其中：\n", "\n", "![注释](../_static/image/braincell_HH_thalamus_oscillations_structure_2.png)\n", "\n", "噪声在本例中为频率为 100 Hz 的泊松输入，建模为 `noise2HTC` 、 `noise2RTC` 、 `noise2IN` 、 `noise2RE` 。\n", "\n", "HTC 神经元和 RE 神经元拥有自身的 gap junction 连接，建模为 `gj_HTC` 和 `gj_RE` ，同时 RTC 神经元与 HTC 神经元间也拥有 gap junction 连接，建模为 `gj_RTC2HTC` 。\n", "\n", "由 AMPA 受体介导的突触电流被建模为 `HTC2IN_ampa` 、 `HTC2RE_ampa` 、 `RTC2RE_ampa` 。\n", "\n", "由 NMDA 受体介导的突触电流被建模为 `HTC2IN_nmda` 、 `HTC2RE_nmda` 、 `RTC2RE_nmda` 。\n", "\n", "由 GABA 受体介导的突触电流被建模为 `IN2RTC` 、 `RE2RTC` 、 `RE2HTC` 、 `RE2IN` 、 `RE2RE` 。\n", "\n", "结合参考文献中的相关参数进行建模，网络模型如下：" ] }, { "cell_type": "code", "execution_count": 146, "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": [ "对于这个丘脑网络，不同的参数可以使其产生不同频率的振荡。\n", "\n", "相关参数与振荡频率的对应关系如下：\n", "\n", "![参数对应表](../_static/image/braincell_HH_thalamus_oscillations_parameters.png)\n", "\n", "为便于模拟仿真，我们直接将参数包装好：" ] }, { "cell_type": "code", "execution_count": 147, "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": [ "为了在神经网络中模拟给神经元注入周期性刺激，我们需要定义一个函数 `rhythm_const_input` ，进而可以生成由节律的方波电流输入。" ] }, { "cell_type": "code", "execution_count": 148, "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": [ "为便于作图，封装好一个画时间序列曲线的小工具函数。" ] }, { "cell_type": "code", "execution_count": 149, "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": [ "同样，封装好一个画神经元脉冲光栅图的绘制函数。" ] }, { "cell_type": "code", "execution_count": 150, "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": [ "对网络进行模拟。" ] }, { "cell_type": "code", "execution_count": 151, "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": [ "此处我们以产生 `alpha` 和 `gamma` 振荡的参数为例，运行仿真程序。" ] }, { "cell_type": "code", "execution_count": 152, "id": "b6cc2fef8d85cac2", "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "2f68a92ece0b4397afe588b85f350d0d", "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": "c41be24c65254872909c77bb4ad84885", "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": [ "观察运行结果，我们构建的丘脑神经网络模型成功的展现了相关的振荡节律。" ] }, { "cell_type": "markdown", "id": "4561b83096fc6a2b", "metadata": {}, "source": [ "本文中举例的模型来自：\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\n" ] } ], "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 }