{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "24e3c9251f9e0754",
   "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_neuron-zh.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_neuron-zh.ipynb)\n",
    "\n",
    "本节介绍使用 `braincell` 对 Hodgkin-Huxley 神经元仿真。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8121ee574361d5d",
   "metadata": {},
   "source": [
    "## HH 模型概述\n",
    "\n",
    "艾伦 · 霍奇金（Allen Hodgkin， 1914 - 1998）和安德鲁 · 赫胥黎（Andrew Huxley， 1917 - 2012）首次描述了动作电位的产生机制。 在 1950 年左右两人通过一系列的实验阐明了动作电位的离子机制，并于 1963 年凭借动作电位机制的阐述获得了诺贝尔生理学或医学奖。 霍奇金 - 赫胥黎模型（Hodgkin - Huxley model， HH 模型）是第一个精准描述神经元动作电位产生过程的数学模型。\n",
    "\n",
    "对于一个 HH 型神经元，我们主要关注它的细胞膜。神经元细胞膜是分隔细胞内外的磷脂双分子层结构，在细胞质或胞外的水性环境下，水溶性的带电分子，如离子和神经递质等，都不能被动地穿过细胞膜的磷脂双分子层，因此需要运输蛋白将它们从一侧运至另一侧。离子通道就是一种跨膜运输蛋白。\n",
    "\n",
    "有些离子通道是由闸门控制的，可以在电、化学或物理刺激下开放或关闭。简单来说，对应的离子通道允许对应的离子通过，如钠离子通道允许钠离子通过，钾离子通道允许钾离子通过。\n",
    "\n",
    "有些离子通道是无闸门的，这些通道始终开放，没有闸门控制。HH 模型将这些通道合并成一个等效的漏电流通道。\n",
    "\n",
    "HH 型神经元细胞膜的典型结构如下图：\n",
    "\n",
    "![层级结构](_static/image/braincell_HH_neuron_structure.png)\n",
    "\n",
    "对 HH 型神经元的机制进行总结：\n",
    "- 神经元细胞由细胞膜包围，膜上分布着各种离子通道\n",
    "- 离子通道控制离子的进出\n",
    "- 离子的流动会引起细胞膜内外电位的变化，从而驱动神经元的活动\n",
    "\n",
    "为了描述细胞膜电位的变化以及动作电位的产生，一种有效的描述方法是构建出神经元的等效电路。\n",
    "通常这种等效电路模型包含三个部分：\n",
    "- 代表离子通道的电阻\n",
    "- 代表离子浓度梯度的电源\n",
    "- 代表细胞膜储存电荷能力的电容\n",
    "\n",
    "神经元各部分与等效电路组件的对应关系如下图：\n",
    "\n",
    "![等效电路组件](_static/image/braincell_HH_neuron_part.png)\n",
    "\n",
    "细胞膜电容由 $ C_m $ 表示，离子通道可等效为一个包含电阻 $ R_X $ 和电源 $ E_X $ 的电路元件。\n",
    "由于电阻是电导的倒数，因此离子通道的电导为 $ g_X = 1 / R_X $ 。此处的 $ X $ 指代不同的离子，如 $ Na $ 、 $ K $ 等。\n",
    "\n",
    "结合以上，HH 模型的等效电路如下图：\n",
    "\n",
    "![等效电路](_static/image/braincell_HH_neuron_circuit.png)\n",
    "\n",
    "电容器存储的电荷量 $ q $ 与膜电位 $ V_m $ 的关系可表示为 $ q = C_m V_m $ 。\n",
    "\n",
    "掌握了以上知识后，我们可以开始对 HH 型神经元建模。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "804bea63109ec209",
   "metadata": {},
   "source": [
    "## 离子通道建模\n",
    "\n",
    "如上述， HH 型神经元的电学特性主要由离子通道控制，所以如果想对神经元甚至网络进行准确的建模，首先要学会如何建模离子通道。\n",
    "\n",
    "在 HH 模型中，每个离子通道都被视为一种跨膜蛋白质，形成一个孔，离子可以通过该孔沿着其浓度梯度扩散。这些孔有许多闸门，可以打开也可以关闭，而每个闸门打开或者关闭的概率取决于膜电位。\n",
    "\n",
    "对于每个闸门，将 $ \\alpha(V) $ 和 $ \\beta(V) $ 分别定义为闸门从关闭状态到打开状态以及从打开状态到关闭状态的电压依赖的速率。令 $ m $ 是闸门打开的比例， $ 1 - m $ 是闸门关闭的比例，则我们可以得到:\n",
    "\n",
    "$$\n",
    "\\frac{dm}{dt} = \\alpha(V)(1 - m) - \\beta(V)m\n",
    "$$\n",
    "\n",
    "$$\n",
    "= \\frac{m_{\\infty}(V) - m}{\\tau(V)},\n",
    "$$\n",
    "\n",
    "其中\n",
    "\n",
    "$$\n",
    "m_{\\infty}(V) = \\frac{\\alpha(V)}{\\alpha(V) + \\beta(V)}, \\quad \\tau(V) = \\frac{1}{\\alpha(V) + \\beta(V)}.\n",
    "$$\n",
    "\n",
    "在 HH 模型中，$ m_{\\infty}(V) $ 和 $ \\tau(V) $ 中的参数 $ \\alpha $、$ \\beta $ 可以根据实验数据拟合计算得到。\n",
    "\n",
    "闸门分为激活闸门和失活闸门。不同的离子通道拥有不同种类和数量的闸门，如典型的钠离子通道：\n",
    "\n",
    "$$\n",
    "g_{\\text{Na}} = \\bar{g}_{\\text{Na}} m^3 h,\n",
    "$$\n",
    "\n",
    "其中 $ \\bar{g}_{\\text{Na}} $ 是钠离子通道的最大电导值，激活门控 $ m \\in [0, 1] $ 和失活门控 $ h \\in [0, 1] $ 是满足一阶动力学方程的变量，$ m^3 h $ 表示钠离子通道具有 3 个独立的激活闸门和 1 个独立的失活闸门，只有每个部位都打开，才能允许 $ \\text{Na}^+ $ 通过。\n",
    "\n",
    "相比之下，典型的钾离子通道只拥有激活闸门而没有失活闸门，漏电流通道则两种闸门都没有。\n",
    "\n",
    "`braincell` 框架中对于离子通道的建模主要就依赖于上式。\n",
    "\n",
    "以典型的钠离子通道为例，其数学形式如下：\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "g_{\\mathrm{Na}} &= g_{\\mathrm{max}} \\cdot p^3 \\cdot q \\\\\n",
    "\\frac{dp}{dt} &= \\phi \\left( \\alpha_p (1-p) - \\beta_p p \\right) \\\\\n",
    "\\frac{dq}{dt} &= \\phi \\left( \\alpha_q (1-h) - \\beta_q h \\right) \\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "其中 $\\phi$ 是温度系数。\n",
    "\n",
    "对其建模："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "1a1a45cc60bad870",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-16T06:37:47.446454Z",
     "start_time": "2025-09-16T06:37:46.537861Z"
    }
   },
   "outputs": [],
   "source": [
    "from typing import Union, Callable, Optional\n",
    "\n",
    "import braincell\n",
    "import braintools\n",
    "import brainstate\n",
    "import brainunit as u\n",
    "\n",
    "from braincell.channel import SodiumChannel\n",
    "from braincell import IonInfo, DiffEqState"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "f14e9e0066cfb758",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-16T06:37:47.467196Z",
     "start_time": "2025-09-16T06:37:47.460148Z"
    }
   },
   "outputs": [],
   "source": [
    "class INa_p3q_markov(SodiumChannel):\n",
    "    def __init__(\n",
    "        self,\n",
    "        size: brainstate.typing.Size,\n",
    "        g_max: Union[brainstate.typing.ArrayLike, Callable] = 90. * (u.mS / u.cm ** 2),\n",
    "        phi: Union[brainstate.typing.ArrayLike, Callable] = 1.,\n",
    "        name: Optional[str] = None,\n",
    "    ):\n",
    "        super().__init__(size=size, name=name, )\n",
    "\n",
    "        # parameters\n",
    "        self.phi = braintools.init.param(phi, self.varshape, allow_none=False)\n",
    "        self.g_max = braintools.init.param(g_max, self.varshape, allow_none=False)\n",
    "\n",
    "    def init_state(self, V, Na: IonInfo, batch_size=None):\n",
    "        self.p = DiffEqState(braintools.init.param(u.math.zeros, self.varshape, batch_size))\n",
    "        self.q = DiffEqState(braintools.init.param(u.math.zeros, self.varshape, batch_size))\n",
    "\n",
    "    def reset_state(self, V, Na: IonInfo, batch_size=None):\n",
    "        alpha = self.f_p_alpha(V)\n",
    "        beta = self.f_p_beta(V)\n",
    "        self.p.value = alpha / (alpha + beta)\n",
    "        alpha = self.f_q_alpha(V)\n",
    "        beta = self.f_q_beta(V)\n",
    "        self.q.value = alpha / (alpha + beta)\n",
    "\n",
    "    def compute_derivative(self, V, Na: IonInfo):\n",
    "        p = self.p.value\n",
    "        q = self.q.value\n",
    "        self.p.derivative = self.phi * (self.f_p_alpha(V) * (1. - p) - self.f_p_beta(V) * p) / u.ms\n",
    "        self.q.derivative = self.phi * (self.f_q_alpha(V) * (1. - q) - self.f_q_beta(V) * q) / u.ms\n",
    "\n",
    "    def current(self, V, Na: IonInfo):\n",
    "        return self.g_max * self.p.value ** 3 * self.q.value * (Na.E - V)\n",
    "\n",
    "    def f_p_alpha(self, V):\n",
    "        raise NotImplementedError\n",
    "\n",
    "    def f_p_beta(self, V):\n",
    "        raise NotImplementedError\n",
    "\n",
    "    def f_q_alpha(self, V):\n",
    "        raise NotImplementedError\n",
    "\n",
    "    def f_q_beta(self, V):\n",
    "        raise NotImplementedError"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce1a270aa5f6a7b3",
   "metadata": {},
   "source": [
    "观察以上代码，不难发现 `braincell` 框架可以很清晰的对离子通道进行建模。\n",
    "\n",
    "- `init_state` 创建需要的状态变量\n",
    "- `reset_state` 构建闸门的方程\n",
    "- `compute_derivative` 求解导数用来更新\n",
    "- `current` 计算通道产生的电流\n",
    "\n",
    "当然，温度系数 $ \\phi $ 也是一个重要的参数。离子通道对温度很敏感，因此开关状态的速率于温度成指数关系，即温度越高，开放关闭的转换速率越快。这里不做进一步讨论。\n",
    "\n",
    "综上，对拥有不同种类和数量闸门的离子通道建模非常便捷。对于这些不同的通道， `init_state` 、 `reset_state` 、 `compute_derivative` 都是相似的，主要的差别就在于 `current` 的数学形式。\n",
    "\n",
    "在建立了这种有 3 个激活闸门和 1 个失活闸门的钠离子通道后，可以根据要实现的具体模型，去确定离子通道的参数。\n",
    "\n",
    "以 1991 年的一篇文献中的具体钠离子通道为例，其数学形式如下：\n",
    "\n",
    "$$\n",
    "\\begin{split}\n",
    "\\begin{aligned}\n",
    "  g_{\\mathrm{Na}} &= g_{\\mathrm{max}} m^3 h \\\\\n",
    "  \\frac {dm} {dt} &= \\phi\\left(\\alpha_m (1-x) - \\beta_m\\right) \\\\\n",
    "  &\\alpha_m(V) = 0.32 \\frac{(13 - V + V_{\\text{sh}})}{\\exp\\left(\\frac{13 - V + V_{\\text{sh}}}{4}\\right) - 1}  \\\\\n",
    "  &\\beta_m(V) = 0.28 \\frac{(V - V_{\\text{sh}} - 40)}{\\exp\\left(\\frac{V - V_{\\text{sh}} - 40}{5}\\right) - 1}  \\\\\n",
    "  \\frac {dh} {dt} &= \\phi\\left(\\alpha_h (1-x) - \\beta_h\\right) \\\\\n",
    "  &\\alpha_h(V) = 0.128 \\cdot \\exp\\left(\\frac{17 - V + V_{\\text{sh}}}{18}\\right)  \\\\\n",
    "  &\\beta_h(V) = \\frac{4}{1 + \\exp\\left(-\\frac{V - V_{\\text{sh}} - 40}{5}\\right)} \\\\\n",
    "\\end{aligned}\n",
    "\\end{split}\n",
    "$$\n",
    "\n",
    "其中 $V_{\\text{sh}}$ 是膜电位偏移，文中值为 -63 mV，$\\phi$ 是温度系数，文中默认值为 1 。\n",
    "\n",
    "对其建模："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "24d0bc9c50fc69f1",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-16T06:37:52.130288Z",
     "start_time": "2025-09-16T06:37:52.126017Z"
    }
   },
   "outputs": [],
   "source": [
    "class INa_TM1991(INa_p3q_markov):\n",
    "    def __init__(\n",
    "        self,\n",
    "        size: brainstate.typing.Size,\n",
    "        g_max: Union[brainstate.typing.ArrayLike, Callable] = 120. * (u.mS / u.cm ** 2),\n",
    "        phi: Union[brainstate.typing.ArrayLike, Callable] = 1.,\n",
    "        V_sh: Union[brainstate.typing.ArrayLike, Callable] = -63. * u.mV,\n",
    "        name: Optional[str] = None,\n",
    "    ):\n",
    "        super().__init__(\n",
    "            size,\n",
    "            name=name,\n",
    "            phi=phi,\n",
    "            g_max=g_max,\n",
    "        )\n",
    "        self.V_sh = braintools.init.param(V_sh, self.varshape, allow_none=False)\n",
    "\n",
    "    def f_p_alpha(self, V):\n",
    "        V = (self.V_sh - V).to_decimal(u.mV)\n",
    "        temp = 13 + V\n",
    "        return 0.32 * 4 / u.math.exprel(temp / 4)\n",
    "\n",
    "    def f_p_beta(self, V):\n",
    "        V = (V - self.V_sh).to_decimal(u.mV)\n",
    "        temp = V - 40\n",
    "        return 0.28 * 5 / u.math.exprel(temp / 5)\n",
    "\n",
    "    def f_q_alpha(self, V):\n",
    "        V = (- V + self.V_sh).to_decimal(u.mV)\n",
    "        return 0.128 * u.math.exp((17 + V) / 18)\n",
    "\n",
    "    def f_q_beta(self, V):\n",
    "        V = (V - self.V_sh).to_decimal(u.mV)\n",
    "        return 4. / (1 + u.math.exp(-(V - 40) / 5))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0a99f891b90b794",
   "metadata": {},
   "source": [
    "观察以上代码，不难发现在原有的钠离子通道基础上进行精细建模非常方便，只需要将对应的数学表达式引入即可。\n",
    "\n",
    "当然，除了自定义具体的离子通道，你也可以使用我们在 `braincell` 中已经内置好的通道，这些通道可以在进行神经元建模的时候直接导入。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "365d88c0ba7bdbd3",
   "metadata": {},
   "source": [
    "## 神经元建模\n",
    "\n",
    "在完成对离子通道的建模后，我们可以很轻松的完成对单个神经元的建模。\n",
    "只需要将对应的离子与通道导入到建立的神经元模型中，并确定参数即可。\n",
    "\n",
    "以经典的 HH 型神经元为例，对其建模："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "a4c59bc1159baa68",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-16T06:37:54.571055Z",
     "start_time": "2025-09-16T06:37:54.567996Z"
    }
   },
   "outputs": [],
   "source": [
    "class HH(braincell.SingleCompartment):\n",
    "    def __init__(self, size, solver='exp_euler', V_th=-20 * u.mV):\n",
    "        super().__init__(size, solver=solver, V_th=V_th, C=1 * u.uF / u.cm ** 2)\n",
    "\n",
    "        self.na = braincell.ion.SodiumFixed(size, E=50. * u.mV)\n",
    "        self.na.add_elem(INa=INa_TM1991(size, g_max=100. * (u.mS / u.cm ** 2), V_sh=-63. * u.mV))\n",
    "\n",
    "        self.k = braincell.ion.PotassiumFixed(size, E=-90 * u.mV)\n",
    "        self.k.add_elem(IK=braincell.channel.IK_TM1991(size, g_max=30. * (u.mS / u.cm ** 2), V_sh=-63. * u.mV))\n",
    "\n",
    "        self.IL = braincell.channel.IL(size,E=-60. * u.mV,g_max=5. * (u.mS / u.cm ** 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "580e0cece4eace46",
   "metadata": {},
   "source": [
    "观察以上代码，不难发现使用 `braincell` 对单神经元的建模非常简单。\n",
    "\n",
    "对其仿真："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "5223cf7648e1dbb4",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-16T06:37:56.677981Z",
     "start_time": "2025-09-16T06:37:56.441028Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "HH(\n",
       "  in_size=(1,),\n",
       "  out_size=(1,),\n",
       "  ion_channels={},\n",
       "  C=1. * ufarad / cmeter2,\n",
       "  V_th=-20 * mvolt,\n",
       "  V_initializer=Uniform(low=-70 * mvolt, high=-60.0 * mvolt),\n",
       "  spk_fun=ReluGrad(alpha=0.3, width=1.0),\n",
       "  solver=<function ind_exp_euler_step at 0x000001D9B8073E20>,\n",
       "  na=SodiumFixed(\n",
       "    size=(1,),\n",
       "    name=None,\n",
       "    channels={\n",
       "      'INa': INa_TM1991(\n",
       "        size=(1,),\n",
       "        name=None,\n",
       "        phi=1.0,\n",
       "        g_max=100. * msiemens / cmeter2,\n",
       "        V_sh=-63. * mvolt,\n",
       "        p=DiffEqState(\n",
       "          value=ShapedArray(float32[1]),\n",
       "          _derivative=None,\n",
       "          _diffusion=None\n",
       "        ),\n",
       "        q=DiffEqState(\n",
       "          value=ShapedArray(float32[1]),\n",
       "          _derivative=None,\n",
       "          _diffusion=None\n",
       "        )\n",
       "      )\n",
       "    },\n",
       "    _external_currents={},\n",
       "    E=50. * mvolt,\n",
       "    C=0.0400811 * mmolar\n",
       "  ),\n",
       "  k=PotassiumFixed(\n",
       "    size=(1,),\n",
       "    name=None,\n",
       "    channels={\n",
       "      'IK': IK_TM1991(\n",
       "        size=(1,),\n",
       "        name=None,\n",
       "        g_max=30. * msiemens / cmeter2,\n",
       "        phi=1.0,\n",
       "        V_sh=-63. * mvolt,\n",
       "        p=DiffEqState(\n",
       "          value=ShapedArray(float32[1]),\n",
       "          _derivative=None,\n",
       "          _diffusion=None\n",
       "        )\n",
       "      )\n",
       "    },\n",
       "    _external_currents={},\n",
       "    E=-90 * mvolt,\n",
       "    C=0.0400811 * mmolar\n",
       "  ),\n",
       "  IL=IL(\n",
       "    size=(1,),\n",
       "    name=None,\n",
       "    E=-60. * mvolt,\n",
       "    g_max=5. * msiemens / cmeter2\n",
       "  ),\n",
       "  V=DiffEqState(\n",
       "    value=float32[1] * mvolt,\n",
       "    _derivative=None,\n",
       "    _diffusion=None\n",
       "  ),\n",
       "  spike=ShortTermState(\n",
       "    value=ShapedArray(float32[1])\n",
       "  )\n",
       ")"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# neuron\n",
    "neuron = HH(1, solver='ind_exp_euler')\n",
    "brainstate.nn.init_all_states(neuron)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5c89948cc11a871",
   "metadata": {},
   "source": [
    "`solver` 可以选择积分方法，不同的积分方法适用不同的步长且有不同的精度。\n",
    "\n",
    "同样，`braincell` 中也内置了丰富的积分方法，如 `exp_euler` 、 `rk3` 、 `rk4` 等。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "aa2b99e926f9ad53",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-16T06:38:01.797639Z",
     "start_time": "2025-09-16T06:38:01.794767Z"
    }
   },
   "outputs": [],
   "source": [
    "def step_run(t):\n",
    "    with brainstate.environ.context(t=t):\n",
    "        spikes = neuron.update(30 * u.uA / u.cm ** 2)\n",
    "    return neuron.V.value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "3dd6dbd65f011395",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-16T06:38:02.799105Z",
     "start_time": "2025-09-16T06:38:02.642466Z"
    }
   },
   "outputs": [],
   "source": [
    "# simulation\n",
    "with brainstate.environ.context(dt=0.1 * u.ms):\n",
    "    times = u.math.arange(0. * u.ms, 100 * u.ms, brainstate.environ.get_dt())\n",
    "    vs = brainstate.transform.for_loop(step_run, times)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "4f4f61195f37e81c",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-16T06:38:03.756398Z",
     "start_time": "2025-09-16T06:38:03.676646Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj4AAAGxCAYAAABiPLw8AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAV71JREFUeJztvQmcXFW1t73S6TFJp7szdRIyAxJCEsAgEBBEwAQ+rl4EuQ5cLiA/EAyDDAKRy+CAQfCCMijDK5GraBCVi0RBIvMQDCRACIEAQkjIPHZn7LG+39pVp+qc7q5O1Tm7qvrUeZ7vq7eqT3WafZd7+O+11l67VywWiwkAAABABCgpdAMAAAAA8gXCBwAAACIDwgcAAAAiA8IHAAAAIgPCBwAAACIDwgcAAAAiA8IHAAAAIgPCBwAAACJDaaEb0NNob2+X1atXS3V1tfTq1avQzQEAAIAM0HrM27Ztk+HDh0tJSXq/DsKnAyp6Ro4cWehmAAAAgA9WrlwpI0aMKD7hc9NNN8nMmTPlkksukZ/97Gfm2e7du+Xyyy+XOXPmSFNTk0yfPl1+8YtfSH19fcZ/Vz09juH69++fs/YDAACAPRobG43jwlnHi0r4vPrqq3LPPffI5MmTPc8vvfRS+etf/yoPP/yw1NTUyIUXXiinnHKKvPTSSxn/bSe8paIH4QMAABAu9pSmErrk5u3bt8vpp58u9913n9TV1SWfNzQ0yK9+9Su59dZb5dhjj5UpU6bI7Nmz5eWXX5ZXXnmloG0GAACAnkHohM+MGTPkpJNOkuOPP97zfOHChdLS0uJ5Pn78eBk1apTMnz8/7d/TkJi6x9wvAAAAKE5CFerS3J1FixaZUFdH1q5dK+Xl5VJbW+t5rvk9+l06Zs2aJd///vdz0l4AAADoWYTG46PJxprI/OCDD0plZaW1v6sJ0homc1763wEAAIDiJDTCR0NZ69evl09/+tNSWlpqXs8995zcfvvt5rN6dpqbm2Xr1q2ef7du3ToZOnRo2r9bUVGRTGQmoRkAAKC4CU2o67jjjpO33nrL8+zss882eTxXXXWVOcJWVlYmTz31lJx66qnm+2XLlsmKFStk6tSpBWo1AAAA9CRCI3z0XP7EiRM9z/r27SsDBw5MPj/nnHPksssukwEDBhjPzUUXXWREz+GHH16gVgMAAEBPIjTCJxNuu+02U6ZaPT7uAoYAAAAASq+YXm4BSfQ4uxY/1ERn8n0AAACKa/0OTXIzAAAAQFAQPgAAABAZED4AAAAQGRA+EGp2t7QJaWrQk/ojAPRsED4QWjZub5L9r3tCzpzd+QoTgHyz+JOtMv7aJ2TW394pdFMApHF3i8xZsEK27mwudFN6HAgfCC1z31wt6ux5/r0NhW4KgPzkiXfN+z3Pf1jopgDIFX94U67+81tywW8XFbopPQ6ED4SWXr16FboJAElK6I/Qg3hy6TrzPv/DTYVuSo8D4QOhhXUGehK9S+iQtiBvD3IJwgcAwAK9UeJWaNjZIkfc9LTM/LP3bkbIDrpjehA+ABFne1Or/OG1lbKjqbXQTQk1JXh8rPCbV5bLmobd8vsFKwrdlFBTVsLyng4sA6GFZcYO97/4kVz5x8XyxTtfLHRTQg26xw6bdnAKyQalvemQ6UD4QFH4ctvayQnwyyuJ5McPN+zAjgEoZYdtha07WwrdhKKgrDf9MR1YBkKLez/T0tZewJaEm3GD+yY/Y0f/EOqyw2Y8PlZA+KQHy0BRwILtn14uCdnOaRrfEFmwAwX37FBGh0wLwgeK4tRCaxsLtg07EuryDx4fOzS1somxAR6f9GAZCC3uNbqlncnSRuE9zOgfChjagcKkdiC5OT0IHwgtba7wVgseHzseH0JdVur44DnzD44zO3CcPT1YBkJLq2txaSXHxzdurdOKy8dKqIucM/9QAdsOZaXYMR0IHwgtbi8PHh871wOge/zjTqkgT8U/hLrsQHmF9GAZCC1uLw87bP+4ozKEuuycjqM/+sft8GknZGjlVBehVy8IHyiSUBcD2y/uI+wsNHbs2IzHx0qSuHuMg38oU+EF4QOhxZ2PwqkuSx4fFhrfuG2Hx8eSx4cF2wqMay8IHwgtbi9PCzvsAKTsSKjLzuKCx8dOjg8eH/+4hzIC0gvCB4oj1MUE6Ru3s4ydoX/ctkNA2vH40B/tgBm9IHwgtJDcbAf3bpCFxj9u8U3k1c5xdvqjHbCjF4QPhJYWT04FA9sv5PjYwe3lIbRgB/qjf9yW49CCF4QPFIXHhwKGlur4sGD7ps0lvrGjfwi92oFxnR6EDxRHaIFx7RtCXfb7I3b0T4xkeytQnys9CB8oilNd7Gj847YcC7alekiY0c6CTQjbSn9kevSC8IHQ4l6kET7+IccnFx5I7Ogb7o6zf8qQce0B4QNF4RJnnbEU6sKQvnEnkJJMastzhh39gvBJD8IHQot7M8gE6R8uKbWD2zuBgPSP23LU5/KPuwvSHb0gfKAoPD7Mj5ZO0TBD+sa9q8aM/iHZ3g7uscy49oLwgdBCSfZcLDS4fPxCaMEO7qGMHf3jDrdiRy8IHwgtFOjKxamuAjYk5JBsbwdOGdrB7eVxh7MB4QNFU6CroE0pGjuy0PiHU112oD/agbvj0oPwgdBCqMsO7rUFO1ry+OA58w2hrlycMixoU3ocCB8ILe4pEVeuf0gmzcGCTX/0Df3RDtwdlx6ED4QWQl12oIBhLupKYUe/uE3HcXb/uPP1GNdeED4QWgjR2IGcilyEaArZknBDQU07UAgyPQgfKI5TXYxr3xCisQMLjX24q8s/nDJMD8IHisJTQWjB0oKNgrQkxLGjXwh15aKOT0Gb0uNA+EBRQIjGP4QWLEHo1Qpe22FHv5DcnB6EDxRJaKGgTSmeXCkM6RsKQdqBEHYuyitgSDcIHwgt1PGxgztMSGjB1ilD7OgXcqVyMT8WsiU9D4QPFMntw4xsv3Cc3Q54zixBf7QCIez0IHwgtHA7ux3wVNiB/mgHdx+kO/qHuwzTg/CB0EIdn1x4fArZknBDWQA7cDrODmxo0oPwgfBCDNsKTJB2IPRqBw4t2IEQdnoQPhBauCLADuSm2IeFxj8cWrAPdvSC8IHQwgRpB3bYdsCOdsBzFpyOdqM/hlT4zJo1Sz7zmc9IdXW1DBkyRE4++WRZtmyZ53d2794tM2bMkIEDB0q/fv3k1FNPlXXr1hWszZDPW5wL2pRQQ65UDoQ4K41vuHw4OB2HMR7IkAqf5557zoiaV155RebNmyctLS0ybdo02bFjR/J3Lr30Unnsscfk4YcfNr+/evVqOeWUUwrabsgdnvquLNi+4eqPXJzqwo5+Ibk5OB2thh29lEpIeOKJJzw///rXvzaen4ULF8rRRx8tDQ0N8qtf/Up+97vfybHHHmt+Z/bs2bL//vsbsXT44YcXqOWQKwh12YFCZ3bgVJcdCBkGp+N8yPwYUo9PR1ToKAMGDDDvKoDUC3T88ccnf2f8+PEyatQomT9/ftq/09TUJI2NjZ4XhANK29uBSrm58EAWsCEhhxyfXIS6CtWSnkkohU97e7t85zvfkSOPPFImTpxonq1du1bKy8ultrbW87v19fXmu+5yh2pqapKvkSNH5rz9YAeOYduBHbb9/khOhX84ZWg37KowPxaB8NFcnyVLlsicOXMC/62ZM2ca75HzWrlypZU2Qr53hoVsSbhhh20HQq+2QIgHpWP3Q0CGNMfH4cILL5S5c+fK888/LyNGjEg+Hzp0qDQ3N8vWrVs9Xh891aXfpaOiosK8IHyQTGoHQl124IoAO3DK0D7knIXU46M7URU9jzzyiDz99NMyduxYz/dTpkyRsrIyeeqpp5LP9Lj7ihUrZOrUqQVoMeQadth2IFfKDhzDtgMh7Bx4fDBjOD0+Gt7SE1uPPvqoqeXj5O1oXk5VVZV5P+ecc+Syyy4zCc/9+/eXiy66yIgeTnQVJ9wxZQc8PnZwW44dtn8Q4jk41YUhwyl8fvnLX5r3Y445xvNcj6yfddZZ5vNtt90mJSUlpnChntaaPn26/OIXvyhIeyH3UH/GDu0u0YgZ/eNeXBCQ/sGOwaGOT5EIn0wWtsrKSrnrrrvMC6IFA9s/nEayAzk+dqAsgP31knEd0hwfgI5QeM8OJJNagv5oBa7+CA4en+5B+EBoITfFDtTxyUGOD4b0DUniwaGAYfcgfCC04BK3A3ee2ReQ2NE/eCCD07H/dSxoGHUQPhBaOPZqB+xoB+7qsoN7kUZA+qOj2TCjF4QPhBaOveZih13IlhRTQc2CNiXUkLuXgxwfDOkB4QOhhSRIO5ArZQf6ox0oTBqcjp4yuqMXhA+EFkI0dnAv0pjRlgcSQ/oFz1lwONXVPQgfCC0sNHbAjrmoh1TQpoQat9ghx8cfHccxdvSC8IHQQi6AHbBjDm655xSNb/DkWoC7uroF4QOhhePDdiDHxw6UV7ADd/AFh1BX9yB8ILTgqbADAtIOeCpykZSLHf3A7ezdg/CBooAJ0tYOGzv6hfIKuag/gyH90DHUih29IHwgtFDa3hJ4zqzAMexchGgK1JCi8/hgSDcIHwgtnP6wA6GuHFwRgB19QajLDh3thoD0gvCBonDnEqLxDyGaHOywScr1Rcf+R3/0Bx6f7kH4QGghtGAHknKDwykaO5CbkhswoxeED4QWPBV2wI7B4YoAO+CpsAN27B6EDxTFYsPO0FLhPezoi45CBzv6g2PYufGcIXy8IHwgtFDHJzgkk9qBhSZXSbnY0Q/kSnUPwgdCC3dMBYekXDvgqbBDR7MxrP3BKcPuQfhAaKGOT3BIys0N2NEfeHxyNK7Z0HhA+EBooY5PLnaGBWtKkVUcLlRLwg2esxz1Ry7N9YDwgaJYtKnjYysXADv6AU+FJTqFXrGjPzhl2B0IHwgtHMMOTsedYBsLti8IGdoBAWkHNjTdg/CB8EKoKzCEFuxAHR87ICDtQOi1exA+EFrckyITpB0QkLZOI2FHP3DHlB0or9A9CB8ILYS6gkNowQ54znLlqcCQfqA/dg/CB0ILd3UFhzo+dqAQpB0IGdqBKyu6B+EDReHOZVz7g5wKO7DDtgP90Q4d7YbnzAvCB0KLe3FhgvQHdXzsQI6PHRCQuQFPrheED4QXhE9gOPZqB0JddsBTYQdCXd2D8IGiCHWxo/EJE6QVuCLADoS6cnWqq2BN6ZEgfCC0kNxsf4LEjP7gdJwdOlZqRkD6g9Nx3YPwgdBCHZ/g4BK3BAXjcgL90R94zroH4QOhhTo+waFgnB1YaHKV41OwpoQaxnX3IHwgtLjHNq5cf7Bg2wHPmR2wox2wY/cgfCCUUOgsVwUMMaQf8FTkxo5cmusX+mN3IHwglLCjsQOnP+yA5yxXdixQQ0IO82P3IHygSI4PM7D9wARpBzyQduA0kh0Q4t2D8IFQwkJjByrl2gEBaQcKQeaoLABm9IDwgVDCjiZXdXywow0woz8oBGkHrlDpHoQPhBIKxtkBT4UduGrBDoxrO+DJ7R6ED4QSBrYdqPdhB/pjrnJ8CtWSYju0gCHdIHygKGCH7Q88PnYg9GoHPD6WQIh3C8IHQgk77NzAOuMPku3tgBC3Q8f+x8bQC8IHQgmuXDuw0NiBZNLcgBn9wfzYPQgfKJIdDYuNHwgt2IFj2HagP+aqInuhWtIzQfhAKOlK5DBHZg+Vcu1A6NUO2NEO5Jx1D8IHQklXw5h7fSx4KlhpfMFCYwc8PrkZ15jRC8IHQklXA5lJMntYsO3AMWw7cBWNHRjXERQ+d911l4wZM0YqKyvlsMMOkwULFhS6SWAZQl124DSSHUgmtQP90Q7knEVM+Dz00ENy2WWXyfXXXy+LFi2SAw88UKZPny7r168vdNPAInh87NCVyUgSz56OyaP0RX9wytAOeCAjJnxuvfVWOffcc+Xss8+WCRMmyN133y19+vSR+++/v9BNA4t0NY7ZHWaPY7JevVLPsKMNj0/BmhJqOtoNO/oDARkh4dPc3CwLFy6U448/PvmspKTE/Dx//vwu/01TU5M0NjZ6XtDz6corweDOHsdkvV3KBzva2GFjQztJudjRSl2pArWjp1JUwmfjxo3S1tYm9fX1nuf689q1a7v8N7NmzZKamprka+TIkXlqLQShq51gjFoVWeOInN4lCB8bOPoRT4U/SMq1gyMYU/0ROxat8PHDzJkzpaGhIflauXJloZsEPkILCoM7exyTlbqED2YM7jmjL/rDsZvTHRGQ/mjv2B/ZFHoolSJi0KBB0rt3b1m3bp3nuf48dOjQLv9NRUWFeUHISAxsnSCdQU4dH/8CsgSPj50FO9EhnUrivdzJU7BnkkK8RJrb2umLvunYH7Fj0Xp8ysvLZcqUKfLUU08ln7W3t5ufp06dWtC2gV2cYVzSqxfuXBueCo/wKVx7wopjMneuFN0xe5y+V5JYmbChLQ9kYdvT0ygqj4+iR9nPPPNMOeSQQ+TQQw+Vn/3sZ7Jjxw5zyguKB0fk6LjuLb2kNRbfZUN2kNxsh1iaXKkSwePjxwOpHh8RPD5+cazmhLCxY5ELn69+9auyYcMGue6660xC80EHHSRPPPFEp4RnCDfOOO6l/58Z2zEGd4CFxr1gkyQewOOD58yOxwcvbiAcszkhbPpikQsf5cILLzQvKF6S47iX+f8NDG47oS5ypWyFDLGjX89Zae94rIsxbee0Jjk+ERA+UPw4d/jouFavj/sZ+MuVcmDBzh5nYXHbETMG8FQk7MiCbWdcM6aLOLkZoocObGeTzdgOlitFkriNUFfqGXYMkuMT74xtbGYC5pzFf8aMXhA+EO4Fm12NnVypXu5ddmHbFEZIEs9NyJAFOxjUleoahA+EPAmS4+zBSIVoSCi1sMPuTXKzlcJ7noKaGNK3gEz0R0zoBeEDoS/Jnjq5wOj2fzpObcku224dHwxppyxAARsU9uRmPD5dgvCBUOJMhrpYU6TLjh2THh8MGaxyc/JZARtUJMfZ488wZPDj7NjQDcIHQorrVBe7muCeM89JmgI3KoyQ42MJdwHDONgxuAcSEe4F4QNF5KkoaJNCiTMfmgrYCEgrdiTnzMaVFZQFsBkyJOzqBeEDIa/3wakuO6e6SBK3VX8Gz5l/HJs5x9kV+mPwSuJ4fLwgfCDk9Weo42Mt1MUk6Rv34szpOLtXqNAffdCpLABGdIPwgdDX8SHHx06Ihmq5dirlcjrOP47N8PjY2xgq+iPjOgXCB8IfWkhWJ2VgBwvRxD+zYAergM3pOAtXf3BprqXkZtczumMShA8UTcVhhI//0IJ7d8g1AT4gx8cKVMC2A5fmdg/CB0K9YOsiw5FNO8nN5KZYqOODEM9Bjg92DDI/OjA/pkD4QMiPs7uODzOyrdx5xjqTPcmux+m4QDglKdzjug07+j8d57lCBTs6IHygCHIqEiEaBrb/pNwSPBVWclM8Hp8CNyrkSeIIcRv9kXpIXYHwgdAn5aaKdBW2TaEkeVcXngpbl+amyitgRytJ4tgxcB0fBTumQPhA6OvPcJw9eC6AN0m8wI0KIXh8LOESkJQF8A9J4t2D8IEicInHP3MaKUBOhavwHp6KADlneM7sCHF3IUjGtZWyAJgxBcIHQklyMtQ7pgh1+SZpMldOBRNkMM8ZHkg7hxbI8bF3Sal5hiGTIHwg9DkVLDSWQjSUt7ec41PQJhVBeQX6Y+AK2J5TXYVrT08D4QMhr1ORcokT6sqe5Cls7piyFFrAU2HvKhrvM/B/O7uCHVMgfCDcO0OhgKGVJHGOD1vrj3gq7B5nZ1xnj9P3dG5EQHYG4QNFc2UFMewgCza5KbaOYbPQ2BDiJNvbyZViQ9MVCB8I/e3DVHj1D6fjcpHjg6fCZn0u7Jg91EPqHoQPhHzBpm5KEJJzIadoLN7VFX+GpyJ7koszp+MsCciUHTFjCoQPhHyhcVduZmQHu6vL+wyyoMvTSIVtUtg9PvRHO1dWYMfOIHwg9LkA5FT4x3W3JpVyLXl86I+2Thn28hTZhMwhx6d7ED5QNPU+2pggs4adYe4WGuwYNLkZO9oJvWLHjiB8IJS4Cje7kiAZ2H7hdJzNqxbYYdsIdeE5s1Xg1fsMED5QFJdCep+BnxwfclPsLzQYMljOGf3R7qW5GNIB4QNFc2UFoa5g9ZBYsO1WbmbBDpJz5r76A0MGSwVwnmFHB4QPhD4XIFW5mYGdLW6TsWAHvzTXLDSJWZX+GLT+DP3RL9gxB8Ln+eefl9bW1k7P9Zl+B5DXnWGiF7OjCVjAEDtaPY2EHYPWn4l/RkBmD5c450D4fP7zn5fNmzd3et7Q0GC+A8jvZYZOqIuBbWdniB2tLDSEXn3DnWf2cyDpjwGFjxrVGdxuNm3aJH379vXzJwH8l7bHlesf1+k46qZYXmhYsAOEDDkdZ6vAKwKyM6WSBaeccop5V9Fz1llnSUVFRfK7trY2Wbx4sRxxxBHZ/EkAX3AXjeVj2DpBJj5z51nQisMIcRshbEJdtupKxT9jRp/Cp6amJrm7qa6ulqqqquR35eXlcvjhh8u5556bzZ8ECAQ7Gns5FTHj9yE3JegdU5yi8Q9JubmoJM78GEj4zJ4927yPGTNGrrjiCsJa0DMmSG5x9k3KZhQ6s59MWtg2hV2IU5jUkgeSU4bBhI/D9ddf7+efAVjDyUNxu3IZ2EFCXfE8HwU72rmdHTsGKFNhkpu9zyDo1R8FblQYhc/BBx/cZUJzVyxatChImwCyOIbtTsplZGeLK0LDBGkBTiPZG9cUJrV9SSn9MWvhc/LJJ2f6qwD5LW1PqMs31PGxgyO63ZWbMWOQEDaeXHs5Ps6zwrYplMKH8Bb0KDynaOKfmSCDucRT9WewY5AdNqeR7FYSR4gHyznDA2kpx8dh4cKF8s4775jPBxxwgAmHAeR/Z8iCbeOuLkJdwXOlvJdCFrhRYfdAYkffUFcqB8Jn/fr18rWvfU2effZZqa2tNc+2bt1qqjbPmTNHBg8e7OfPAvjYYbPQ2Ct05n0GPvqjKykXOwY7rYnnzO7GEDMGrNx80UUXybZt2+Ttt982V1foa8mSJdLY2CgXX3yxnz8JYGGHzci24RLHjEF32IRofOM6zs6Gxs5xdur4WPL4PPHEE/KPf/xD9t9//+SzCRMmyF133SXTpk3z8ycBAu+wqTgcICmXHba1HTZ1fCx5Kki2t7ChcYe6Ctqk8Ht82tvbpaysrNNzfabfAeQcZ4ddkip0xvwYINRVwp1n9q6siH9GQNrKOcOOdur4YMdAwufYY4+VSy65RFavXp18tmrVKrn00kvluOOO8/MnAQKcoiG52S+c/rADOWf2Pbnccm83dw/PWUDhc+edd5p8Hr26Yu+99zavsWPHmmd33HGHnz8J4L+OD6EuOxWHndL2rNjBcnwI0VipJI7nzPbGsLBtCn2Oz8iRI011Zs3zeffdd80zzfc5/vjjbbcPYI+hBUJd/nFEjtqQ3BQ7x7DxQNq5q4tke/9whUoOPD4rV640g/sLX/iCOeGlr1yKnuXLl8s555xjvEp6I7x6mLSgYnNzs+f3Fi9eLEcddZRUVlYacXbzzTfnrE3Qk469EqKxU9reeYYds8UROZ66UpjRyl1d9MegOWf0RyvCR0Ncn/vc5+S+++6TLVu2SK5Rr5ImTd9zzz3mCP1tt90md999t3zve99L/o6G2fRE2ejRo01hxVtuuUVuuOEGuffee3PePigc7hh2GyM7a5zwoCY2cwzbVo6P8ww7BjmNhAfSbrkPxnVA4fPaa6/JoYceKj/4wQ9k2LBh5h6vP/7xj9LU1CS54IQTTpDZs2cbYTNu3Dj50pe+JFdccYX8+c9/Tv7Ogw8+aDxA999/v6kirQUWtabQrbfempM2Qc/J8eE0kt36M9gxaGiBEI3dU4YYMlucfB7vFSoFbVL4hY9eTaEelRUrVsjjjz9uKjWfd955Ul9fL9/85jclHzQ0NMiAAQOSP8+fP1+OPvpoKS8vTz6bPn26LFu2rFuvlIo19Ra5XxCmY69a74MdjZ36M95nkD3xHJ/4Z+wY8DSSk2yPHbOGCtg5ED4OOlnqNRUa8tJEZ83BeeCBByTXfPDBB+b02Le+9a3ks7Vr1xrh5cb5Wb9Lx6xZs6Smpib50twgCJtLPP6ZUFf2tCV2hprcjMfH1ilD7OiXtq6S7TFk1lCmIofC55NPPjEJxAcddJAJffXr189Ub86Uq6++OnncLt3LOTXmrhekoa/TTjtNzj33XAnKzJkzjffIeWniNoRrR0Ooy+5lhnjOgt6N5H0GQRfswrapeK5QKXCjwn6cXZOMf/e738lLL70k48ePl9NPP10effRRk1icDZdffrmcddZZ3f6O5vQ4aMFE9TAdccQRnZKWhw4dKuvWrfM8c37W79JRUVFhXhBOvHdMMbL97rAdGyos2HZO0dAfg12hgoC0G8J2Ep7Bp/D50Y9+JF//+tfl9ttvlwMPPND3f1xzgzK9yV09PSp6pkyZYhKdS5wAcIKpU6fKNddcIy0tLcnrNObNmyf77bef1NXV+W4j9PTjw65QFxOk/x12Sa/k4s0O207oFTsGOGXoCr0yrIPVlXLs6IS1wafw0aRmTQrWujlz587tdD+XnrqyiYqeY445xniUfvrTn8qGDRuS3znenG984xvy/e9/39T7ueqqq8xt8T//+c/N0Xco3oFN3RR7p5Ec87HDDnY3EnWl7IS6SMq1I8SdAq/YMaDwefLJJ+WMM86QjRs3dvpOB31bW5vYRD03mtCsrxEjRnQ54WhisrZrxowZxis0aNAgue6668xpMyjuBZuB7R/HZponFUvm+BS2TcWT41PYNoU/1MWGJniOj6s/Yshgyc0XXnihSS5es2aN8fa4X7ZFj6J5QPo/ZFcvN5MnT5YXXnhBdu/ebRKv1fMDxX+ZIQPbVi6A4xLHjjZyfFCQAYS4CXV5n4G/wx8ISEseH00avuyyyzodHwfIGy6PD6EFO6EFJ/kROwYLLbDQ+McR3e4QNkniwQoYJnN8sGMwj89XvvIVefbZZ/38UwDrd0ylQl2FbVO4LyllwbaX4xN/xkITzOPDlRV2UwEQkAE9PnfeeacJdWlYadKkSclTVA56VQRALnG8E+67kRjYAXNTEjbFjsEqDifrSrFiZ41jsvjdcc4z7Bgk9EqBV0vC5/e//71JJNZb0NXz4yhzRT8jfCC/pz/ITbFhR+dsJgtNsFOGJNvbCHXhgbR1SSkecUvCR+vl6NFxrbzcsZ4OQD7gklL7oS5px47Bk+1TOWfUTbGT3IwHMlgqQFJAMrCT+FItegv6V7/6VUQPFA5X4T0uM7R0KSShBUs5Fd5n4FP4JDokntyg4xoPZEd8KZczzzxTHnroIT//FCCHl0IysG3cjYQZAwpx+qPVU13onuzxnjKMfybZPmCoS2v16OWkf//7303tnI7JzbfeequfPwuQMbEuXbmFbVMYcSZDnRwTkS4W7MAFDBE+NpObCXUFO2WYOtVV4EaFXfi89dZbcvDBB5vPejWEG3eiM0DuY9h4fKxUeFVPBXd1WfVAkuNjq3IzHTJbKEyaA+HzzDPP+PlnADnJqSA3xT+Olywe6qKAoY2QoZPjg6ciwLimjo+1cU3OWWfITobQ3zHlJEEyQQYJdaXsyIIdpG5KfNFW2GFnTxt1fKzgWMzjOaM/JkH4QChxFpXeJSW4xK1cZug+ho0dg10KyRUBgUNdrkrimNFOf2RYBwx1AfQc4eMKdTGyfdtRvRS9yfHxTVJ0azIpC7ZvKK+Qu0tKEeIp8PhA+D0+hLosHWePfybUFSS0wBUBdsa1O8cHO9rIOcOOKRA+UAQeHybI4AXjuCLAVuVmrqywc8cU/dFyAUMMmQThA+FPymWHbWWCdDwVLNjZ4ywqpb25q8v2uMYDGaTOWSrZnukxBcIHQkmrs9Bobgo5FYGPvVIpNxitCUN6QjTU8QlWx8dZsLGjlXIfbAxTIHwg5JdrkgtgryxA/Bk77AChV82pIJnUyl1deCDtVLZPbQyxowPCB4rgOHviGQM74M4QARnUA+m+VZycigChLtedZ4zroCFs7NgRhA+EPrmZu2hsXP1BiMZWjk8qp4IOGajiMCFsKyHDVM5ZgRvVg0D4QChxdi/q8SHUZSe0QN0UGx6fVEFNpwox+Au9EuqykQOZ8ojjgUyB8IGiKWBI8p6dSyFZZ7KnzZ1sT66UhYKalFew0h/xQHYJwgdCn+NDqMtuAUMmSP87bM+VFazYWUN/zEXOmdMfC9yoHgTCB8I9sD3HsJkg7dyGjR0D5fhgRzuh1+RxduwYyAPJqa5OIHygCI6zx5+xww4QWvDc6VPgRoV8h51MJmWHHeiUIULcVl2p+DNOdaVA+EDxLDSM66xx5kJTf4bcFCs7bEI0NoR4L2PL+LMCNyrU/TGVCsD8mALhA6HEWVTiCw2u3OC3OBPqsrHD9uT4YMdAntxkIUhcZ1ZyfAgZpkD4QChpbXMXOos/Y6EJcjeS6xQN60ygUzSpUBf9MVBys+PxwYxZof3OmQrNxpBTXZ1A+ED4r1pgRxM81EWIxmKoi9CCjcrNTuiVce3P26P0Nsn28c/kQKZA+EDRuHJZr7PHmQy9l5RiyEAFDJO5KdgxW5xwdTznLL48YUd/YdeOp7owYwqED4Q+F4CcCjunulKVcgvbprBfUornzE5/TOX4YEffHh+PBxI7OiB8IORXVqRuFWdg+98dlvUuMSdAFBaaAB4fzfFhofHt7Unm+LhCXWxosqPNlRRV5vJA0h9TIHwg1MnN3h1NgRsVYjuq8EkuNBjSfwFDTzJpgRsVMtzrMhWwg4twNZ/n8Ad2TILwgaI5zk4SZPa0JIqkuCsOM0Fm76ng+HBw3J4d9Zppn1TwVPhPtFe40qczCB8ogruR4s+YIP3bUV3iTqgLO2aHW9+4LyklRJMdbsHtvqTU8UpC9lWbFac+FxuaFAgfCP/dSIQWfOMsKnE7Jp5hSN+naLjzLLj3MRV6xY5BqzYr5Jx1BuEDRXMbtkJ4ITtaEos2hff8495Je44PUwgyK9yeHSN88FT4wh12VfCIdwbhA+FPJk3pHgZ3FuiC4pgrHuqiLIAf3Aszx4eDe3zMUXbXHXz0x2A5PnjEO4PwgVDiSSZ1KR8Gt7/Qgie5mZyKAB4fPT7c+TnsmeZEf1Rvj4IHMviJV4VDC51B+EDoS9t7Ql3sDgPnVLDDzg53TpTjrVBYZ/yXVlCcDQ05Z0FPdcWfc4lzCoQPFMHdSKnnCB9/ORXx00jsDANVbU4kNhPqCibEyxLH2FO5UtjRT96eFtP0nOqiPyZB+EDo0J1LqrR9R49PARsW0gmyV8ecCozoW/goCB+7oS4W7GCnulJ1pQrarB4FwgdCXzeFUFfA0EJJidkVJk/RYMNgyaRUyrUS6koJ8YI2K/Q5Phxn7wzCB0Kdm6LuXE+oi8XGVw0fd06Fzo/kA/g/PkxSrqVQV1L4oHxsVG5GiKdA+EBoXeJKuSspV2FsZ29HZ4J03hUmycxxFubOoa6CNit0tCSFuDdEQ1/0V1DT2dAgfDqD8IHQ0dyaEj4VpfEwjQPuXH83syvusgCcpMmc1jR1UwgZ+vX4lHjsSVf0m3NGrlQ6ED4QWuGj3h5H9CSrkzJL+g51ObkACgLSX39UOI0UTPiU46mwKsQdO3LnWQqED4R3oSlNdV9qpwS4mb3DzlBhscmcpkR/rCjrbd7ZYVsKdSF8rCQ3pzxn2NEB4QOhzU1xCx8uhgxwM3uHHbbCYpO9ENewq9uO2hXx+gSv44OAzI7mtjZPf0zeck9fTILwgdCHFhSOEAfw+HQI0SjYMXOaWts8QtwjIFm0fef4EOryR1OLI8R7e0LZiPAUCB8IbWjBE+rC4+M/tOBKynW0Dwu2n4WGkGEu6vgoLNqZs7sl4fEpw+NTNMKnqalJDjroIBPaeOONNzzfLV68WI466iiprKyUkSNHys0331ywdkL+dtgKu8Psae0iZJgML2DH7HN8nB02p+MCVm7unGyPHf30xw6n47BheIXPlVdeKcOHD+/0vLGxUaZNmyajR4+WhQsXyi233CI33HCD3HvvvQVpJ+Q31OWEa1iw/Xt8FBJKgyfbuyuJY0f/oVfnlnsFT65/IZ481UVfDKfwefzxx+XJJ5+Un/70p52+e/DBB6W5uVnuv/9+OeCAA+RrX/uaXHzxxXLrrbcWpK2Q31NduHP977CdCdK7OyxYs0Lrgey4w1YQPtmHupJlAbCjlf5IXakQC59169bJueeeK7/5zW+kT58+nb6fP3++HH300VJeXp58Nn36dFm2bJls2bKl29CZeovcLwjLgu3y+OCpyJrdzfEJsqo8JXw4SRM8tODJlaI/+q4kTpJ4sJyzykR5BebGkAofvTforLPOkvPPP18OOeSQLn9n7dq1Ul9f73nm/KzfpWPWrFlSU1OTfGluEIS3jg+DO3N2JZIgqxITpDfUhcvHRrI9/TFzmhL90VmwPacMKb6XMbs7enxcfZE7+HqA8Ln66qtNknJ3r3fffVfuuOMO2bZtm8ycOdN6G/RvNjQ0JF8rV660/t+A3NZNUSga51/4OAuNd3dYsGaFPqdCoT9mz86EB7JPwgOJxyfgKcOyzqFXdHicUikgl19+ufHkdMe4cePk6aefNqGsiooKz3fq/Tn99NPlgQcekKFDh5pwmBvnZ/0uHfo3O/5dCF8BQ9y52bMrGepy5Uphx8A5FR7hg6ciY3Y6HsiE8IlvfikEGVSIu+/g03Hd2/VzVCmo8Bk8eLB57Ynbb79dfvSjHyV/Xr16tcnfeeihh+Swww4zz6ZOnSrXXHONtLS0SFlZmXk2b9482W+//aSuri6H/1dAoXY0ngKG3Efju96HO9RFiMZy6BVPRfZCvIMHUk8fYkc7yfacjusBwidTRo0a5fm5X79+5n3vvfeWESNGmM/f+MY35Pvf/76cc845ctVVV8mSJUvk5z//udx2220FaTPkY2B3cRqJgR0ox4cF23Koi1yprIWPE+pK5afEEOJZsLtDcrPbw8Op1xAJn0zQxGQ96j5jxgyZMmWKDBo0SK677jo577zzCt00sMz2pvgE2bci1X05zu5/oal0n+oi1BU4p0IhVypIqCs1rumPwT0+lAUoEuEzZsyYLrPTJ0+eLC+88EJB2gT5Y0dTq3nvV+Hy+HAfjV2PD3bMuj92LcRRPpmyq7m1k8eH0Kt/j09S+FBQM5zH2QG6FD6VnXeGeHyC5fhw2Wv2bE/0x2qX8KEQpP9TXZ66Us6GhtCrj41hvD9SV6ozCB8IHdu62GGndoasNFl7fFwLTWningAWmuz7o7PQeJLt6Y/BPJDJcV2wZoWObbsTQrwyfshHwXPmBeEDod/ReEM0BWtWaHfY7jo+HGfPnu27Wzp5IEm2t5TcjIDMmm2J/ljdhUecQwtxED4Q3pyKLpIgmSAzp2FXfILs794ZJmYEhE/2oa4uPT6UV8h6XHeV48OwzgwdtzsSArJL4UN/NCB8ILQLjSfUxQ47axoTwqe2j1v4cMt9tmxPhhY6e3ywY2aonZyQYU1V6r5FPBX+5sZ0OZDYMQ7CB0LHjsRxdvcO21lo2GFnhp6K3LozLnxqqty5APF3Jsjsd9gejw+XvWYtwh1TefojAtJXmEuLaVJXKj0IHwhviKaKeh9B8nucE3AsNHZ22G4PpFNegVOGmbE1MaZVPHL5cPDE5v4ub49CXSkvCB8IFTubW5OnPwb2S92xhivXn3jUCdGTU8FCkxVbdjSbd7WhO0ncCRlSVyoztu5s7iTCFcor+Atfu090KdSV8oLwgVCxaXtzsjhXXxZs3zhhLs3v0csgHciVyo6N25vM+yCXCHeHDPH4ZNcf6/p6F2zKK2THxsT8OLhDf6SulBeED4SKTYkdti403gWbpNxs2JywY22fVCKpZ2dIrlSWwqe86wWb/pgRWxIen1pXYrNCeYXs2LBtt3kfXO0VPpQF8ILwgVCxKbHQDOy00DBBZsPaxvgEObR/pec5IcPs2JDYYXf0+CR0Dx6fDFnXGB/XQzos2JRXyI4NaYQ4nlwvCB8IFRu2JYRP3zSeCibIjFjnCJ+ayjQuceyYCRsT/XFQhwWbEE12rG3Y1WV/pOKwv/mxo8cnWecMT64B4QOhYuWWneZ9RF0fz3M8PtmxtqFrjw/HsLNj1db4gj2sox1ZaLJidaI/Dqut8jzHA5kdaxJ2HFKdRkBiRwPCB0LFys3xhWbkAO8ESS6AvwW70w4bO2bFis1xIT5qYBohzkKTEWsauhaQyRAN/TEjPt4U74+jO/RHxrUXhA+E0uMzEo9PID7csN28jxvc1/OcCTI7ViQWmlED+nTtOcOOGRXT/GjDji4XbELYmdPS1p7c0IweyLjuDoQPhGqC/GBdfMEey4Ltm6bWNlm5JT5B7jO4n+c77JhdlVwnSXxMh4UGIZ5deEarX6vN0i3Y5Epl5n3U/lZZVtIpSZxLc70gfCA0fLJll7nPp7x3ieydZsFmZ7hn3l+33UyQWt21UxIknoqMeWfNNvM+rKZS6jok2yMgM2fZ2rgdxwzq66narGDHzFmyqsG8jx/aPxn6dyDnzAvCB0LDm59sNe/71veTMuecawJ2NJnzxsq4HQ8cWeuphaSQTJo5byUWmgOG9+/0HQt25ry+Yot5nzyiptN32DFz3vok3h8n7dXZjsyPXhA+EBpe+XCTeT907IBO37GjyZwFH2027wePrO30HcmkmTP/X/H+eMiYzv2RBTtz/pnoj58eVdfpO8cDyYK9Z+Yn5scpozvbkVwpLwgfCE1+z9PvrDefj9h7UKfvUzkVVCbtjubWdnnuvQ3m89GfGtzp+1SFVybIPV1O+tIHG83nI7voj4ReM6NhZ4u89nHc43P0vvTHIOUp3l7daD4fsc/ATt87l+YixOMgfCA0uxmt9aGXQR61b/qFhhBN9/z97bXmglJNfjyoC48PBQwz429vrTGX5eqpuIl7dRHqwlOREY8tXm0W4/FDqzuVBFDoj5nxf2+sSnp7OtbwUThl6AXhA6Hgnuc+NO+nfHovzy3YDiTlZuY1+/XLy83nrx86Sko75EkpFDDcM7oI/+/8uB2/MmVEpzwppXdih03otfvj1/e98GHSjl1Bfa49s7ulTX714kfm81cPGdnl73DK0AvCB3o8z7y73oRndOyee9S4Ln/HWWgY2On5y5urZeHHW8zJGRU+XZHKTclz40LEwwtXypJVjdKvolT+Y08LDQIyLQ+8vNwU3NN7pdL2x6QQz3PjQsTdz/3LXFWxV22VnHzwXl3+DjlnXhA+0KP5aOMO+e4f3zSfzzpibKc6Hx0nSHIBuuadNY3y348sMZ8v+vw+nSo2O5Ar1T2LP9kqN/xlqfl88XH7dLqctHNoATumO6jwkyfeNZ8v/cKnpG9F6R4WbOzYFc8sWy93PP2B+XzlCft1KgfgQCqAl657G0APQJNHL5nzumzc3iz7D+tvBnY62NGk5x9L18mlD71haiDpibhvfW7vtL+bCi3ksYEh4fG31sh3/7jY5PZoYvg3jxyb9ndTAjKPDQxJuPVPi1bJf//fW9LSFpMTJw6Vb6Tx9ih4INPb8aFXV8p1f3nbzHmaAvDvB3Xt7VGYH70gfKDH8e7aRvnFM/8yoRlFRc9vzjm0y9weB2LYnXlz5Vb5xbMfyN/fXmd+PmR0ndx3xiFpd4UKSbmdF5jFnzTI7U+9L0+9Gz9VeNjYAXLnNw7uMkfKAU9FZzu+unyLseOLidNwn99vsNz21YO6zJFyoD92tqPa746nPpAFy+NlAFQ8/uTUyd3+u94l8b7K/BgH4QM9YjC/v367PP3uenliydpkgT2d8zT2f92/TehW9CjsaOJ2XLqm0Rz7//vStSYPxbHNWUeMkatOGN+t6PEeH47ugq19SKvgvvD+Bpm7eI28m6gsrHa84HN7y0XH7SMVpZn2R4ksrW3tpsjj8+9tNKeONGytlPXuJd85/lNy/uf2TtopHSQ3p+w4b+k6Mz9+mLCjXk1x6fGfMnmPHSs1dySRAhlpO7pB+EDeB7HeKfOvDTtM3olWbX3zkwbZvKPZ4735woR6ufDYfeSA4Z2rkHZF1HY0rYkLCfX6Cc07URvq+5adLcnf0as9Tpo8TC44Zm/5VH11Rn83aiEarWukN4NrX1ShuGR1gxHeWzvY8d8mD5Nvf34f2WeI96qUdETN46PjbvmmHbJ0daMR31pT5vWPt5jwqoOWovjSgcNlxuf3kZEdLnVNh+NUi0runvZHteMH63VcN8iiFVtMRWYNrzpoUr2egjv36HEmoTkTojY/7gmED1j1OOhEt75xt6xvbJL12/S1W9Y1NsnqrbvMYNZBrbH9jqgnYuq4gXLs+CFy4qShXdaiyGSCLKaBbS7BbNht6hep/ZZv3GF2e3qzuorHruxYVdZbjtxnkBy3/xAjHtMl36ajmCo3a39UQb16q9pwl7Fh/LXbiEYVPNpHu4qiVFeUyuF7D5TP7zdETpo0TGr6lGX13y6mZFK1Y+Ou1pQNE/1xTcKW+lz7aVfipKaqzIxr7YsnTByaNok5HaWJBbsY+qNjS62jpf1v1Za4PfXzRxt1Mxgf113NYdWVpSavbNqEejNHVldm2x+jJSD3BMIHMh6sukisc4ka/bzBJW70fXfLnne46qLVS0b3HdLPFNE7aFSd7D+seo/hg0x2NGEZ2DuaWs3CqzdTr3EtHrqo6IKin9275a5QsThuUF/jFTtwZI1MHlEb2I5hKm2/q7nN2M3YL7GAmAW5IfVzU+ue+6PaUfvixOE1phjhxL1qzH1H3eXw7Ikw1ZXSOjDxfhi3mWM/t8DR29MzGdd6QeaE4f1lwrD+cuCIWvN5T+Gs7ghbXSkd12sbd5vxqy8d46s69M+de7ClenT2HtJP9h9aba7xOHhUrZkv9xTOisqGxgYInwijg2DLzuaEZ8YlZBpTQsb5Tl2wmaK75SH9K4zXRt/r+1ea196D+5pQwfCaqkCDuCt6yiV8zu54TeMu1+SXmAgTE6Iu1tt2dy9qHPQG9eG1Veb4+ZiBfU2l4LGD4q9itqMKho3bm5KLheOpiS/I8c/u8Gh3aJVqteHw2kpjs/jnKhMmGFZbKQP7lnebYOuHnlJXSv/7OqYdD5fbjo7A2ZShHQf0LTc30Rv7Jd6HuT7rGA8icrrC0Z6FXrD1v6920jnSjGedI533xLN1GWxWHLR2kfY/py+OGtDHzI360v5qvT/2kKs/WtvaZcP2Jtm4rVkmdXEpbb5A+BQh7kHa0SOj3pp125pkg363vanLcEk61G2tg1InOH0frKImIW5U5NQn3qvK/Xsc/JKvS0p1IVG7qkt6ZeL1yZb4IuJMgO54fHeo+1oXkmE1Van32krPs2xDA0HJVzKpToBqK7XdJ1viNnQ+6yKdLnTSkb7lvWWvOrWVI2YSC7NZkKukvqYikAfML/mqK9XU2mY8Xk7oxNgy+XmX6ZOZtEFDpEYYJuw2zPVZn6t9Czquc2xH3dipCHT3x5WJd+2LOndmOleqx0Y3KkP7V5p3FTiOyIn31co9HtawTb5Ox+1uaTM2c/dH57O+a3/U/yl1g/Xej060vnHLFIRPiNDFQgWNES8acnKJmg0J74w+17o32SxcupMzQsYlapICJyFm9Lt8D9ZssF0pVxcUzUl6b902c6rnvbXbZPkmnQi7zq3pSF2fMhmaEC9qR33XSdC868+1VWaC7GkkT39YtOM7a7bJ++u2GXvqS0/w6SS4pz6qu1S1VXJBTnpq4guxflaPmO3dsQ1shxZ0QdEEbO2L/1q/3eSDaL6XCu89/Sc62lFt57ahPtdNTU+0o+0yFTubW40d9VCA9kOnT6oXcU9dXs2jOXOOoHHenfFdn/i5R45ry6kAO5tbzWGAZetc/XHDDjOu99yWXsZWmj5R17dcCkHP+18ooujAdnYajvfACY84noRN25v2OMm5B6m68J1wk8dT43hn+lfK4H4VezziHAaCXsKnovKfH22W59/fIK8tj5+kaE5ztEkn4xF1VeZkir70s+6MHWGjdu7JIrE7eidiC20+PWe6c9YLZZ9/b4O5HkNP+aSzo/a7EYld8Ii6uB31ZXbIdVWm79oOneSLoKEFFTpa3fjZZXE7am2rdIJbc2vUZmrDvVw2NP2yNuR2DOip0DwwHdNaDFVPSKkITzdHqB3d/XBkwp6O91U3f2UB8r4KSdCQ4famVmNDvTpo0cdbzIYw3Z9yvLCmPyb6of7seL4G9qsoeH9E+BQIHXz/eGedPPXOOnl9xVZzZ026BcJNSWLXYXJnEqImKWT054TAGdivPLSDNMjOMNuFRmuL/OrFD+Wvi9d4joIrugveb2i17FdfLZ8aWm1ylDQWrxNhoQdurkjdjZSdHXVhvv/Fj+Txt9Z2ynNQAT5+WLXsMziew7DPkGqTq6Siu1Cu7p7qqdCikw/+82NTP6hjEqza8YC9aowd9x7SV8YN6mf6pC7IPdFbU8jQq5bJuP+l5TJv6dpOBy50jtRxve+Qatm3Pt4nNX9O826K1Y5+PD6xWEzm/2uTPLhghcx7e12n9Uk9Xpq8rvbTfqgJ2OMG9zPe7p5uR4RPAVBPznm/ec3UaXBTUVpiPAhuV2q941I1CcIVPUIt90ScZNJMdzSaGHvjX9+RP7/+SdLFrQP22PH1cvi4AfKZMQNk9MA+PX4AF3pnqHkRP3hsqTzx9trkM12Ij99/iLkeQ0+lqFiMmh2zXbDV26t3gOlmyEG9DMfsN1iO2HuQOf2oO+eo2THb03Fa6uH6v7wtL7wfrw6tqJdB++MhYwbIlNF1xgsWNbIt97Fs7Ta57tElxgvuoPOhlnc4fNxA0x/T3fcXBhA+eUZd2Gf86p8mvqz5CacdMlKO3GegKTCXi1M6USGbZFLdVZ/zwKsmF0rRuhhnHznG1BsJcoQ5ajtDdX2f/9uF5oSadtsTJw6TM48YY67GiHo/zsbjo3epXfaHN6QxYUct8nf64aONHaMmdIKcjvvbW2vMnXRawkCrQ+vdVf81dbQpTRB5Ozr1kDLw5P7h1ZVy7aNLjB01HP0fh4yQr31mlCnzUCwgfPLMnxZ9YkSPhqse+fYRGVcwhQyPYe9hglzw0WY5a/YCE0bQ2i03f2WyHDyqLk+tDJHHZw8TpC7WKnpUIB04slZ+cuokU8MFsqs/oxefXvj7183Crrvon5422YQCIbvQqy7WV/15sfHe6kZy1pcny6iBzK3Zbgx//dJHcsNjS81n9Tbe+OVJGVeHDhMInzzzx4WfmPfzPzcO0WORTO6Y0jomF/x2oRE9n91nkNx9xpQeeQKjpyeJf7B+m1z0+9fNJPrFA4ebxboQR8Z7MqUZeCo0L+rSP7xhfufUT4+Qm06dFKm8PFun415bvlm+98hbRvSccfhoueFLB5AOkKY/dmfHlz7YKD+YGxc9es3Nd6ftV7SeW2b9PLJ1Z7MJsyh6hxLkN7Rww1/eNuUADhjeX/7fmYeE9uRVIU8j6cR51Z/eMrWKdGd9238cGPnwoJ+7kYwd/7jYJN4ete8g43lkse7Oc5a+XMJVf1ps+qvep/aDfz8g8mEtPxXZd7e0yZV/XGxOamlo68rp+xW1HZmx8ojWPdCONWZg/GQQ2GNPdyPp0eC/v73O/N5tXz0I0eMzZPj3t9ea49V6ZPWWryB6/IYWHlu82lwsqx7H//mPAxE9e/Scde3JnbNgpbnwWE9kaVimmBfrXI7r+1/6yNTg0YT6739pYtHbkVkrj2ixJ4VciBwKnzRbQz1qrXz1MyMzvqk8inSXm6LHW+95/kPz+ewjx0bydIytEM3/eyHeH7919LisL+SNEt2FXvXZrxLj+uLj9jXlJyD7Ctgtbe3y65eWm8+XT9uvIBW68w3CJ49oboTyqfp+hW5KpDw+euTaOSb8zSPH5L1txbJga/G3N1ZuNSc99PQW+AsZarj7rVUNxo56egsyGNddOHy0MKFeHaOC5ytTRuS/cSGiu8r2T72z3lT9V6+ZniiMAgifPKK39Cpa0RJyNUHGugzP6ONDxwzgxEyABdup1XPMpwabWj3g77LXuYtXm/cTDhhqrosBf5Wbn0z0xy8eOEz6lJOumlFdqS484o+9Ge+Pp04ZURRV/DMhGv9X9hD09l5FixJC/io3P7Nsg3n/woT6vLerqATkkvhCc8LEoXlvV9jo7tJcLfuv0B/3TDohrh7JeUvXm8/TJtAf/Xp8Wtva5YX3N0TOjgifPKL3bilahRnyc4pGL9PTxGbl8+MHF6RtxbDD1vviNEdN58/jxrNg+/X4aNj1vXXbjR31NBf4C72+vbpRNm5vkuqKUlNJGPxVwH7zkwZTOFPDhQeOKJ4ChXsC4ZMn9Lig3karIHzyN7A1J0UvzhxeU2nukgF/O2y9T07RxPCaPiSR+j0+rHcfKVr0sbYPYS6/V3+8sXKLef/06LrIhGdy4cl9+YP41R5a1yxKJzSj839pgdmyszm5E+xfRTw6X6Xt317VaN4nj6gt+iOaudxhv55YaKhyHaxg3JLV8fv5tEoz+K/c/MbKhqSABP/C561VcTsePCpadkT45InGXfEbq/tX9fyba4vJ4+MsNJMi5Ma1ssOOde3xidoEadvj4wjxA4bTH4Ms2I7H56CR2DGI52zpmmj2R4RPnnDCXNSayG+IxtnRaLVmyCIJ0pWUq/V7liTsiKfCfyVx9f44C83EveiPfoXPtt0tpmihcuAI+qPf/tiws0U+2bLLfJ4wLFr9MVTC569//ascdthhUlVVJXV1dXLyySd7vl+xYoWcdNJJ0qdPHxkyZIh897vfldbWuKel0DQmhI/eyA65G9ju8MKOplb5aGN8giymm4XzXcBwXWOTud9MF6Gxg/oWsHXhXrA/3rxTtje1mpwU8s38X5rrjGktqTCwH2UV/I7rpQkRPqKuKnJ5e6FZhf/0pz/JueeeKz/+8Y/l2GOPNYJmyZIlye/b2tqM6Bk6dKi8/PLLsmbNGvmv//ovKSsrM/+mp3h8NNQF9nFfpqden/KSXvLxpp3m4sK6PmUyiAnSd8E4Z6EZWVfFJZoBPJAfrN9u3vcd0g87Bqjc7PTHsQMR4UHKfbyXuElg/4h5e0IjfFTkXHLJJXLLLbfIOeeck3w+YcKE5Ocnn3xSli5dKv/4xz+kvr5eDjroIPnhD38oV111ldxwww1SXl7YExSNuxE+efP4JHY1WtVVGcUEmX1ycxc7bLw9wey40umPAyhgmimlXZSpWL4xbscxg7BjkEMLKxL9cXQE+2Moth2LFi2SVatWSUlJiRx88MEybNgwOfHEEz0en/nz58ukSZOM6HGYPn26NDY2yttvvy09Jrm5EuGTC9yXPDq7mhWb4ws2C40PT4XL5fPRxrinYuwgwjNB7JgU4vTHjEnoHq/w2RQf12MQ4oE8kCsT/XFkBPtjKITPhx/GL0ZUz81///d/y9y5c02OzzHHHCObN282361du9YjehTnZ/0uHU1NTUYcuV+5YEdzXPhUk+OTc+HjJOZGeUcTvIChdPb4DGahyb6AYerZJ1uiu9DYPM5OqMuOx2dlIrF55IDoXTZcUOFz9dVXm6Pd3b3effddaW+P75quueYaOfXUU2XKlCkye/Zs8/3DDz8cqA2zZs2Smpqa5GvkyJGSC3Y1t5n3yrLiv/m2kBOke5JcsTk+sNlhB0vKXb4pEVoYiB2zP87e2eOD8Am2YDsen9EIH9+XOMdiMfnE6Y8RvDuyoO6Hyy+/XM4666xuf2fcuHEmUbljTk9FRYX5Tk9yKZrUvGDBAs+/XbduXfK7dMycOVMuu+yy5M/q8cmF+NnVEhc+VQifnCU361qj49pZbKLsyrUpfJw75obXRm9nGLyAoSQXmpUI8cAhGj2puXVnS2Q9FbbGdcOuFtnW1BrZS7MLKnwGDx5sXntCPTwqdJYtWyaf/exnzbOWlhZZvny5jB492vw8depUufHGG2X9+vXmKLsyb9486d+/v0cwdUT/rr5yTUr4hCK6GNrwQktbzCw2ukNMhRaYIP3uDHWhcSbIeq5aydoD6YjwjdubzRygj4fXYke/Hp/125rMe9/y3lJNvqRv4bMyIcL1tGtVefQ246FIOFHxcv7558v1119vvDEqdvSEl3LaaaeZ92nTphmBc8YZZ8jNN99s8no0H2jGjBl5ETZ7Ynci1BXFTpbf8ELMLDZ6RYiKIIUF29/xYfVS6OWkSr+KUvOCbE91iceOutBUlDIH+K0/49iRMR2ssv2qrbuSNXyiSGhmMhU6paWlRtjs2rXLFDJ8+umnTZKz0rt3b5P0fMEFFxjvT9++feXMM8+UH/zgB9IT2N1Kjk8+PD5NicG9YXt8Zzigbzk1U3yXBYgXL1SG9C/85iG0yfbaHxOeiiHV2DFIXSnH40N/DObx2bA92v0xNMJHCxH+9Kc/Na90qCfob3/7m/RESG7O7+B2FprBFC70XQhS7ejssIeyw/YvfGKu/hjRhSb46bhEqCvRH4dU0x+D5EptjHh/ZCucJ3a1xLcsJDfnWfhEdGDbWLB1sSG0YMnjk9hhI8T9CXGnHlKqP2LHIAJyQ6I/RrWiPcInT+x2kpvJ8ckZvZ0qr64ddlRduTZCXbo7JNRlN9SFEA9WV8rpjwjxYAJyY8T7I8Inz6EuPD65w0nlaW1joQmaTKoQ6rJUV8rt8aE/BspNWb8t3h+xY7CCmhvw+EA+cI6zk+OTn3t9nCRIJsgAoS48FdauUMGOdsorrMfjY6Wg5saIC3GET55I1vEh1JWXSZKFxh+u9drYUcsCKAP6FPaS37ChVeUdW6qAdEILUd1h2/L4OAs2dvRfUDPmnh8jasfQnOoKM9rRzjh8tBE//bmrKz/JzRHf0QRZsNWOakN9bd0Vr5Jb04dicX48kM1t7XGPD/0xcF0pzU9p3B0vpllHf/RdUHNHc5vsThy2GVQdzQ0Nq3CeFpNr/y199WiwL3w2sTMMNEm2SUL4JDw+dXh8ssZEXttEdja3ybbEgk1/9B8ydES4UlOF8PGT3KyOsw2u6td9yqMpAQh1QdHtalra2s1dNEotO8OsSaRKSePulmT1a4SP/5wzJzyjCaZ4fP0Ln03b4yK8urJUSilK6vu05vrEgYWBERbh9B4ouklSLzF0Ti/UVrFg+12wnYWmvLREKrljznd/dOyoIly9v5A5ZYncFLeAZDMTrDDp5h2OFze6dmQ2g6JL4HMmSHXl6qIN2eHMkZtcEyQLdgDhs8NZsBHhfkW44oRo8D4G8/hsTIzrmgjbkVUBii4RMrXDju7AtuOpSCzYeM0C2VFvZldqyUsJ5PFxhA/5PcHqc21KjGs8PgBFtKtJ7bCjO7BtVMAmtGAn5ywpIBHiWaOeRmdcOyfj8PgE8/hsQogjfKB4KOmww2aC9Ee5EzLchh1teHycnAoEZLAQtuPxwY7BksQ3E+pC+EARenwSO0Nqz/ijLJEX5VwPwEJjKbk5wjvsIJQlTnClhE90F2wbBTU3JkPY0e2PCB8owmRSTi3YWGicaz9YaIIJ8Y2J0GtdX+xoRfhEeMG2Oj/2ja4dET5QxDtsFho7wie6E6SN0KvTH0nKDZbgnMzxifCCbTX0WhXd+RHhA0W3w97eFK+Sy4LtD6cEQOr4MHYM0h+dYprkSgU70s6CbSfZfnMyxye64xrhA0V5ZFNhoQmW3OxQw0JjpT8ixP3RsRYXdgye4Bz1+RHhA0V3+sOBCTJYqMsBj48/6I/2j2Ir5JzZET61EQ69Inyg6OrPODBBWhI+JOVa8vhgRz8gxHNjx/4IH4Dw02GDzQ7bVmghwhOkTU8FC3bw6s2qJasrsaMfKlz37dVUlXXyAEUJhA8UrccnyjHsIJR32BlGOQnS1g5bxWRVWe+CtqcY7Bj1BTsIlaWp/lcX8TGN8IGiocN6Lf0rSwvVlKLZYfcp7y0VrgkTMqeq3LvQcNFr8FwpvI/+qXQJ77qIh68RPlCUO0MVPaUdlRBkbUe8Zv5xe3iwo53+SJ6Ufypdoa66iNuRlQGKhr4VKQ9P1Hc0Nq6sUCi6Z2eHjR1tCXHs6Be357Y24nZE+EDRoGEZB1zidnJ8qJLrn6pydti2k8Tx+Njx+AyIuB0RPlA09C1PeXyYIO2c6sKOlkJdCEgrHsioeyqCUEGOTxKEDxRlMikTpJ3kZjxndoQPAtI/ZW6PD1XErZzqqo34/IjwgaKhbwXJpLZzKgZEfGcYhMoOp7rAQo4PnjMrdXwGRHx+RPhA0dDHFeoimdROqKu+f2VB2xJm8PjYwT2WB/atKGhbisfjUy5RBuEDRUNfl/Bhh+2fatfpOISPfzjObodhtVXJz8Nr6Y9WjrP3jfb8iPCBoqGPK9S1V12fgrYlzOw3tH/yc31/dth+KXHlpgzqh/Dxy5DqVB/cyyWCwL8nd0DEhTjCB4qGMteVFZ8ZU1fQtoSZ/YdVJz8PqWaH7Rd35fBJe9UUtC3FUqZiUD+EuF/2HRIf11PHDZTBLjEZRajpD0XD+GHVcuiYAbJvfb/Ix7CDoJdA3vjlidK4q1WG1iB8/PKFCUPlsi98So7adxBVxAPw2X0HyYEja2XCsP4eLxpkx/83aag8e8UxMnJAn8hfn9IrFovFCt2InkRjY6PU1NRIQ0OD9O+fcvkDAABA+NdvtiEAAAAQGRA+AAAAEBkQPgAAABAZED4AAAAQGRA+AAAAEBkQPgAAABAZED4AAAAQGRA+AAAAEBkQPgAAABAZED4AAAAQGRA+AAAAEBkQPgAAABAZED4AAAAQGRA+AAAAEBlKC92AnkYsFktebw8AAADhwFm3nXU8HQifDmzbts28jxw5stBNAQAAAB/reE1NTdrve8X2JI0iRnt7u6xevVqqq6ulV69eVpWoiqmVK1dK//79rf1d8IKd8we2zg/YOT9g5/DbWeWMip7hw4dLSUn6TB48Ph1QY40YMSJnf1//h2ZQ5R7snD+wdX7AzvkBO4fbzt15ehxIbgYAAIDIgPABAACAyIDwyRMVFRVy/fXXm3fIHdg5f2Dr/ICd8wN2jo6dSW4GAACAyIDHBwAAACIDwgcAAAAiA8IHAAAAIgPCBwAAACIDwidP3HXXXTJmzBiprKyUww47TBYsWFDoJoWaWbNmyWc+8xlTYXvIkCFy8skny7Jlyzy/s3v3bpkxY4YMHDhQ+vXrJ6eeeqqsW7euYG0uBm666SZT0fw73/lO8hl2tsOqVavkP//zP40dq6qqZNKkSfLaa68lv9dzKNddd50MGzbMfH/88cfL+++/X9A2h422tja59tprZezYscaGe++9t/zwhz/03O2Enf3x/PPPyxe/+EVTNVnniP/7v//zfJ+JXTdv3iynn366KWxYW1sr55xzjmzfvl1sg/DJAw899JBcdtll5gjfokWL5MADD5Tp06fL+vXrC9200PLcc8+ZxfaVV16RefPmSUtLi0ybNk127NiR/J1LL71UHnvsMXn44YfN7+tVJKecckpB2x1mXn31Vbnnnntk8uTJnufYOThbtmyRI488UsrKyuTxxx+XpUuXyv/8z/9IXV1d8nduvvlmuf322+Xuu++Wf/7zn9K3b18zj6jwhMz4yU9+Ir/85S/lzjvvlHfeecf8rHa94447kr+Dnf2hc6+ubbrJ74pM7Kqi5+233zZz+ty5c42YOu+888Q6epwdcsuhhx4amzFjRvLntra22PDhw2OzZs0qaLuKifXr1+uWLfbcc8+Zn7du3RorKyuLPfzww8nfeeedd8zvzJ8/v4AtDSfbtm2L7bvvvrF58+bFPve5z8UuueQS8xw72+Gqq66Kffazn037fXt7e2zo0KGxW265JflMbV9RURH7/e9/n6dWhp+TTjop9s1vftPz7JRTTomdfvrp5jN2toOO/0ceeST5cyZ2Xbp0qfl3r776avJ3Hn/88VivXr1iq1atitkEj0+OaW5uloULFxq3nvs+MP15/vz5BW1bMdHQ0GDeBwwYYN7V5uoFctt9/PjxMmrUKOzuA/WunXTSSR57KtjZDn/5y1/kkEMOkdNOO82Ebg8++GC57777kt9/9NFHsnbtWo+d9U4iDZtj58w54ogj5KmnnpL33nvP/Pzmm2/Kiy++KCeeeKL5GTvnhkzsqu8a3tJx4KC/r+uleohswiWlOWbjxo0mrlxfX+95rj+/++67BWtXMdHe3m5yTjRUMHHiRPNMB1l5ebkZSB3trt9B5syZM8eEaDXU1RHsbIcPP/zQhGA0JP69733P2Priiy82tj3zzDOTtuxqHsHOmXP11Veb28FVnPfu3dvMzTfeeKMJsSjYOTdkYld9V9HvprS01Gxmbdse4QNF4Y1YsmSJ2bmBXVauXCmXXHKJiblrYj7kTrzrTvfHP/6x+Vk9PtqnNR9ChQ/Y4Q9/+IM8+OCD8rvf/U4OOOAAeeONN8ymSRNysXN0INSVYwYNGmR2Fh1PuejPQ4cOLVi7ioULL7zQJME988wzMmLEiORzta2GGbdu3er5feyeHRrK0iT8T3/602b3pS9NYNYkRf2sOzbsHBw96TJhwgTPs/33319WrFhhPju2ZB4Jxne/+13j9fna175mTs2dccYZJjlfT4kq2Dk3ZGJXfe944Ke1tdWc9LJte4RPjlFX9ZQpU0xc2b2705+nTp1a0LaFGc2fU9HzyCOPyNNPP22Op7pRm+sJGbfd9bi7LiTYPXOOO+44eeutt8zO2HmpZ0JDA85n7BwcDdN2LMegeSijR482n7V/6+TvtrOGbDT3ATtnzs6dO03OiBvdmOqcrGDn3JCJXfVdN1C62XLQuV3/t9FcIKtYTZWGLpkzZ47JXv/1r39tMtfPO++8WG1tbWzt2rWFblpoueCCC2I1NTWxZ599NrZmzZrka+fOncnfOf/882OjRo2KPf3007HXXnstNnXqVPOCYLhPdSnYOTgLFiyIlZaWxm688cbY+++/H3vwwQdjffr0if32t79N/s5NN91k5o1HH300tnjx4ti///u/x8aOHRvbtWtXQdseJs4888zYXnvtFZs7d27so48+iv35z3+ODRo0KHbllVcmfwc7+z/5+frrr5uXSotbb73VfP74448ztusJJ5wQO/jgg2P//Oc/Yy+++KI5Sfr1r389ZhuET5644447zOJQXl5ujre/8sorhW5SqNGB1dVr9uzZyd/RAfXtb387VldXZxaRL3/5y0YcgV3hg53t8Nhjj8UmTpxoNknjx4+P3XvvvZ7v9UjwtddeG6uvrze/c9xxx8WWLVtWsPaGkcbGRtN3dS6urKyMjRs3LnbNNdfEmpqakr+Dnf3xzDPPdDknq9jM1K6bNm0yQqdfv36x/v37x84++2wjqGzTS/8fuz4kAAAAgJ4JOT4AAAAQGRA+AAAAEBkQPgAAABAZED4AAAAQGRA+AAAAEBkQPgAAABAZED4AAAAQGRA+AAAAEBkQPgAAABAZED4AAAAQGRA+AFAUHHPMMXLRRRfJd77zHamrq5P6+nq57777ZMeOHXL22WdLdXW17LPPPvL444+b39+yZYu5ZX7w4MFSVVUl++67r8yePbvQ/2cAQI5B+ABA0fDAAw/IoEGDZMGCBUYEXXDBBXLaaafJEUccIYsWLZJp06bJGWecITt37pRrr71Wli5daoTQO++8I7/85S/NvwWA4oZLSgGgaDw+bW1t8sILL5if9XNNTY2ccsop8r//+7/m2dq1a2XYsGEyf/58+fGPf2yEzv3331/glgNAPsHjAwBFw+TJk5Ofe/fuLQMHDpRJkyYln2n4S1m/fr3xBs2ZM0cOOuggufLKK+Xll18uSJsBIL8gfACgaCgrK/P83KtXL88z/Vlpb2+XE088UT7++GO59NJLZfXq1XLcccfJFVdckfc2A0B+QfgAQGTRxOYzzzxTfvvb38rPfvYzuffeewvdJADIMaW5/g8AAPRErrvuOpkyZYoccMAB0tTUJHPnzpX999+/0M0CgByD8AGASFJeXi4zZ86U5cuXm+PsRx11lMn5AYDihlNdAAAAEBnI8QEAAIDIgPABAACAyIDwAQAAgMiA8AEAAIDIgPABAACAyIDwAQAAgMiA8AEAAIDIgPABAACAyIDwAQAAgMiA8AEAAIDIgPABAACAyIDwAQAAAIkK/z97wTZ6k5zOxAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# visualization\n",
    "plt.plot(times, u.math.squeeze(vs))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14ee745b190f115d",
   "metadata": {},
   "source": [
    "观察运行结果，很明显，我们建立的 HH 型神经元产生了 spike ，精确地表达了它的电生理特性。\n",
    "\n",
    "当然，除了对单房室神经元进行建模，我们的 `braincell` 框架也支持对具有一定空间结构的多房室神经元进行建模，这里不做进一步讨论。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c8bb77e401f6640",
   "metadata": {},
   "source": [
    "本文中举例的模型来自：\n",
    "\n",
    "- Brette, R., Rudolph, M., Carnevale, T., Hines, M., Beeman, D., Bower, J. M., et al. (2007), Simulation of networks of spiking neurons: a review of tools and strategies., J. Comput. Neurosci., 23, 3, 349–98"
   ]
  }
 ],
 "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
}