Chemical Reaction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Markovian stochastic dynamics\n",
    "\n",
    "|State|Discrete|Continuous\n",
    "|---|---|---|\n",
    "|Trajectory view|Recursion equation(e.g. Ranodm walk)|Stochastic differential (e.g. Ornstein-Uhlenbeck process)\n",
    "|Distribution view|Master equation(e.g. random walk, birth-death, chemical master equation)| Fokker-Planck equation(e.g. diffusion equation)|\n",
    "\n",
    "Taken Brownian motion as an example:\n",
    "\n",
    "Trajectory view: $x(t+\\Delta t) = x(t) + \\sqrt{2D\\Delta t}\\xi(t)$\n",
    "\n",
    "Distribution view: $\\frac{\\partial P}{\\partial t} =D \\frac{\\partial^2 P}{\\partial x^2} $\n",
    "\n",
    "Solving both, we can arrive at the final distribution:\n",
    "\n",
    "$$P(x,t)=\\frac{1}{\\sqrt{4\\pi Dt}}e^{-\\frac{x^2}{4\\pi Dt}}$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Chemical master equation\n",
    "\n",
    "Chemical master equation is a stochastic model for chemical reactions. It is stochastic with fluctuations, many chemical species beyond a single \"walker\".\n",
    "\n",
    "One specis., one reactin: $x\\to^k \\emptyset$, which means that the species $x$ decays to nothing with rate $k$.\n",
    "\n",
    "Continuous: Reaction Rate Equation\n",
    "\n",
    "$$\n",
    "\\frac{dx}{dt} = -kx + \\text{noise}\n",
    "$$\n",
    "\n",
    "which is the rate of change of the number of molecules.\n",
    "\n",
    "Discrete: Chemical Master Equation\n",
    "\n",
    "$$\n",
    "\\partial_t P_t(n) = -k(n+1)P_t(n+1) - knP_t(n)\n",
    "$$\n",
    "\n",
    "where $P_t(n)$ is the probability distribution, $t$ is the time, $n$ is the number of molecules, and $k$ is the reaction rate. This function is intuitively the probability of having $n$ molecules at time $t$.\n",
    "\n",
    "The discre description is more accurate, but it is hard to solve for large $n$.\n",
    "\n",
    "#### A solvable example: birth-death process\n",
    "\n",
    "Define\n",
    "\n",
    "$$\n",
    "X\\in [0,N], \\emptyset \\to^{k_2} X, X \\to^{k_1} \\emptyset\n",
    "$$\n",
    "\n",
    "We have the equation\n",
    "\n",
    "$$\n",
    "\\frac{dP_t(n)}{dt} = B(n+1)P_t(n+1) + F(n-1)P_t(n-1) - (B(n) + F(n))P_t(n)\\\\\n",
    "F(n) = k_2\\\\\n",
    "B(n) = k_1 n\n",
    "$$\n",
    "\n",
    "where $P_t(n)$ is the probability distribution, $t$ is the time, $n$ is the number of molecules, $B(n)$ is the birth rate, $F(n)$ is the death rate, $k_1$ is the birth rate constant, and $k_2$ is the death rate constant.\n",
    "\n",
    "> How to solve it?\n",
    "\n",
    "Steady state: $\\frac{dP_t(n)}{dt} = 0$, or\n",
    "\n",
    "$$\n",
    "B(n+1)P_t(n+1) + F(n-1)P_t(n-1) = (B(n) + F(n))P_t(n)\n",
    "$$\n",
    "\n",
    "We can use a generating function to solve it:\n",
    "\n",
    "- Define $G(z) = \\sum_{n=0}^N z^n P_t(n)$\n",
    "- Then we have $G'(z) = \\sum_{n=0}^N n z^{n-1} P_t(n)$\n",
    "- And $G''(z) = \\sum_{n=0}^N n(n-1) z^{n-2} P_t(n)$\n",
    "- So we have $G''(z) = B(z) G(z) + F(z) G(z)$\n",
    "- Then we have $G''(z) = (B(z) + F(z)) G(z)$\n",
    "- Then we have $G''(z) = (k_1 z + k_2) G(z)$\n",
    "- Then we have $G''(z) - k_1 z G(z) - k_2 G(z) = 0$\n",
    "- Then we have $G(z) = A e^{r_1 z} + B e^{r_2 z}$\n",
    "\n",
    "Solve the equation ,we have\n",
    "\n",
    "$$\n",
    "A=\\frac{1}{r_2 - r_1} \\frac{r_2}{r_2 - k_1} \\\\\n",
    "B=\\frac{1}{r_2 - r_1} \\frac{-r_1}{r_2 - k_1}\\\\\n",
    "r_1 = \\frac{k_1 - \\sqrt{k_1^2 + 4k_2}}{2}\\\\\n",
    "r_2 = \\frac{k_1 + \\sqrt{k_1^2 + 4k_2}}{2}\n",
    "$$\n",
    "\n",
    "\n",
    "> This is only for stationary state, how about the time evolution?\n",
    ">\n",
    "> Omitted.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Gillespie algorithm\n",
    "\n",
    "Gillespie algorithm is a stochastic simulation algorithm for chemical master equation. It is a Monte Carlo method that simulates the time evolution of a chemical system.\n",
    "\n",
    "A pseudo-code of the algorithm:\n",
    "\n",
    "```\n",
    "# Define the system with the initial number of molecules and the reaction rate constants\n",
    "# Initialize the time to zero\n",
    "# Loop until the end time or condition is reached\n",
    "  # Calculate the total propensity function by summing all the reaction propensities\n",
    "  # Generate two random numbers from a uniform distribution between 0 and 1\n",
    "  # Use one random number to determine the time interval until the next reaction event\n",
    "  # Use another random number to determine which reaction will occur next\n",
    "  # Update the system by increasing the time by the time interval and changing the number of molecules according to the chosen reaction\n",
    "# End loop\n",
    "# Plot or output the results\n",
    "```\n",
    "\n",
    "This algorithm is based on the following theory:\n",
    "\n",
    "- The time interval until the next reaction event is an exponential random variable with the rate parameter equal to the total propensity function: $t \\sim \\text{Exp}(\\alpha_0)$\n",
    "- The probability of each reaction event is proportional to its propensity function: $\\text{Pr}(R_i) = \\frac{\\alpha_i}{\\alpha_0}$\n",
    "- The number of molecules after each reaction event is a multinomial random variable with the number of trials equal to the number of molecules before the reaction event and the probability of each outcome equal to the probability of each reaction event: $X \\sim \\text{Multinomial}(n, \\text{Pr}(R_i))$\n",
    "- The number of molecules after each reaction event is a Poisson random variable with the rate parameter equal to the propensity function: $X \\sim \\text{Poisson}(\\alpha_i)$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABuP0lEQVR4nO2dd3gb15W33wOw9ypKJCWRqlYvllXc5B73lsSxHZckLmmb9mWzabspu0k2m96d2HHiLttxbMdx4hZ3q9iSLNnqXRSL2HsngPv9cWcIEARAsIBgue/z4AEwGMxcAIMzZ373FFFKYTAYDIbJgyPaAzAYDAbD6GIMv8FgMEwyjOE3GAyGSYYx/AaDwTDJMIbfYDAYJhnG8BsMBsMkwxj+EIjI70Xkv6zH54hImc9rx0Xkggjue4aItIqIMwLbViIyZ6S3a237ORG5NcTr94nI9yKxb2v73xORWhGpjNQ+BouInCUiByK07e+IyEMjvM2PichbI7nNEPvq/Y8NYxt9/puGgYmJ9gCiiYhcD3wJWAy0AceA+4G7lOZT0RqbUuoEkBKt/Q8VpdQl9mMR+Rhwu1LqzNHYt4jMAL4MzFRKVY/GPsNBKfUmMD/a4xiLRPM/NhREpAhtJ2KVUq4oD2fITFqPX0S+DPwS+DEwFcgDPgWcAcRFcWjjEtFE+3iaAdQNxeiLyKR2ggyTi2j/UaOCiKQD/w18Rin1hFKqxfLwdyilPqqU6rLWC0uWEBGHiHxNRI6ISJ2IPC4iWdZrCSLykLW8UUS2ikie9dprIvK/IvKOiDSLyN983ldkSTIx9phF5F4ROSki5ZakEVAGEpHVIrLZ2t9JEfmNiAQ8mYlItoj83dr/Vmu7b/m8frq1vMm6P93ntddE5PsishFoB2ZZy24XkQXA74F1lmTV6LPbTBH5h4i0iMjbIjLbZ5tKRD4jIoes1/9HRGaLyCZrjI8H+iyW7PYSkG/t7z5r+ZUissf6Ll6zxmW/57iIfFVE3gfaAhl/EfmliJRa+94uImcF+h6tdS8Vkb3WuMtF5N+t5YFkwq+IyPsi0mb9rnmiZbIWEfmXiGRa69rHwZ0iUmH9nv8eYgxrre+qUUTeE5FzQqw7XUSeFJEa6/j8jd/rPxGRBhE5JiK+V3JBj0XRMtFGEfm5NYaj1jH0Met7rBYfKVD8/mMicpWI7LS+7yMicrG1/OMiss/6fo6KyCeDfS6/zyDWWKqtbe4SkcUicpqIVInPf0hErhWR96zHq0Vkm/WeKhH5mbXaG9Z9o3WcrbPW/4Q1vgYReUFEZvpsd0jHdERRSk26G3Ax4AJiBljvPuB71uNzgDKf144DF1iPvwBsAQqBeOAPwAbrtU8CfweSACdwKpBmvfYaUI6WmpKBvwIPWa8VAcoeI/CUtd1kYArwDvDJIOM+FViLlvKKgH3AF31eV8Ac6/Gj1i0JWAiUAm9Zr2UBDcDN1rZusJ5n+4z/BLDIej3WWna79frH7G35fad1wGrrPQ8Dj/qN7W9AmrXdLuBlYBaQDuwFbg3yuf1/o3loCe9Ca2z/ARwG4nx+w53AdCAxyDZvArKtsX4ZqAQSgqx7EjjLepwJrAxx7GxBX2UWANXAu8AKIAF4Bfi233GwwfrtlwA1eI+97+A9Zgqs7/ZStFN3ofU8N8BYncB7wM+t7SYAZ/r8bj3AHdZ6nwYqABnoWLTe6wI+br33e+hj5Lfo/8ZFQAuQEuA/thpossbtsD7PKdZrlwGzAQHWox2NgN+v3+f8ALAdyLDeuwCYZr22F7jEZ92ngC9bjzcDN1uPU4C1gf6X1rKr0MfVAvRx8p/AppE4piNmA0dzZ2Plhv4zV/ot2wQ0Ah3A2QEOyj4HF30N/z7gfJ/Xpll/nBjgE9a2lwYYx2vAD32eLwS6rT9M7wGGNhBd+BgntBF+NczP+0XgKb8DcY61nx5gvs9r38Nr+G8G3vHb1mbgYz7j/+8An2kgw/9Hn+eXAvv9xnaGz/PtwFd9nv8U+EWQz+n/G/0X8LjPcwf6RHuOz2/4iUEeOw3AsiCvnUCf6NMGGNdx4KM+z/+Knleyn38OeNp6bB8Hp/i8/iPgXuvxd/Aa/q8CD/rt+wUCGBVgHfoE0s/5sX63wz7Pk6wx2JJo0GPReu8hn9eWWO/N81lWBywP8B/7A/DzMH+Hp4EvBPp+/dY7DziIdoQcfq99FXjYepyFPpnYJ4U3gO8COX7vsX8PX8P/HHCb33HWjp5rgmEc05G6TUqpB33g5fhe2iulTldKZVivDfZ7mQk8ZV3aNqJPBG70n+RB9J/vUetS/UciEuvz3lKfxyVozzQnwPZjgZM++/gD2tvqh4jME5FnRaRSRJqBHwTYJkAu+sTiOwbfx/nWmHwpQXtigdYPF9+Im3b6T2JX+TzuCPA83EnvPuNXSnnQ4w17/CLy79YlfJP1vacT+LsE+CD6RFYiIq/bMkAQBvsZ/Y+T/ADbnAl82D5GrPGeiXZE/JkOlKjgE5S9v5FSqt16mEJ4x6L/Z0EpFc5vOB04EmgwInKJiGwRkXprn5cS/HfoRSn1CvAb9BVHtYjcLSJp1ssPAVeISDJwHfCmUuqk9dpt6CvG/aIlzstD7GYm8Euf76MefXXhe5yN1DE9IkxWw78Z7bVcNULbK0VfMmb43BKUUuVKqR6l1HeVUguB04HLgVt83jvd5/EMtAdeG2D7XWjvw95+mlJqUZDx3AXsB+YqpdKAb6APRH9q0JflhUHGU4E+qH2ZgfaabVSQMQz02mjQZ/wiIujPF9b4Rev5/4E2CpmWY9BE4O8SpdRWpdRVaCP4NPD48IbfB//jpCLAOqVoj9/3OExWSv0wyLozZPCT2oM9Fge77dn+C0UkHn1V9BP0lUMG8E+C/A7+KKV+pZQ6FX1FPQ/4irW8HG0LrkVf3T7o855DSqkb0L/l/wFPWCeIQMdLKVrq8v3eE5VSm8L72KPPpDT8SqlG9GXc70TkQyKSKnqCdjlatxwsvwe+b0/oiEiuiFxlPT5XRJZYk0jNaMPu8XnvTSKyUESS0BPOTyil3H7jPQm8CPxURNKssc4WkfVBxpNq7atVRE5Ba7SBvgc38CTwHRFJstb1PSn9E5gnIjeKSIyIfAT953k2zO+lCigc9YkrL48Dl4nI+dZV1pfRRivcP2Qq+sRYA8SIyLfQOm0/RCRORD4qIulKqR709+8JtO4Q+S/rN1qE1s8fC7CO7cF+QEScogMLzhGRwgDrvoOek/ihiCRb654x0CCGcCwOhnuBj1u/l0NECqxjMg49P1ADuERPNF8UzgatSdw11u/fBnTS93d5AH1yX4L+L9jvu0lEcq2rxEZrsccagwetz9v8Hvi69dvYk98fHuRnH1UmpeEHUEr9CPh/6B+9yrr9Aa37DfZM/UvgGeBFEWlBT9ytsV6bCjyBNgT7gNfx8Sysx/dhTRoCnw+yj1vQf4C9aJ35CQJfwgP8O3AjehLtHgIbCZt/Q8sXldZYNqCNI0qpOvQVypfREth/AJcrpfyvSILxCrAHqBSRcN8zYiilDqDnc36Nvoq6ArhCKdUd5iZeAJ5Ha8QlaKMRShq6GThuyWufAj46xKEH4nX0BOLLwE+UUi/6r6CUKkVfxX4DbaBK0d5tv/+5ddK/Aj3XcwIoAz4S5lgGcyyGjVLqHfRJ7efoK6vX0Tp5C/p/8bi1vxvR/7dwSEP/BxrQv2EdOoTb5iksqdZH0gIdALJHRFrR/+/rlVId1jrfBzZa0s5apdRT6KuCR63ffjdwCWMYe5beEAVE5DX0xNwfoz0WGxH5P2CqUurWaI/FMHEShsYyInIELdX8K9pjGS0mrcdv0IjIKSKyVDSr0ZNaT0V7XAbDaCAiH0Tr9q9EeyyjiclWNKSi5Z18tNz1U3TMscEwobGuuBei4/VHcj5mzGOkHoPBYJhkGKnHYDAYJhnjQurJyclRRUVF0R6GwWAwjCu2b99eq5TK9V8+Lgx/UVER27Zti/YwDAaDYVwhIv6Z94CRegwGg2HSYQy/wWAwTDKM4TcYDIZJhjH8BoPBMMkwht9gMBgmGREz/CLyJ9Htznb7LMsSkZdEtyB7Saz2cgaDwWAYPSLp8d+HrnDny9eAl5VSc9FVBr8Wwf0bDAaDIQARM/xKqTfQnWh8uQq433p8P3B1pPZvMBgMo4rHDe8+CC2VA68bZUZb48/zaW1WiW5NGBARuVN0l/ttNTU1ozM6g8FgGCr/+g4882/w50ugqXzA1aNJ1CZ3la4OF7RCnFLqbqXUKqXUqtzcfhnHBoPBMHbYuQE2/QrmXwatNXD/5WPa+I+24a8SkWkA1n31KO/fYDAYRpbSd+Dvn4eis+C6++Hmp8a88R9tw/8MYHd2uhVT991gMIxnmsrg0Y9CWj5c9wA4Y2H6aWPe+EcynHMDuoP9fBEpE5HbgB8CF4rIIeAC67nBYDCMP7rb4dEboacDbngMkrK8r41x4x+x6pxKqRuCvHR+pPZpMBgMo4JS8PSn4eT7cONjMOWU/uvYxv/Ba7Txv/VZSC8Y/bEGwGTuGgwGw2B5/Uew92m48Lsw7wPB1xujnr8x/AaDwTAY9v4NXvsBLL0eTv/8wOuPQeNvDL/BYDCEy8n34alPQeFpcMUvQSS8940x428Mv8FgMIRDazVsuAESM+EjD0NswuDe72/8mysiM84wMIbfYDAYBsLVBY/dBO11cP0jkBq06EBofI3/fZdFzfgbw28wGAyhUAqe/X9Q+jZc/TvIXz687Y0B428Mv8FgMIRiy+9g50Ow/quw+NqR2eb00+DmJ6Nm/I3hNxgMhmAcegle/E9YcCWsH+Eq8tNXR834G8NvMBgMgag5CE98AqYsgmt+D44ImMsoGX9j+A0Gg8Gf9nrY8BGIiYcbHoG45MjtKwrG3xh+g8Fg8MXtgic+Do2l8JGHIGNG5Pc5ysbfGH6DwWDw5YVvwNHX4IpfwIy1o7ffUTT+xvAbDAaDzfb74J0/wLp/gxU3jf7+R8n4T2jD/6Pn9/OfT++K9jAM4dBUBk9/FjqbBvc+Vxc883loOB56PbdLx2KffC/w6zUH4ZnP6e2NF965B/Y9O/r73fEwbL139PcbaSp2wj++DHMugAv/O3rjGAXjP6EN/6HqVrYdb4j2MAzhULJJx0pv/NXg3le+Hd69X4fdhWL/s7DtXtjzdODXdzwA7z4Ax94c3P6jydt/gB0Pju4+XV3w4jdh829Gd7+jwbHXweOCa/4ADmd0x2Ib/7hk3cR9hJnQhj8lPoa2ble0h2EIB1envt9yl66JEi62B99WG3q9d+7W98GuDGyDf/D58PcdbTw9g79CGi6HXoSOBj3xGQGDFFXa68AZD0nZ0R6JZvpquPMNyJg+4pue0IY/Od5JW9cEOzgnKrbE0tMGb/4s/PfZhr89hOGv3AUlGwGBhmP9X+9ohMr39eODL+gU/fGA2zX6hn/nBn3v6YGWk6O770jTXqe7aIVbcXM0iETuABPd8MfF0NZlPP5xgW34F16tJZnG0vDeF47H//YfICZRp9sH8vhPbAblgSXXQdMJqN43mJFHD88oG/62Wjj0AuQt1s8bSkZv36NBe8PY8fYjzMQ2/PExdLk8uNyeaA/FMBC21HP+t/T96/838Hu626Fmv37cXhd4nfZ62PUXWHod5K/QMkVHY991jr2pL/HP/YZ+Pl7kntGWenb/VZ9s1n9VP2+caIa/rm/f3AnMhDf8gJF7xgOuLkAgaxasug12PgK1h0K/p3qv9tRjEoN7/O/er08qaz4JmUV6mb/Xf/wNradmFcO0ZVruGQ+4XdDdqu9Hg52PwNSlVqtBgcYTo7Pf0aK9DhKN4R/3JMfpmXkzwTsOcHVCTILWV8/6sn786g9Cv+fkTn1fdEZgjd/t0mGHRWdB3iIfw++j87fXQ+VuKD5bP593CZS9A21BriDGEh7ruO5qjvy+qvfp73vZDbqMQeq0iSf1dNQbqWci4PX4jeEf87i6tEEBSMmFdZ+BPU/qVnfBOPm+9tCmLdfemsdP0jv4HDSVwuo79fNAHn/JJkBB0Zn6+bwP6KuIw/8a/meKNJ4efd/ZGPl9vbcBHDGw5MP6ecaMiSX1eDxaBjRSz/gnxTL8rcbwj31sj99m3b9BQga88r3g7zn5npZmknO0se7wy9l4+w+QPh3mX6qfx6dCUg7U+3j8x9/UUlHBqfr5tOWQkjc+dH63bfgjrPN73PD+4zDnQn1SBsicObE8/s5GfQwZj3/8k2RLPUbjH/v4evwAiRlw5hd1FMmJLQHW79Ya/7Sl2phDX7mnao826qfdBs4Y7/Ks4r4e/7E3YcYa774dDph7ERx+2WtYxyIeN2CFnUba8B99TYduLrveuyxjJjSX699hItBer++N4R//9Eo9RuMf+/h7/KAlmpQ8+Nd3+8fW1+wHd7fl8Vt/Vt8J3nfu1ttbeWvf92UWeQ1/Wx1U7/HKPDbzLoaupsAnnLGCx+eYjrThf28DJKTD/Eu8yzJnAkpLaRMBOyrMTO6Of1KMxj9+8Pf4Qaern/0VOLEJjrzc9zU7fn/acki25Afb42+vh/ce03q0v2abWazrArl7oOQtvazo7L7rzDoHnHFjW+7xvRrxD08dSTqbdT2gxR/s+/tkzNT3E0Xn77A9fmP4xz1J8bbUYwz/mCeQxw/aY8+YAS//d1+v/+R7EJeqDbkt9dge/46HwNWhQzj9ySwC5dae6rE3ITYJClb2XSc+RUcCjeWwTo+P4Y+kx7/3b/q7XHZj3+WZtuGfICGdtsdvDP/4xzu5azT+MU8gjx8gJg7O+YY29Pue8S4/+Z7W9x0Ory7bXqe17633wMwzYOqS/tvLKtb39cfg+Fu63roztv968y6GukNQd2T4ny0S+NbJiaThf+9RyJ4Dhav6Lk/N11E+E2WC12j8E4fEWCci0G40/rGPO4jhB511mzNfR/h43PpWuUvr+6BPDvHp2uM/+Lz2Qu0QTn/skM6ybVCzT3v2gZh3kb4fq16/exQ8/objWg5bdn3/+jXOGEgrmDhST3udlvfiUqI9klFhQht+ESE5LsaEc44HXF2BpR7QJXLP+0+oPQjvP6Yzel0dXsMPeoK3rUaHcKYVwCmXB95WylRdnmHnQ/p58dmB18ssgtwFY1fnHw2p573H9P3S6wO/PpFCOu2s3bFUoC2CTGjDD3aFTmP4xzyuzuAeP8CCK/RE7qv/qzNroa/hT8qBsq26prp/CKcvDoc26o0ntHfnuw1/5n1AV/XsHIXM2MES6agepXQ0T9FZwcsCZ8ycOB5/x+Qp0AZRMvwi8iUR2SMiu0Vkg4gEcfWGT3J8DG3dRuMf84Ty+EF7Yud/S1fPfOX7Oukqe6739eRcPWHrjIeVHwu9L1vumbEusL5vM+9ibWCPvBLupxg93BE2/KVv69IWy28Mvk7mTH2V1d028vsfbSZRgTaIguEXkQLg88AqpdRiwAkEuZYcPqY08zhhII8fYPZ5MPNMaK3UtXd8vXo7ln/Jh72Pg2FP8BYH0fdtCk+DxMyxqfNHWup5b4OOeFpwZfB1Mor0/USI7DGGf1SIARJFJAZIAiLWTt5IPcNkx0Pwp0uG1pzkTxfr9w+EUtDTEdrjB6/XD/0lmuQp+n5NkEldXzItw++fuOWPM0aXKTjwj4F7+o42f/+Cvk/MgpaK/uUqhoO7B3Y/pY1+fIjJTjuk07cExnB59X/hsSg0OW+rnTTJWxAFw6+UKgd+ApwATgJNSqkX/dcTkTtFZJuIbKupqRny/tISYmnuMIZ/yLz7gE6gGqxh6WrVDU62/H7gdRtP6PLCWbMGXnfGGrj2j3DG5/suP/VjcO09oTV7m2UfgSt/DfkrB173zC8BAvddMbYmMsu26vtVH9d9CR6/ZeRKTFTs1JnLvpm6gchbBI5YKB3BDOcjr8C+v8OJt0dumwPRVqsTuLLnjN4+o0w0pJ5M4CqgGMgHkkWk3yleKXW3UmqVUmpVbm7ukPeXlhhLU8cYrrkylmmv9xqYwU7itVkn66pdUHMw9Lrl2/R94WnhbXvph706vU3GdB32GQ6JmbDylvAiOPIWwi1/06WP77t87MkaM8+AK38Fx96AZ780Mm0jSzZa2z499Hpxyfo3O/r68PdpY5eA2PSrkdvmQFTu0veB8j4mKNGQei4AjimlapRSPcCTwABH2NBJN4Z/6Bx5RVcshMF7u20+V2l7ngy9btk2LfPkLRrcPkaL/OVwy9PaC77vsugbf1/j7ozVE7Bn/TvseHBkDOaJzXriPGXKwOvOWq+T6UZCanJ1Q0ulrgu0/x9Qe3j42wwHY/hHhRPAWhFJEhEBzgci1uQ0PTGWjh433S7TfnHQHP4XxKfpx0P1+JOyYfeToT3Rsm26LWKoCJtok78Cbn4aOposzz+Kxcm6WryPHdYE97nfhEXXwEvfhr3PBH5fOHg82vDPXBfe+sVnA0pnQQ+X5jK9rTO+qJOpNv96+NsMh6rdOvfDTO5GDqXU28ATwLvALmsMd0dqf+mJ2pg0dxqvf1B4PHDoJV2iODFz8JObrdX6/tSPQe0BXUI5EK5u7THa9fDHMgUr4ZandFG0+y/Xxd6igW/jFYd1snQ44Oq79Pf45J1Q/u7Qtl29V0cJzTwjvPULVunon2NvDG1/vtgn04JTdbbwzg3e4yiSVO7yNpCfJEQlqkcp9W2l1ClKqcVKqZuVUl2R2pdt+I3cM0hO7tTVLudepBN1hir1rLoNxKEbdQeiapcu1xCuvh9tCk6Fm5/S8x/3XQ5N5aM/Bt9qnA6n93FsItywQec0bLhhaCemkk36fkaYHn9MnJ4LGAmd3x5vxnQ4/XO67PY79wx/u6Ho6dQZ4VON4Z9QGMM/RA69BAjMOd/KdB2C4U9Ih/QCLQcEk3vK7IndVf1fG6sU2sa/Tmv+zRGLRg6Mr8fvL4+lTIEbH9NJVY9c31cWCocTmyCtUFdEDZfis/VVXfPJwe3Ln6ZSQPT+c+bqzmlb74lsgljNfp2kN4n0fRik4RcRh4ikRWowkSDNGP6hcfglLW0k5+h47cYT/XvahqK12lsnf/EHdRao3Rzdl7JtunF3WsGIDHvUKFwFNz2pQwHvu3x0jX8fjz/AvEjeQrjuPi3bPHFb30qeoVBKe/wz1w2uZk3xen1//M3w3xOIxlJInaqvIkCH7HY0wI6Hh7fdUFTt1vd5xvD3QUQeEZE0EUkGdgN7ReQrkR/ayNCr8RvDHz5tddogz7lQP8+YqS+7WysHsY1ab1LVKZfrSchAck/ZVi2fjMfiWNNPg5ufhNYqy/gP0+MNlz4af5CaRHMugEt/pFtXvvDN8LZbf1R/loHCOP2ZukT3Rx6u3NN0QvdItpmxFgpXw+bf9C1RMZJU7oLYZG829yQhHI9/oVKqGbgaeA4df39zJAc1khipZwgceQVQWt8Hb4bmYHT+tmpvY+6kLF1uYc/TfeWetjp9JTBe9P1ATF+tPf/WKj3h2zKIk+NQ8Q2dDFaMDuC022HtZ+Dtu8LTym19P9yJXRuHU5e/OPb68PIIGkv7F4Q74/NaZtw3jEilUFTu1ldIvnMlk4BwDH+siMSiDf8zVuz9CGSJjA5pifqPYTz+QXDoRR2Gmb9CP7drsgwmssdX6gEt9zSVehPCwCdxaxzp+4GYsQZu+qv2+O8bBePfR+oJYfgBLvqeLjb33Ffh0L9Cr1uySf/uOfMGP6bi9fr3bRhi+QaPRzdvT/cz/PMvhazZOj9hJJLTfFFKe/yTTN+H8Az/H4DjQDLwhojMBMZgndrAxMc4SYh1GI8/XDxu3d92zgU6RBAsL0zCn+B1dWs5ItknAWj+pbpy5m6fZK6ybTrixz7BjGdmrIWbntBa//1XQEtV5PYVKJwzGA4nfPCPMGUh/OVjUBUkrBb0xO6MQer7NrbOP9SwztYqLSemF/Zd7nDC6f8GFTtGJlfAl6ZSnZQ3yUI5IQzDr5T6lVKqQCl1qdKUAOeOwthGDJO9OwgqduhoFVvfB101My0/fKnHbnqenONdlpAGcy+EPU95JxvLtsKURTr1fyIw83T46F90iOf9V0QuBt3X4w8n6S0+FW58VH/Pj3wk8LiaK/QV3WD1fZucuXqSfqg6v12qIVA00bIbdL+FkS7j0Juxu3RktzsOCGdyN09E7hWR56znC4FbIz6yEcQY/kHgG8bpy2CabtiGxT/lf9E1eoL4xGZ9aV/+7viXefwpOgM++rg2ZJEy/n08/jC16fRCbfzba3WMf09H39d79f0hGn4RHdZ57I3BRX/Z2IbfX+oBnZ+w+k4tQVaPYJJ/5W5AtMY/yQhH6rkPeAFdUA3gIPDFCI0nIhjDPwgOv6SNsX/6+mDa7LXZHr+f4Z9/ic7y3P2kbmTe1TTxDD/ocs83Pq6/r/uvhNahV5cNyEDhnMHIX6ErmJZvh6c+1ddAl2yCuNThhTUWr9cnlpohGGc7azdYt6/TbtfNdzaNYBmHyvche/bEueIcBOEY/hyl1OOAB0Ap5QLGVUsrbfhNaeYBaavVXrgdzeNLxkw9+ebqDmM7lpfrK/WA/oPN+wDs/RucsEr5jueInlAUn6U9/4bj8MCV3pPhSNAnqmeQ9Y0WXA4Xfhf2Pg2vft+7vGSTjlAKFSU0EHb/4qHIPU2lOiQ0PjXw68nZsOImeP/xkQubrdo9KfV9CM/wt4lINlYkj4isBSLU3TkypCXGmqiecDj8MqD0xK4/mTP1a/YleSjscg2BqjsuulZ7hZt+DfHpfdsnTjSKz9ZZtPVHtec/UsY/nDj+UJz+eV2W+s2fwM5HdPmJmn1Dl3lsMqbrngpDmeANFMrpz7rPgnLD22H0eBiIzmZ9Up6EET0QnuH/f8AzwGwR2Qg8AHwuoqMaYdKN4Q+PQy/qEMxpy/u/lmHH8h8feDut1fqyPC5A96a5F+rldYd0ZrAjnENwHDNrvWX8j8ADV+ncheHg8fRttTiUCBwRuOxn+sT0zOfhzZ/q5cM1/KDlnpKNg0+4aiqF9AHKRGQV665g2/48+FIU/lTt0ffG8AdGKfUusB5dM/+TwCKl1PuRHthIkp4YS0uXC7dn3KQfjBzhxpQHCuP0xW58EmyC1+P2atltNfoEEsgoxSbCKZfpxxNR3w/ErHPghkeh7rA2/u31fV9XKvzfqbvF2yNhODhj4boHtDHd/BsdahtOR7KBKD5bN60JVJ4DdFE0/9r9SoXn8YNO6Opqgu33hz+mzub+9X4mYQ1+X4IafhG51r4BVwLzgXnAFdaycUN2sq79UVIXwWJPY5GSzfDT+VC6deB1K3bqP2QgmQd0qJ4zPnh/1a1/hF+vBFeXZfhzAq8HsMTqlBVuBciJwOxz4fpHdFGwN37c97X9z8LPFoZ3NTWSvX8TM/XVSFK2TkKLHaDncTjYOv+RVwO//sr/6F7MvnQ06BNaoIgefwpOhRmnw7uDMPyP3qi7k/lSs0/PKaROC387E4hQHv8VIW6XR35oI8cFC/NwOoRHt0axeUY02Gk1Om8Ko2NUiZUcY/9x/XE4dEZnsHC6429qT6+jUecBhDL8cy+AT76pyzhMJuacD2nT+nu8B57X2rXthYbi3Qf0CfiLu+ALI3DhnTULPvmG7mM8EiTn6Oihwy8Ffr32ENQc6BskUGd12gq35+2sc/R2wqnaqZTu91BzoO/ylkp9ohmPNaJGgKAzQ0qpj4/mQCLJtPRELl40lce2lvKlC+aRGDcJ6nL0dHg7MYWjhx7fOHC7vbyFwbMnK3bq+84maG/QmaKhmDb5kmYAHX7p2xRdKV3jBrQxC0VXC7z3KCy+dnBlkwfCP1t2uMy9SF/VtNf3DwtuqwGU7raVNUsvq7V6MueEOdE/dbHeRvW+geXC9nrtkPhLaS2V4bWWnKCEk8D1rUC30RjcSHLLupk0dfTwzHtRaJwRDQ4+rw940BpnKDxuHV5ZNEBxrikLdUinv8faWuON9uls0h5/4uRpYzcoHDG6/rtN/VHvd2d7vsF4/zHobtUx7WOZORfqeYgjr/R/zY748u1bXHtQt1q0AwgGwu7NHM4Vkl07qK26b3nq1mpIyQtvfxOQsMI5fW5u4BKgKIJjigiri7M4ZWoq920qQY10saexyPuPW/qleE8AwajarSfMBqrKaMc8+9d78Z3Ia62CnrZJ1b90UDhj+xr+o6/p+/TpXs83EErB1nth2rKx36ayYKU+8R/yk3uUCmL4D+lCbOHmEGTM1L2g7Vr6obDnpJTHu2+l9HFqPP7gKKV+6nP7PnAOMCviIxthRIRb1hWx72Qz20oaBn7DeKatTodmLvmQ/oMMJPWEm65vp7b798+t2OF9bHtYSdnhj3cy4XD2lXqOva6b0My9UBvAYE7Jic36ez/t9rGvSzucOkjg8L/6Zgd3tYCrUz/2bVZfezB8mQf0589b5A3JDIVvtdAWK/GrowE8PbrpyyRlKEHUScAIi4Kjw9Ur8klLiOH+TcejPZTIsvcp7VUu/YjOhBxI6jn+lvaiBtJ6U6fpSAj/P1zFDt1mEbR0AcbwB8Ph4/F7PDrZqXi9nl/pbNQyWSC2/lEnvC3+0KgNdVjMvVAn6p30cQpsjxu8Hr+rW3vlgy0FbRv+ga7e7eMRvBVTW6174/EHR0R2icj71m0PcAD4RcRHFgGS4mK4btV0nt9dSVVzZ7SHEznef1zr8XmLdVXMUFJPb7u9MJpv2J5WII/fjgayL62N1BMYZ6z2NkHXiulo0FEqtscbaIK3pUpP1K/4KMQljdpQh8Xs8wHpK/fYhl8cXsPfcExHNA3a8C/Wx7WvZBSI+mOQaXXXsj3+XsNvNP5QXI43jPMiIF8p9ZuIjiqC3LxuJm6lePjtMEIcxyP1R6H0be3ti2ippzNEhY2aA9BRH37W5pSFWuO3Pa3mk/oPNWOd9maN1BMaR4x3ktGO5ik+2xvKGEjn3/GAPlms+sTojHEkSM7WcxGBDH/uAq/BHmxEj03vfNMAOn/DMZi+Rj+2DX5v9Vhj+EMxDahXSpUopcqBRBFZE+FxRYyZ2cmcO38Kj7x9gm7XCGRAjjXe/wsgWt8Hy+MPofHb8fsDRfTY5C3SyTb2H9ee2M1fqeUeW7s1hj8wjhivxn/0dciZr2P7M2bo+Pw6P4/f7YJt9/W9KhgvzL1IVwK1axTZhr/wVGip0DLPkA3/QkCs0spB6G7Txj5nrs4k7+fxG6knFHcBrT7P26xl45Zb1s2ktrWL53aPUnPs0UIpHfJXdKZXr49PDS31lGzS2r19OTwQdiidLfdU7NCX7lOXaMOvLG82MXNon2Gi44jR3rurS0/YzrI6VzmcukRwrV9I56EXdMz7aXeM/liHy9wLAWUV/8Nb0iN/pY6yaS7X0lZqfvCqnMGIS9Z5AKE8fjvLOasYUqZ6Y/lbKiEmQV8NT1LCMfyifOIflVIeQiR+jQfOnptLUXYSD2weRPPw8UD5u7oY2NKPeJfFpwWf3FVKJ27NPCP8SJEpC/S9PcFbsUN7rfEp3gne+PTBlwueLNjhnGVboadde/I22XP6e/xb/6ijfub5lTkYD0xbrj1tO4u3rUY7BNmz9fPGE4OP6PElb1Fow2/PN2UW6wge2/DbMfxjPToqgoRj+I+KyOdFJNa6fQE4OuC7xjAOh3DzuiK2lzSwu3xcVZgOzfuPablg4ZXeZaEmd+uP6o5Yg6nKGJ+qZQk7oqJih7dnboLlQZmJ3eA4nFq+Ofq6vlLynVTPmauNlV3OoO6IToI69ePDq5MfLRwOn7BOt06iSs71Zh03ntAe/1Cau4O+yqw/Bl2tgV+355uy/A1/1aTW9yE8w/8pdGXOcqAMWAPcGclBjQYfOrWQpDjnxAntdPfA7r/qLle25w3a43d3a2nBn974/TD1fZspVmRPc7n24noNv7Vfo+8Hxw7nPPa6/t4SM7yvZc/VUpktUWz7k5aGVt4SjZGODHMv1JFLttafPEVfwYhDX/V0NQ/d8Octord0QyDqj+krjMRMbfjt7N3W6kmt70N4CVzVSqnrlVJTlFJ5SqkblVIR6iI9eqQnxnLNigL+9l4FDW1hdJUa6xx5VcdN+8o84DXGgeSeko3aSOfOH9y+8hZpT83uomUMf/g4Y7UhLNum4/d9sQ1g3SHobocdD+n686nj2Dudda428ode0gY3OUd/B2kFXu1/yFKPHdkTpHRD/VHv3FXqVG/2rvH4w4rjnyciL4vIbuv5UhH5z8gPLfLcsq6IbpeHx7ZNgKqd7z+mPRv/ssr2pFkguadko5Z5Bqt15i3Unumuv4A4raJZ+Bh+I/UExRGjT9DK3VffB8ixQzoPwZ4ndULXWK/LMxBJWVC4WmeSt9V4Pe2MGXrSGobu8WfMsEo3BMngbTimZR7Qk7ug5aWOemP4w1jnHuDrQA+A1YTl+kgOarSYPzWVdbOyeXBzyfhu0tLVAvv/odsaxsT1fc2OXPA3/I2l+k8w88zB72+KFdlz6CUd1x+bqJ8bj39g7FaJMQne+HKbhHQthdQe0pO6uQtGpitWtJl7gQ777WzUGj94df7YZEjLH9p27YTCQCGd7h59jPd6/Fbd/ZPv6Xsj9QxIklLqHb9lw+pcLiIZIvKEiOwXkX0iErWOHLeePpPyxg5e3lcVrSEMn31/B1dHf5kHvBOu/lJPuPV5ApE9W1dTVG7IX+5dHm88/gGxDf/0II1PcubCwef0pPlpt02MyJO5F3kf+xv+nDnD+4x5i7XH7/HLyWkq1cdnlo/UAzpb2vf5JCUcw18rIrPxNlv/EDDcAPhfAs8rpU4BlgFBZmcizwUL8piWnjC+Qzvff0y3Rpy+uv9rwTz+kre0obbj8geDM1aHcIJX3wevx29KMgfHDnOdtT7w6zlzdb2euJTAJ/LxyNSlXqmln+EfosxjYycU+jcbskM57Zr/KVMAMR6/RTiG/7PAH4BTRKQc+CLw6aHuUETSgbOBewGUUt1Kqcahbm+4xDgd3LR2Jm8druVwtV+G64HndEmCgy8OfsMn3obmCu/jRr95BI8b9v+D5o5uXt5XNfRS0c0ndWigXaLBH1vjD+Txz1ynwwuHgn3CCGT4jdQTHNvj99f3bbKtic6lH/FerY13RLxzT74aPwzf8Ns9c/3lHjuU05Z6nLF6YrnX8BuNPyRKqaNKqQuAXOAUpdSZSqnjw9hnMVAD/FlEdojIH0Uk2X8lEblTRLaJyLaampr+WxlBPrxKZ7m+ut9nP5W7YcP18NhH4ZEP61jkwfDIdbDxl/rxozfCC1/v+/rm38CjN7Lj+fu47f5tbDla338b4bD7r4CCJR8O/LptjH09/vZ63fRjxtqh7RO04Uor7HvFkDNXX2EM1H1rMpM9W9een7Y88Osz1kFcKqwe9xHTfVn6YX2FaXvgUxbqY3O4cxh2QmHN/r7L645CTGJfScfX2E/yzPKgWSEi8v+CLAdAKfWzYexzJfA5pdTbIvJL4GvAf/mupJS6G7gbYNWqVRGdec1NiScx1kl1i0/FTvtAsmvSPPdV+PTm/pOngehu0xNZHQ3Q06mjOI6+riec7Et962rA01QBFHDX60dYN3sInvLuv+rmHMFC4nrDOX0S1WyvZ9qywe/PZvkN+uZL9mz4+gSIkIokqz6hE7KC6dqFp+rvcCJo+77MOge+VuL9XMk58LURKJQYl6yvMJvK+i6v2Q+58/p+j77GXiZB+9UQhPL4Uwe4DZUyoEwp9bb1/An0iSBqiAg5qXHUtPgkOdnlcdvr9OV53WHY8rvwNmhX/+tq0ZmxoD3usm3edaxL/q5unUPwxsEa9lQMMou47ghUvBu6RrvDqb3wjkbvMnuCa+owDL9h6Axk1Cea0beJ1OdKL+xv+GsPeuehbHwTGyc5oZqtfzcSO1RKVYpIqYjMV0odAM4H9g70vkiTmxJPTauv4beqBiqPvizNmAmv/wiWXjdw+JldhbCrxdv8AeDIy1pXh15tvau7mymp8bR3u/n960f59Q0rCJs9T+r7xdeGXi8hXV+B2FTu0gk0yUaLN0wA0qdrJ8imq1VH9eR+rO96vlnSk5wBC4CISCHwa8DO638T+IJSqiz4uwbkc8DDIhKHrvvz8WFsa0TITY3nWG2bd4FvsazETPjA9+G3a+Clb8EH/xh6Y74ev10KNjFTZyqeZ+W+WR5/T3c3BZmJrC7K4p43j/LvF81jZna/KY/A7Pqr1oQH6pyVkNHX4z/5vo60MBgmAunTde9ipfRVhe20+WekD0HX7+npoaysjM7Osd24KSEhgcLCQmJjwyuOGE7lpz8DjwD27OFN1rILhzRCQCm1E1g11PdHgtzUeLYet3rxejx9y+MmZup44DO/CK//n9ZoQ9Wvt+t9d7d6C0MtuQ7euVv3w03O7jX83T3dZGbF8Ykzi/nzxuPc8+ZRvnf1koEHXLUHavbBpT8ZeN3EDK/G392uT2qLrh74fQbDeCC9UP/XOhv1f7XmgF6ee0rf9RIyvI/DlJ3KyspITU2lqKiod35zrKGUoq6ujrKyMoqLwyuvHk44Z65S6s9KKZd1uw8d4TOhyE1JoL6tmx63R6eSuzq8L9qewhlfhPQZ8M+v6AqLweiVelq1x++ItaJuFBx9Vb9mG/7uHjISY8lLS+DalQU8vq2s71xDMHb/VU9QLbx64HV9pZ6qPVq+Mh6/YaJgX/HaOn/tAf2f8+8xMQSPv7Ozk+zs7DFr9EHPUWZnZw/qqiQcw18nIjeJiNO63QQE6Qg9fslNjQegrrW7f/s7OxM1Lgku/gFU74Ft9wbfmO3xd7Xox6lToWClV+6BXo3f5eohI0lHCt159ix63B7+vPFY6MEqpQ3/rPWQEsY52FfqqbQjeozhN0wQ0qfre9vw1xzQvQ38S1kPUeMfy0bfZrBjDMfwfwK4DqhEZ+x+iDGgyY80tuGvaenq3wXJ11M45XKYfR688n1vRyF/bI2/p00fjKlTtaGfda6ur65Ur8fvcbvITNK63KzcFC5ZPJUHt5TQ0tkTfLDl23Xp3lDRPH3Gn+GVek6+r08E9p/FYBjvZAQw/LkBEsN8pZ5xxtNPP42IsH///oFXDoNwErhKlFJXKqVyrdLMVyulJlyn8l7D39qpPf5YnwlW3xIEInDx/2mj/vJ3Am/MNvygS8PaiSNzztfhnVV7eg1/DG4ykr25AZ9aP5uWThePhGoGv/uvulbOgsvD+3AJ6Xq87h4d0TNt6cQNGTRMPpJydAOiplKdN9NwrL++D35Sz/g6/jds2MCZZ57Jhg0bRmR74ZRlLhaRn4nIkyLyjH0bkb2PIfp6/Af7RgT4a4O582DtZ3S99LLt/TfW5mP4m8u9lQFnn6fvj7zca/iduHs9foClhRmcMSebe986RpfL3X/bHjfsflIXvgo3Ltn2dNpqdQMVo+8bJhIOB6QX6LIo9Uf0HFagUhDjNJyztbWVt956i3vvvZdHH310RLYZTlTP0+i6On8HPKFXHb/kpGivu6alSydrzTpXJ0dB4GqT6/8D3n8c/vlluP0VffDZtNZoY2+Hctpp42n5Oifg8Mu97RFj8JCR2Dcb+NPr53DTvW/z1LvlXL96Rt/9lmzSVw2LPxj+h7MP+LKt4OocXsauwTAWsZO47Iz7QB7/MKWe7/59D3srgrQxHSIL89P49hWhCyX+7W9/4+KLL2bevHlkZ2ezfft2Tj311GHtNxyNv1Mp9Sul1KtKqdft27D2OgaJj3GSnhhLc2O9Nti587xp3YGiAeJT4aLv6fK5Ox70Lu9q1bJK1mzvMtvjB+31n9jc21fViZuMpL6xt2fMyWZxQRp/eONo/z4Bu5/QMtRgmm/bVwbH39T3U8MIFzUYxhPpMyzDf0B3/Mqe03+dcZq5u2HDBq6/XrdAuf7660dE7gnH4/+liHwbeBHojTNUSr077L2PMXJT43E2WBO72XOtxtju4GWGl3xI90X913dgwRX6ysCWebJn6dLH0Ld13pwLdIG2E5sBcOIhM7mvxy8ifHr9HD77yLu8uKeSS5ZYJw5XN+z9G5xyqY4wChfb0zn2hm4Akh2kro/BMF5JL9QOW9UenWUfqNeBbyXaIcxxDeSZR4L6+npeeeUVdu3ahYjgdrsREX784x8PK9ooHI9/CXAH8EPgp9YtjKyh8UduSjxJzUf1k5x53hK6weJ/ReDSH+sY+Vd/oJfZE7t2FULo6/HPWKerBh7TF01O8fTR+G0uXjyVouwk7nr9iLdk89HXdOG3cKN5bGypp2a/rqbpH+ZmMIx30gsBBcfeDCzzjFOeeOIJbr75ZkpKSjh+/DilpaUUFxfz5ptvDmu74Rj+DwOzlFLrlVLnWrfzhrXXMUpuajwZ7SXa4GcVa6knLiV0Rc6pi3Vf1G336oiZXsMfROqJTYCiM3vDK+PEQ2Js/0qBTodw59mzeb+siU1HrLSJ3U9o7332IL9+30tcM7FrmIjYSVxdTYFDOccpGzZs4Jprrumz7IMf/OCw5Z5wDP9uIGNYexkn5KbGk9d9QnezcsbqS8Nwsv3O/YZe759f8SZvZVuG3xnXfxtzzu99mOhUQS/Zrl1ZQG5qPL9//YgutbD/H3pSOJzS0L74TmqZxC3DRMQ3L2UCefyvvvoqF1/cdz7v85//PHfdddewthvONX8GsF9EttJX479yWHseg+SmxjNTlePKWqy/mHANf2ImXPAdeOZzWopB9MkDdMs5f8M+22v4E5zBWw0kxDq57cxifvjcfkrefpqZ3a2Dl3lAX2U448HdZUoxGyYm6QXex/7lmAMyvuL4R5pwDP+3Iz6KMUJOcgxFUklb6pWkg5Z8wq3vsfwm2PZnHQKalKMbRDhiAzd1zpmrPZSm0pCGH+DGNTP47SuHqd/yCDNT8rRMNBQSM3Qcf57pjmWYgMQm6n6+bTUTSuqJFOFk7r4e6DYagxttpksN8eKiPrEIgHaXok6lhPdmhwMu+wkKoTlGnyw6ncnUOwJEBInQWng2APGO0KkRaQmxfOK0bBa2bqF51uVD7pHbIsm0ps7SfxCDYQLSFJdHR+JUb5/pkES0qd+YJxyNf9JQgK69U45uCP1Ez5k80xV+osTx+FP4pftDPNGj+4g+5lrPUz2rA667waOrWncycP3s25PeIF56+Btnhz0Wf/7SfQYPui8Y8vsNhrHOYz1n8aAr9DF+m/ub/MO9mh41uU2fievzIc+hI20OtidzJvBrx0c5JSY17Ip0//PsXl7uuYZcVzyfAH7QcwPLVAa3BVj38fIslnvmkaw6Arzqg6uL1B13syd+OfceSecmFXwyOBR3q2uob+zmdreHWOfkPugNE5M/dV9IVWsnt/S4SQgQKQfwcs8iXmYRhyd5raqgFkBEXrbu/2/0hhNdYju0x7+nSSd/eDyKpo4QVTJ9eHlfFS/vr2ZKajz1bd10udx0uTxUtfSvkV3R2MGh6lYaVSqpngFSwN9/DFpOUrvsMxyva2dX+SD78lo0dfTQ7fZwuLp1SO83GMYyPW79X1MKjta0DfyGSU4o12+aiJwOXCkiK0Rkpe9ttAY4qrRW0S1x7KnT+p/Lo2ho7x7wbZ09bv772b3Mzk3mjrNm4fYoyhq0J1/V3OlNwLJ446A+wTSoFJJDGX6PGzb+EqYuZfn6a4h1Cs/srBj0x+p2eejo0QXf9oxwrRGDYSxQ2aSNPsCh6pboDmaEcTqdLF++nGXLlrFy5Uo2bdo07G2Gknq+BfwXUAj8zO81BUy8JK7Watpjszla147bo/B4FI3tA3v89751jJK6dh68bXXvFcJxq39vZ4+Hli4XaQleLf+NQ9rwN5JCijuEId7/rC4Y96E/k54cx/p5U3j2/ZN849IFOBzhX6o2+9T231PRxIdOHaBHr8Ewzihv9EqmR8K4qh1PU7uJiYns3LkTgBdeeIGvf/3rvP768OJrgnr8SqknlFKXAD/yydid0Jm7tFbjSsql2+WhvKEDl0fR0unC5Q4eeVPe2MGvXznEJYunctbcXHJTdHln38bt1c1eucfl9vDWoVoSYh00qhTiVBf0BND5lYK3fqHbxy28CoArlk2jsrmTrcfrB/WxfOUq4/EbJiIVluGPj3FwuGbiypnNzc1kZg6+haQ/A07uKqX+R0SuhN6QkteUUs8Oe89jkdZqnGkFUAlHaltxW9eOTR09ZFsG3Z8f/GMfAN+8bAEAOVZd/+N1XsNf1dzFnCk6xOy9siaaO12cMz+XhsNWqGhHQ/8wy2Nv6JyAy3/eG8J54cI8EmOdPPNeBWtmZYf9sZotw1+Qkci+imY8HjWoKwaDYaxjG/7VxVkcqoqQ4X/ua7osy0gydQlc8sOQq3R0dLB8+XI6Ozs5efIkr7zyyrB3G04jlv8FvgDstW5fEJEfDHvPY5HWKhIz8wF9uWiXRG4IIvdsPFzLP3ad5LPnzKEwU1fLzLFOEMdr23vXq/aZ4H3jYA0OgbWzsmlQVrxxewAP/q2fQ/IUWHZj76KkuBguWJjHc7srdVP4MLE9/tNnZ9PS5aK0oX2AdxgM44vyxg5yUuJYXJDO8bq2gP8Pj0+JczWOtB5b6tm/fz/PP/88t9xyS795w8ESTjjnZcBypZQHQETuB3YA3xjWnsca7h5oryMhcxqZSbEcqfEa/qaO/hO8PW4P335mDzOykrjjbG8lzrSEGOJiHP08fps3DtWwtDCDlPgYmrDaO3b4Gf6KnXD0VV0Gwq+87BVLp/H39yrYeLiWc+ZPCeuj9Rr+Odn8ZXsZeyqamZmdPMC7DIbxQ3ljJ/kZiczJTaHHrSipa2fOlL7Jl63druHtZADPfDRYt24dtbW11NTUMGVKeP//QIQb0J3h83h8djMYiLZaQEHKFGbnpvQJe2xo6+/x37/pOIerW/nW5Qv7xAyLCLkp8X0mm6osjb+xvZv3Shs5e14ubo/yevwdDX03vvEXEJ8Gqz7Rb7/r5+eSlhDDM++FH93T3KkP+FUzs3A6hD0VQwsJNRjGKhWNHRRkJDI3Txv7QGHLTWEEaox19u/fj9vtJjs7fKk3EOF4/P8L7BCRV9GVjc4GvjasvY5F7KqaKXnMzk3h+T2VvS81+sXyV7d08ot/HeLc+bmcv6D/WTcn1Wv4c1LiqG7RHv9bh2vxKFg/L5edpY002OUgfKWeuiO62crpnw/YMSg+xsnFi6fyz12VdIZIVPHF1vhzU+OZOyXFTPAaJhRKKSoaO1g/L5fZubbhbwH61skKNydnrGFr/KA/6/3334/TObTSLTbhTO5uEJHXgNOsRV9VSlWGeMv4xK6jn5LHrNzkPgdJo18s/w+f20+3y8O3rlgUMIs2N8VbNrk4J7k3queNgzWkJcSwrDCd7SX1NGJP7voY/k2/1sXd1n466FCvWJbP49vKeO1ANRcvnhZ0PZumjh7iYxwkxDpZmJ/GW4dqB3yPwTBeaGzvob3bTX5GIsnxMRRkJAb2+H3+02ocBXS63e4R32ZYUo9S6qRS6hnrNvGMPvh4/FN6vQYb31j+bcfrefLdcu44u5jinMA6ea4V2ZMQ6yA/I5Gq5i6UUrxxsJYz5+YQ43Tg8ii6iKPHEe/1+FsqYefDsPyGwFU9LdbNyiYnJS5suae5o4f0RJ1HsHBaGtUtXbqpvMEwAbCvrgsy9HzYnCkpHBrA8E92TNEWG9vwJ09htt+kkJ296/YovvW3PUxLT+Cz5wZo5mxhR/akxMcwJTWequZODlW3Utncydlzc/W23Nrj6IhJ92r8W+4Cj0vLPCGIcTq4dMk0Xt5XTWvXwBNWTR09pFmGf1G+lo+Mzm+YKNihnPkZOiR6zpQUjtS09oniAcJKxpwsGMNv01qtNfXYBKZnJhLr9Eo4tsa/4Z0T7D3ZzDcvW0BSXHCVzNfw56Ul0OXy8KzlnZ89Txt+l3VQdsZahr+zSTduX3iVt3tXCK5clk+Xy8NLewe+AGvy9fjz0wCTyGWYOFT0evza8M+dkkJnj6dPgAUYj9+XkIZfRJwisn+0BhNVWqsgJQ/QHrVvuGNjezcNbd385MUDrJuVzWVLQuvqttSTHB/DlDR9+fnE9jLmTEnp9UrsUNGumAwt9Wz7E3Q1wxlfDGu4K2dkUpCRGFbtnuZOr+FPT4xlelYie43hN0wQyhs7iI9xkJWs59bsME5/nb+Pxj8IiX+4MfOjwWDHGNLwK6XcwAERmTGcQY0LWqt7DT/A7Fxfw9/Dj188QEuni+9eFXhC1xd/qQegoqmT9Za3D16Pvys2DVpOwubfwaxzIX95WMN1OITLl07jzUO1NLSFLiTX1NFDWoL3CmXRtHQj9RgmDBWNnRRkJPb+L23D71+sLVA+zkAkJCRQV1c3po2/Uoq6ujoSEhIGXtkinHDOTGCPiLwD9GYlDbfnrog4gW1AuVLq8uFsa0RorepjdPUEr9b9j9a0sfdkM584o5h5eQN397E9flvqsTnbx/C7PTqzsDsuA2pL9MIz/zioIV+xLJ8/vHGU53ZXcuOa4Ofmpnavxw+wKD+N5/dU0tLZQ2rCwI1gDIaxTHljR++VNEBGUhw5KfEhPf5wKSwspKysjJqammGPM5IkJCRQWBh+8cVwDP9/DX04IfkCsA9Ii9D2B4efxz/LiuzJTIqlob2HnJR4vnDB3LA2lWOFc6YkeD3++BgHa4q9bRhtj787LkMvyF8BxYPrsLUoP41Zuck88155UMPv8Shaulx9DL+t8+872cJqnzEZDOORisYOzpmf22fZ3ACRPUOZ3I2NjaW4uHhY4xuLhNVzFzgOxFqPtwLvDmenIlKILgUxOBc3UnQ2QXeLbtZsYUs9tmzz9UtO6VNaORQp8TEkxDpIiY8hOT6G1PgYVhdn9Um2SrQeq0TL8J75JRhkVyAR4Yql+bx9rJ7Kpv4NXwBaulwoRW9UD3gje3YPsakLwKsHqrn53rd75ypAVyG9+rcbg47FAM/vruT2+7cFlQ4OV7dy1W/e4mTTAJ3Zxhm7y5u45ncbabFKhLd09nDt7zay7+Tw5pq6XG5qWrv6ePyg5Z7D1a19vmczueslnCJtdwBPAH+wFhUATw9zv78A/gMIWmlMRO4UkW0isi3il1m7/qLvi87sXbS0MIOvXXIKP//Icr556QKuWVEQ9uZEhB9eu5Sb180E4H+uXsx/fOCUPut87ry5/MfF81ly0a1w4X/DKUNTu65ano9S8PcgMf32H833pJWXFs+09AS2n2gI+J5weGF3JW8equWITwncPRXN7Cxt5P2yxiFvd6KzvaSef+2rChhnDrrw33tlTTy0pWSURxZZnthexo4TjRyxumPtr2zh3RONvH5weP/tI9VtKEW/3JvinGRaOl3U+8x/1bZ6c1eaJ/lJIJxwzs8CZwDNAEqpQ8CQqwOJyOVAtVJqe6j1lFJ3K6VWKaVW5ebmhlp1eCgF79wD05ZD4Wm9i50O4VPrZ7O4IJ07zp416DLGV68o4JSpab2PlxT2Lb+QGOfkM+fMISZzOpzxhd7Sy4NlVm4Ky6Zn8NSO8oCvt1h1elJ9JndFhNXFWbx9tH7Ik1b7K/XE2ftl3qsGu+GLiZcOji3xbTwcOHv6qHUifWxrKV2ukc/YjBavHtCZ8XWW8S2zKsSG0zQlFPYVw4JpfefeZmbrarkl9Xo/bo+ipqWLudbEb/UkT2AMx/B3KaV6T5siEsPwGticgW7neBx4FDhPRB4axvaGx/E3oWY/rL5z0FLLWOGa5fnsPdnMgcr+Lee8hr+vTLWmOJva1q4+DWPCxeNRHKzS+9rl493bxeDqw2hXOVlxuW3DXxfw9aO1bSTEOqht7eb53RMjSf5YbRslddoA2153Wb2Wso4Ms2nKvpPNxMc4KPKrNmsb/hPWfuvauvAoWFygHbDJnrkejuF/XUS+ASSKyIXAX4C/D3WHSqmvK6UKlVJFwPXAK0qpm4a6vWHzzt2QmAWLr43aEIbL5cvycTqEp3f29/ptqScloe88vj2p+/axwXXzAihtaKe9W3uj7/vME9iXzwOFl05mXFY019tH6wJ2djtS3coHFk2lKDtpwsg9r+6v7n1c26qPDbsn9ZGatmGFSu6vbGFeXioxzr6mrDAzCRF6TzjVVmn0RVZggzH8A/M1oAbYBXwS+Cfwn5Ec1KjRWAr7/wErb+nfAWsckZMSz9lzc/jbjvJ+aep2SYdUP8M/OzeZnJQ43hmC4bdlntOKMtlb0dzb9MKWeuqN4Q+K7fG3dLn6nDQB2rtdVDR1Mjs3hZvWzmTr8YZhT36OBV47WMOs3GRS4mO8Hn+jNshNHT3UDfF4UUqx72RzP5kHICHWybS0BErq9RWt3QzJDmyoaTWGPyRWA5b7gf8Bvgvcr0Yom0Ep9VpUY/i3/Unfn3Zb1IYwUly9ooCKpk7e8evH2xxA4wevzj8Uw3+gsgUR+ODKQrpcnt5Wd7asFKxjmUFr/PZvsclP57dlt1m5yXzo1ELiYxzj3utv73ax5Wgd586fQk5KHHU+Hn9mkpYfh6rz17R0UdfWzYJpgSPCZ2Qn9Xr8djOkmdlJpCbEGI9/oBVE5DLgCPAr4DfAYRG5JNIDizg9nfDu/TD/UsgY/4nJFy2cSnKck6f9JnltqSc1vn8o6pribMobOyitH1wrxgOVLczISurt+7urvBHwkXqMxh8Ul0eRmxLPgmlp/XT+o1bEy+zcFDKS4rhyWT5P7Sjv/Q3HI5uP1NHt8nDu/Clkp8RT29qF26Pr59sJjXakz2DZZ1152kEU/szMSvYx/Nrjz0mJJzc13hj+MNb5KXCuUuocpdR64Fzg55Ed1iiw5ylor4PVd0R7JCNCYpyTDyyeyj92naSzxxsN0tLpIsYhJMT2/6ltnX+wXv++ymbm56UyM0t7T3Zkj311YTT+4LjcHpwO4YzZ2WwvaaCj2/tbHa1pQ4Tect83r5tJe7c7aMTWeOC1AzUkxTk5rTiz1+Ovbumkx61YVZRFQqxjyBO8tgy2MITHX9vaRVuXi+qWLrKT44iLcZCbYgx/OIa/RSl12Of5UaB/+Mh4Qil45w+QMx+K10d7NCPGNSsKaOl09ZlMa+10kZoQE7C+0Py8VNITYwdl+Dt73ByvbeOUqak4HMKSgnR2WVq17fGbqJ7guDyKGKeDM+bk0O32sK3E+90fqWklPz2xN9FvaWEGywrTeXBzyZiuFRMMpRSvHqjm9Nk5xMc4ez1+e2J3RlYSs3JShmz4959sJj89gfSkwImVvZE99e1UN3f2FkzMTY3v1fwnK0ENv4hcKyLXAttE5J8i8jERuRUd0bN11EYYCcq3Q8UO7e2P0xDOQJw+O4cpqfE86eMhhqrH43AIpxVl8faxwKGFgThc3YpHwSmWl7WkMJ19J5vpcrl7J3ebOnr6ZPQavLjcHmKden4lxiF95J6jta39ekHctHYmh6pb2XJ08HMx0eZITStlDR2ce4qWdHKS46hv7+6VXwozE5k9ZeiGf9/Jlt7jMBAzs/SVU0ldO1XNXb3lU6akJhiPP8RrV1i3BHS1svXAOegIn/EbAgM6hDMuFZZdH+2RjChOh3DV8nxeO1DdK7e0dLpIiQ9ekmlNcRbH69p7NdCBsCN65k/VkRRLCzLocSsOVLb0Tu4qZdLjg+HyKJwOITk+hhUzMth0RE/wKqU4VtPGLL+ublcsyyc9MXZcTvK+dkBn5Z4zX+d75qTGo5S3VEhBRiKzc5Mpa+joI0+GQ5fLzZGa1oARPTYzej3+NqpbOslL04Y/NzWetm43bWE0MZqoBDX8SqmPh7qN5iBHlNZq2P0kLL8R4geutDneuHpFAT1uxT92nQS04feP6PFlzazBxfMfqOybMLPUykh+v6yJ5o4epqXry2kT0hkYl1sR69B/u9Nn57CrvImm9h6qmrto63b3KQcOOizxulWFvLCnMuyT81jh1QPVzMtL6W2Qkp2sDe+O0kZyU+NJiHUyOzcFpRh0IuHh6lZcHhV0Yhd074nMpFiO1bZR09LFlFSv1AN9SzhMNsKJ6ikWkZ+JyJMi8ox9G43BRYTt94OnZ8JM6vqzcFoa8/JSeqN7mgcovbxwWhrJcU7eCVPu2V/Zwty8FJxWCYvCzEQykmLZdryeLpeHGVnay/JvUG/QuDye3u/uzLk5KAWbj3prHs3yqzkD8NE1M3F5FI++UzqqYx0ObV0u3jlWz7nzvdVd7Kq1+yqaKczUJwO7xs5g5Z59J/WVZ7BQTpsZ2cm8W9KIR9HH44fJncQVzuTu0+jqnL9GR/jYt/GHu0fH7s8+D3LCK7E83hARrl5RwLaSBk7UtdPaFdrjj3E6OLUo/Hj+/ZUtfbwsET3Bu+mIPnHYVwLG4w+MntzVhn9ZYQZJcU42Hq7rrdHjX2wMoCgnmbPn5fLIOyW9yXJjnY2Ha+lxK9b7lEvOtirddrs9FGZqB6E4JxkRbyhruOy3SjUU+0lj/szMSuKAVV7EntydYgx/WIa/Uyn1K6XUq0qp1+1bxEcWCfb/A1oqdF2eCcxVy3Ul0b/tLB9Q6gGt8x+sah3QWNe3dVPT0sUpU/tKZEsL03uLXtm6qonlD4zLrYi1ygvExThYXZzFxiO1HKlpIznO2euV+nPz2plUNXfx8r6q0RzukHn1QA0p8TGsmunt95Cb4v1stvyTGOekICNx8B5/ZTPzp6b2Xj0Fw47sAa/B7/X4jdQTkl+KyLdFZJ2IrLRvER9ZJHjnbp2sNfeiaI8kohRkJLKmOIundpQP6PEDvQ1iBvL691fquOn5foZ/SUFG72Ovx28mdwNhT+7anDE7h6M1bWw+UkdxbnLQtp7nnTKFgoxEHhwHk7xKKV4/UM2Zc3KIi/GamLTEGGKdXonQZnbu4CJ7dKmGFhaE0PdtbOkR6O2Gl5kUh9MhvfV7JiPhGP4lwB3AD/HKPD+J5KAiQuVuKNkIp90+5BLI44lrVhRwtLYNt0cN2F5xaWEG8TGOAQ3/Ab+IHu/7vSWn89LiiY9xGI0/CHY4p80Zc3IAOFDVwqyc/jKPjdMh3LhmBhsP1/VrKTjWOFjVSkVTZ7+uWCLSO8Hbz/BXt/WrMxWMmpYu6tu6Q0b02Mz0qdppe/pOh5CdHGekngH4MDBLKbVeKXWudTsv0gMbcbbeAzEJsOLmaI9kVLhkybRebytUOCdoyWHljMwB4/n3n2whKzmuzyU7wLT0hN6Ju7TEWLKS44zGHwTt8Xv/dqdMTSUrWX93gfR9X65bNZ1Yp/Dw22Pb63/Nqr1/zvz+bTuyrePE1vgBZk9JpqPHTWWYUUt7rYzdUDH8NrbUk50c1yuxgT4JGKknNLuBjAiPI7J0NMD7j8OSD0PS5Ogxm54YywUL9B9vIKkHdFjn3pPNvUlYgdhf1cIpU1P7yRH2BC/oTl8ZSXEDavyTNYba5fEQ6yP1OBzCutm65tGs3NATlbmp8VyyeBpPbC+jtctFt2tkJno7e9wjti3QYZwLpqUx1Qrt9cVuZerv8UP4kT12Lkk4Us+U1HgSYh29E7u+y43HH5oMYL+IvDBuwzn3PQs97ROiCudguGZFIeD9s4VidXGWrmQRJENUKcWhKl37PBArZmQS6xTSE2PJSo4NWWp3Z2kjS7/7IjtLGwf+EBOIg1Ut1LZ0kxjXV2o8e66We/wltEB8dM0MWjpdLPvui8z7z+eGPabG9m4u+9WbfH7DjmFvC3QZ8G3HG/rJPDb5GYnkpyf06T9tn/DCrdK57Xg9BRmJQUs1+CIizJmS0udEAzp79+Qk7g09sCsI3474KCJNi05mIm9xdMcxylywYAob7ljbO3kbilNnZpIU5+T1gzVcsDCv3+u1rd20d7uDhs/dflYx6+flkhjnJDs5nrKGxqD7+uv2MtwexcbDtSyfnhHuxxnXHKpq4cZ7tpCaEMPtZ83q89oHVxZSlJ0c9KTqi31yGImSGN0uD596aDtHatqoau7C5fb0a2gyWN45VofLozjLmrvw58sXzeO2M4v6LMtN0V55eePADeaP1bbx8v5qPnPO7LDH9JsbVvaZZAbdjP2xbaXUt3X3Sm2TiQEN/7gN3fSloxFik8E5sIcwkRDxyggDER/j5PTZ2bx6oBqlVD85p9TqkTo9K3C1jqS4GJZZRjw7JY7aIJfRLreHf1pZxTtONIY1tvHOoaoWbrhnCyLChjvX9jt5xjgdvSWuB2KgifpwUUrxzad2seVoPRcvmsrzeyrZd7KlX2/owbLpcJ2eM5qZGfD1nJT4flegIkJ+RmJYhv+Pbx4l1uHg1tOLwh5TUQBnxU782neyuXeCfTIRTuZui4g0W7dOEXGLyPhqC9TZBIkZ0R7FmOec+VMoa+gIWB/drtk/3WdSLhg5KboWim/JYZstR+upa+smJyWenaUN47Lq5GA4XN3KDfe8rY3+HWsHnMAdCKdDSPOZsxnq93fX60f4y/YyvnD+XL595UKAfk18hsLGI3WsmpnZR8oJh4KMRMobQhv+2tYunthexrUrC3rLLwwVOyJob8X4MmUjRTgduFKVUmlKqTR0cbYPAr+L+MhGks5GSMiI9ijGPLYua0dl+GIb/sIwDL8d9ROoFsqz71eQHOfkU+tnUdva3VuidyKijf4WADbcsYY5U4Zn9G18tW3XECSff+46yY+eP8CVy/L54gVzmZaeyIyspLDLdgSjrrVryB50YebAHv8Dm0vocnn6SWVDITslnry0+AnR2nIoDErQU5qngQ9EZjgRoqMREoZ3CTsZKMxMYu6UlN6qir6U1neQkxLfb2IyEHbInr/h73Z5eH5PJRcuzGOtJW3smKATvEdqtNFXSllGf+QKAqYneg3/YLX+naWNfOmxnayckcGPPrS0V9I7rSiLrceHdwVml44OV170pSAjkdrW7qBVOju63Ty4+TgXLMgbsRPowmlpvaGhk41wpJ5rfW4fEpEfAuNrOryz0Ug9YXLO/FzeOVbfL9yytKGdGUH0fX9sDdfur2qz8XAtje09XL40n/lTU0mIdbBzAur8R2taueHuLXg8ig13rGVuGJO2g8HX8A+mdk95Ywe337+NKWnx3HPLqj5yzJriLOrbuodcGx9g45FaUuJjWFoweCerwIq6Ceb1/2V7KQ3tPXxy/fC9fZsF09I4XN1Kl2twJaEnAuF4/Ff43D6A7r51VSQHNeJ0NBqpJ0zOnT+Fbrent+iaTWlDO9OzBpZ5ILjH//f3K0hNiOGseTnEOh0sKUhnR2nDyAx8jHCsto0b7tmC26N4JAJGH/oafpc7PA+9pbOH2+7bSpfLzZ9uPa23YJrNacWDK88diM1H6lhTnDWkyKD8dMvwB5D+3B7FH988xooZGawKMmk8FBbmp+HyKA5Vje1M6EgQjsbvW4f/DqXU95VS/UXgsYzx+MNmVVEWyXFOXvXR+V1uDxWNnWFN7ILX4/c1/J09bl7aU8XFi6YSH6M9zRUzMtlT0TxhPK7jtW3ccPcWetza6IcTlz8U+hj+MKQel9vD5zbs4FB1K3d99NSAJ6Oi7CRyU+PZOkTDX9HYwbHaNk4fYoRMKI//+d2VnKhv55Nnzwpay2go+Eb2TDaChnOKyLdCvE8ppf4nAuMZedw90N1qPP4wiYvR/WBfP1DTG9Z5sqkTt0cFDeX0JyHWSWp8DLU+Us8bB2to6XJx+bL83mXLp2fQ7fKw72TLuI/nL6nTnn6Xy82GOyNn9AHSE71x5y7PwFLP/zy7l9cO1PCDa5Zw5tzAhllEWD2I8tz+bDysO4mdPgR9H2BqWgJOh1DhZ/iVUtz9xhGKspO4cOHUIW07GEXZySTGOielzh/K428LcAO4DfhqhMc1cnTqNm/G4w+fc+ZPobyxg0NWJuVgQjltslPi+nj8f3//JJlJsX0Mw4oZGQDsPDG+5Z6Sujauv3sLnT1uHr59bciuUCPBYKSe+zYe4/7NJdxxVjE3rpkRct3VxVlUNHVSZuVsDIbNR+rITo5j/hClrRing6lpCf2knreP1fNeWRO3nTVrwBLMg8XpEOZPTZ2UHn+o1os/tW/A3ehQzo8DjwIjN8MSaToa9b3x+MPGP6zTm7wVvuHPSYnvndzt6Hbz8r4qLl48rU+hrGnpieSlxY/ryJ4Tde3ccPcWOiyjvzA/skYfwpd6Xt1fzX8/u5cLF+bxtUsWDLjd1WGW5/ZHKcXGI7Wsm52NYxjGuSAjkTI/j/+eN46SlRzHh08tHPJ2Q7EwP419J1smfD6JPyE1fhHJEpHvAe+jZaGVSqmvjiuNv7NR35twzrDJz0hkfl4qr+7XYZ2l9R04HdLbTzccfD3+V/ZX097t5opl0/qtt2J65rit2VNa384N92yhrdvNw7evGRWjD/4ef2CpZ9/JZv7tkXdZMC2NX16/PCxveX5eKmkJMYM2/HbJh9NnDy8DtiCzbxLXoaoWXt5fzS3rZg46ISxcFkxLo6mjh4pJVrcnqOEXkR8DW9FRPEuUUt9RSo2/a3Lb8BupZ1CcMz+XbSX1tHT2UNrQTn5GwqCiNXJS4nsN/9/fqyA3NZ41xf313+UzMiipa6dunJXILa1v5/q7t9Da5eLh29ewKH/0HIuBPP7q5k5uu28rqQmx3HvraSTFhVOSS1cKPa0oa9AZvJuPaH3/jDlD0/dtCjISqWzu7D2Z3fPmURJiHdyyrmhY2w3FQnuCd5Jl8Ib6J38ZyAf+E6jwKdvQMq5KNhipZ0icM38KPW7FxsN1lNa3D0rfB234G9p7aGrv4dUD1Vy6eGpAr3OFNan7XlnjCIx6dChr0J5+S2cPD9++hsVDiFsfDqE0/o5uN7c/sI2G9h7+eOuqgKWRQ3FacRZHa9oGVbJ44+E6CjIS+3S7GgoFmYm4PYqqli6qmzt5ekcFHz51ekSLqOky40y6Cd5QGr9DKZXoW7LBuqVa5RvGB8bjHxKrijJJiY/h9YPVlDZ0DMHw6z/r49tK6XJ5uMInmseXJYXpOB0ybgq2lTd2cMM9W2ju6OHh29eOutEHf4/fK/V4PIr/9/hOdpU38asbVgxpbLbOvy1Mr9/jUWw+Wsfps7OHHWpp9+Etb+jgz5uO4/J4uP2s4mFtcyCS42Moyk6edBO8w6vBOh4wHv+QiHU6OHNODi/traampSvsUE4bO5b//s3HmZaewMoZgRNvkuJimJ+XOi50/orGDq6/ezON7T08dPuaYVeyHCrBavX8+MUDPLe7km9euoALA5TWDofF+ekkxDrCTuTae7KZpo4eTh+mzAN6bgl0G8qHtpRw8eKpfVonRooF01KNxx9pRGS6iLwqIntFZI+IfCGiO+xs1C0XY4dXzW8ycs783F6dfjARPQA5Vn/TsoYOLlsyLWS0x/IZGew80Rh2z9VooI3+FhrbenjotjUsLcyI2lhSfVpp2iUbHt9ayl2vHeGja2Zw25lD95LtNpzhTvB64/eHX9rY9vjvevUwLZ0u7jw7/Jr7w2HhtDRK6tppnURd4aLh8buALyulFgJrgc+KyMKI7c2Uaxgyvj1TB2v4s3102WAyj82K6Rm0dLk4Wjs2U+dPNml5p6GtmwduW93bdyBa+J5E3R7FpiO1fOOpXZw1N4fvXLlo2JLL6uIs9lWGbsNps+lIHXOmpJCXNnzHSjfxiaOiqZPVxVmjltRnZ/Dun0Re/6gbfqXUSaXUu9bjFmAfUBCxHXY2mlDOITI1PYFTrAzUQWv8lsc/IyuJpQNIInYi17th6vzvlTbS1N7XKDV19PBemHJRl8vN1jA17KrmTm64ewt1rd3cf9tqVgSRrKLFwapWPvXgdopzkvntR1f2yZMYKquLdBvO7cdDB/F1uzy8c6x+yNm6gbBLN3zy7NFLFbLDcCeTzh9VjV9EioAVwNsBXrtTRLaJyLaamv5lgsOms8kY/mFwxbJ8pqUn9E7WhktqfAxT0xK4blXhgB7orJwUEmIdHKpqGXC7h6tbueZ3G/nta4f7LH9g03Gu+d1GTtQNnHX61LvlfPj3mzlcPfD+fvXyIU42dXL/J1YHnaeIBnlp+sR69xtH6HEr/vSx00gboe5cK2ZkEhfj6FOvKRAHKlvo6HH3TgiPBIvy01iUn8a5PlebkWZqWgKpCTG9meqTgagZfhFJAf4KfFEp1e9Uq5S6Wym1Sim1Kjc3cOPmsOhqhfjI1U2Z6Hx6/Wxe/fdzBi0fiAivfeUcPnPOnAHXdTiE5LgYOoLUYvfl5y8dxKPo592fbO7Eo+CBzccH3IZdeti/Aqk/Sile3lfNufOncOoIVoUcCe699TQAalq6KMhMHLQUF4rEOCcXLszj7+9V0O0KXgvohFXKI1gf5qHwvauX8NdPnz6sDODBIiLkpMRT39Y98MoThKgYfhGJRRv9h5VST0Z0Z91tED8yjRsmIw6HDDlrMiHWGfYfOD7GQVdP6IJju8ub+MeukyTHOdlb0dxnMrjeKg/x+LZS2rtDT9LZBmvzAIZ/78lmKps7OW/B6Hmf4ZIQq/+6HkWfVowjxYdWFtJg5WAEYyilPAbCOYzjbThkJsXS2D7wnMZEIRpRPQLcC+xTSv0s4jvsboU4Y/jHOnExDroHaCrykxcPkJEUyxcvmEdLl6vX8ADUt3WTkRRLc6eLp3aUh9xOiSUHvX2sPmQk0Sv7tNEbTdkhXGIc3r9uWuLISDy+nDU3h5yUeJ58tyzoOqX17aQnxo6YxBRNspLjjMcfYc4AbgbOE5Gd1u3SiO2tyxj+8UB8jDOkx//OsXpeO1DDp9fP7m3tt7vcqxDWtXVx+uxsFuWncf+m40GLbimlKK1v7/2jHwyh87+8v5pl0zPITY0Puk608M2CjoThjXE6uHp5Pq/sr6YhiEEsa+gYdH7HWCUjKY6GdmP4I4ZS6i2llCilliqlllu3f0ZoZ9rjN1LPmCc+1hG0IYtSih+/sJ+8tHhuPb2IuXkpxDqF3RVNvevUtXWTlRzHx04v4mBVa1AZp76tm7ZuN9eu0IFkwdaraenivbJGLjhl7Hn7QJ/onbTEkZd6AD54aiE9bsUz71UEfL20YfClPMYqWcnG8E8cXJ2g3BAX+ew/w/CIj3HQFWQi8bWDNWw93sDnzptLQqyT+Bgn8/JS2V2uDb/L7aGxvYes5HiuWJZPVnIcf950POC2Six9//Q52RRmJrLlaGDD/9qBapRiTOr7ADHOyHr8oOPbF0xLCyj3eDyKsoYOCjMniscfS2ePh47uidENbiAmtuHvssKz4kxUz1gnPsYZ0PB7PIqfvHCAGVlJXLdqeu/yxfnp7KloRilFgzUpl50cR0Ksk+tPm87L+6p6G8j4Yi+bkZXEulnZQXX+V/ZXMzUtobd641gjNsIav80HVxbwXllTv9DXmtYuul2eEZ3YjSZZSTpcebJ4/RPb8Hdbht9IPWMe7fH397ae213JnopmvnThXOJivIfr4oI06tu6OdnU2TspZzd5v2ntTESEh7aU9NuePbFbmJnEutnZNLb3sK+ybzRxt8vDGwdrOG/BlBHt8TqSOEfB4we4ankBTofw13f7TpgPpSvbWCbTyjSfLBO8k8Pwm8ndMU98bP9wTpfbw09fOsC8vBSuXNY3uXuRVXlyd3kTdW26npBdvjc/I5EPLMrj0a2l/S7dT9S3MzUtgYRYJ2tn6UniLUf7ZvG+c6yetm43549RfR8gxndyN0IaP0Buajzr5+Xy1LvluH2ujLyhnBND6sk0Hv8EolfqMRr/WCeQ1PPkjnKO1rTx5Yvm96vlv2BqGg6B3RXNXo8/2Rt9c+u6Ipo6enh6Z19P9URdOzOytZean5HIzOykfhO8L++vIj7GMSKFxyJFn8ndCIdTfnBlIZXNnWyyGq6A7soG+sppIpCVrL/DhkkSyz+xDX+31R/eZO6Oefylni6Xm1/+6xDLCtO5KECJ4cQ4J3OmpLC3oqnX8Ps27FhdnMWCaWnct7FvaOeJ+vY+DUO0zl/X683a2bpnzMkhMW70E4nCxfc8GEmNH+D8BVNIS4jhSR+5p6yhndzU+KgkW0WCDNvjN1LPBKDbmpAyUs+Yxz+q55G3T1De2MFXPnBKUJ19cX46u8ubqbWydjN96tSLCB87fSYHqlp6pZzOHjeVzZ19DP/aWdm0dLp6C3QdqWnjRH07541hmQf054u1dP5IZO76khDr5PJl+Ty/u7K3dHFpfQfTJ0hED0BGou3xG8M//ukyk7vjhfhYbwJXW5eL3756mHWzskP2cV2Yn0ZlcyeHqlrISIrt1xP4quUFZCTFcr8V2llm6dIzs308fisZzJZ7XtlfBTDmDT94s3cj7fGDlns6etw8t+skoDX+iSLzgE5YS0+MNR7/hMCWeozHP+axpR6lFPdtOk5tazf//oH5IaNq7NaCbx2uDdiXVYd2zuDFvZWUNbT31ujxDUHMS0tgVk4ym614/pf3VbNgWlpvN6ixjB3Lnxphjx9g5YwMinOS+eu7ZbjcHk42dU6YiV2bzKRY6o3GPwEwUs+4Ic7pwKN0Bu7vXz/CBQsGrohp11Fv6XSRkxy4rMJNa2cA8NCWE72hnDP9Ys/Xzs5m67F66tu62VbSMKajeXyJcQgJsQ7iYyKvs4sI164oYMvRerYeb8DtURMmlNMmMzmORiP1TAC6WsEZBzGDqyVvGH3irWqTv3nlMK1dLr580fwB35OWEEuRJdsE8vhBR51ctHAqj249wcGqVpLjnP3WXTsrm5YuF7979TBujxqz2br+xDgdo1og7WqrzMWvXj4EjGxVzrFAVtLkKdQ2sQ1/d5sJ5Rwn2F7rQ1tKuHJZfm87vIGw4/mzQjSKufX0Ihrbe/jru2XMyE7uJx+tnaUbiTywuYTs5DiWRbGf7mCIdcio6Ps207OSWDsrq1cWm2gef0ZSnNH4JwTdraZcwzgh3srKVcCXLpgX9vsW52vDnx3E4wdt2OfnpdLt8jAjgC49JTWBOVNS6HZ7OGf+lH45A2MVp1MiHtHjzwdXFgI6nHRaxvD77I4lspJjTRz/hKCrxUT0jBNsqee6VdMpGkRHp8UF+sogmNQDVmjnGUUAfUI5fbG9/vPHicwDul7PaHr8AJcsmUZirJNp6Ykj0t93LJGRFEdHj3vARj6jRWl9O5/bsKM3hHYkmVi/nD9xyZAWuT7uhpFjUX46K2Zk8IXz5w7qfStmZLK6KIvTikL3fb16eQFrZ2Wxfl5gw37NigJWzczk7HnDaPM5ypx3ypRRbxKTEh/DnWfP4vJl00Z1v6OBLfF97x/7gvZzGC1K69u5/u4tvHGwhvKGjhHfvkT7A4bDqlWr1LZt26I9DIPBMMH5v+f3c9drR/jvqxZxy7qiqIzBNvqtXS4evn1Nb9jyUBCR7UqpVf7LJ7bHbzAYDIPgKxfN54IFU/ju3/ey8XDtwG8YYUbS6IfCGH6DwWCwcDiEX1y/gjm5KXzm4Xc5Vts2avseLaMPxvAbDAZDH1LiY/jjratwCNx+/1aaOyMf6TOaRh+M4TcYDIZ+TM9K4q6bTqWkrp3PPbKjTy+CkWa0jT4Yw28wGAwBWTsrm+9etYjXD9bww+f2RWQf0TD6AKOb/WEwGAzjiI+umcnByhbuefMY8/JS+bBP3+fhEi2jD8bjNxgMhpD81+ULOWNONt98ajfbS+oHfkMYRNPogzH8BoPBEJIYp4Pf3riS/IwEPvngdsobh5dQFW2jD8bwGwwGw4BkJMXxx1tX0dXj4Y77tw25rMNYMPpgDL/BYDCExZwpqfzqxhXsr2zmy4+/h2eQkT5jxeiDMfwGg8EQNufOn8LXL1nAc7sr+aXVlyAcxpLRBxPVYzAYDIPi9rOKOVDVwi9fPsS8vFQuWxq6YN1YM/pgPH6DwWAYFCLC969ZzMoZGXz5LzvZXd4UdN2xaPTBGH6DwWAYNPExTv5w8yqykuK444FtVLd09ltnrBp9MIbfYDAYhkRuajx337KKxvYePvngdjp73L2vjWWjD1Ey/CJysYgcEJHDIvK1aIzBYDAYhsvignR+et0ydpxo5BtP7UIpNeaNPkRhcldEnMBvgQuBMmCriDyjlNo72mMxGAyG4XLpkml88YK5/OJfh8hJiecf758c00YfohPVsxo4rJQ6CiAijwJXAcbwGwyGccnnz5vLwaoW7n7jKOmJsWPa6EN0DH8BUOrzvAxY47+SiNwJ3AkwY8aM0RmZwWAwDAGHQ/jJh5dRmJnEVcvzWZQ/do0+jOHJXaXU3UqpVUqpVbm546cBtsFgmJwkxcXwjUsXjHmjD9Ex/OWAb23TQmuZwWAwGEaBaBj+rcBcESkWkTjgeuCZKIzDYDAYJiWjrvErpVwi8m/AC4AT+JNSas9oj8NgMBgmK1Gp1aOU+ifwz2js22AwGCY7Y3Zy12AwGAyRwRh+g8FgmGQYw28wGAyTDGP4DQaDYZIhSg2ufVg0EJEaoGQQb8kBaiM0nLHGZPmsk+VzwuT5rJPlc0L0PutMpVS/DNhxYfgHi4hsU0qtivY4RoPJ8lkny+eEyfNZJ8vnhLH3WY3UYzAYDJMMY/gNBoNhkjFRDf/d0R7AKDJZPutk+ZwweT7rZPmcMMY+64TU+A0Gg8EQnInq8RsMBoMhCMbwGwwGwyRjwhn+ydDIXUSmi8irIrJXRPaIyBeiPaZIIiJOEdkhIs9GeyyRREQyROQJEdkvIvtEZF20xxQpRORL1rG7W0Q2iEhCtMc0EojIn0SkWkR2+yzLEpGXROSQdZ8ZzTHCBDP8Po3cLwEWAjeIyMLojioiuIAvK6UWAmuBz07Qz2nzBWBftAcxCvwSeF4pdQqwjAn6mUWkAPg8sEoptRhdnv366I5qxLgPuNhv2deAl5VSc4GXredRZUIZfnwauSulugG7kfuEQil1Uin1rvW4BW0gCqI7qsggIoXAZcAfoz2WSCIi6cDZwL0ASqlupVRjVAcVWWKARBGJAZKAiiiPZ0RQSr0B1Pstvgq433p8P3D1aI4pEBPN8Adq5D4hDaKNiBQBK4C3ozyUSPEL4D8AT5THEWmKgRrgz5as9UcRSY72oCKBUqoc+AlwAjgJNCmlXozuqCJKnlLqpPW4EsiL5mBg4hn+SYWIpAB/Bb6olGqO9nhGGhG5HKhWSm2P9lhGgRhgJXCXUmoF0MYYkAQigaVxX4U+2eUDySJyU3RHNTooHT8f9Rj6iWb4J00jdxGJRRv9h5VST0Z7PBHiDOBKETmOlu3OE5GHojukiFEGlCml7Cu3J9AngonIBcAxpVSNUqoHeBI4PcpjiiRVIjINwLqvjvJ4JpzhnxSN3EVE0FrwPqXUz6I9nkihlPq6UqpQKVWE/i1fUUpNSM9QKVUJlIrIfGvR+cDeKA4pkpwA1opIknUsn88Enci2eAa41Xp8K/C3KI4FiFLP3UgxiRq5nwHcDOwSkZ3Wsm9YvYwN45fPAQ9bTstR4ONRHk9EUEq9LSJPAO+iI9R2MMZKGgwVEdkAnAPkiEgZ8G3gh8DjInIburz8ddEbocaUbDAYDIZJxkSTegwGg8EwAMbwGwwGwyTDGH6DwWCYZBjDbzAYDJMMY/gNBoNhkmEMv8Hgg4hki8hO61YpIuXW41YR+V20x2cwjAQmnNNgCIKIfAdoVUr9JNpjMRhGEuPxGwxhICLn2P0AROQ7InK/iLwpIiUicq2I/EhEdonI81Y5DUTkVBF5XUS2i8gLdtq+wRBtjOE3GIbGbOA84ErgIeBVpdQSoAO4zDL+vwY+pJQ6FfgT8P1oDdZg8GVClWwwGEaR55RSPSKyC10e5Hlr+S6gCJgPLAZe0uVocKJLEBsMUccYfoNhaHQBKKU8ItKjvJNlHvT/SoA9SqkJ2z7RMH4xUo/BEBkOALl231wRiRWRRVEek8EAGMNvMEQEq/Xnh4D/E5H3gJ1M7JrzhnGECec0GAyGSYbx+A0Gg2GSYQy/wWAwTDKM4TcYDIZJhjH8BoPBMMkwht9gMBgmGcbwGwwGwyTDGH6DwWCYZPx/ug0CrPgc+1UAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Import libraries\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Define the system with two reactions: A -> B and B -> A\n",
    "# The initial number of molecules for A and B are 10 and 0, respectively\n",
    "# The reaction rate constants for A -> B and B -> A are 1.0 and 0.5, respectively\n",
    "x_A = 10 # number of A molecules\n",
    "x_B = 0 # number of B molecules\n",
    "k_1 = 1.0 # rate constant for A -> B\n",
    "k_2 = 0.5 # rate constant for B -> A\n",
    "\n",
    "# Initialize the time to zero\n",
    "t = 0\n",
    "\n",
    "# Create empty lists to store the time and molecule values\n",
    "t_list = []\n",
    "x_A_list = []\n",
    "x_B_list = []\n",
    "\n",
    "# Loop until the end time of 10 is reached\n",
    "while t < 10:\n",
    "    # Calculate the total propensity function\n",
    "    a_total = k_1 * x_A + k_2 * x_B\n",
    "    \n",
    "    # Generate two random numbers from a uniform distribution between 0 and 1\n",
    "    r_1 = np.random.uniform(0,1)\n",
    "    r_2 = np.random.uniform(0,1)\n",
    "    \n",
    "    # Use one random number to determine the time interval until the next reaction event\n",
    "    tau = (1 / a_total) * np.log(1 / r_1)\n",
    "    \n",
    "    # Use another random number to determine which reaction will occur next\n",
    "    if r_2 < (k_1 * x_A) / a_total:\n",
    "        # A -> B occurs\n",
    "        x_A -= 1 # decrease A by 1\n",
    "        x_B += 1 # increase B by 1\n",
    "    else:\n",
    "        # B -> A occurs\n",
    "        x_A += 1 # increase A by 1\n",
    "        x_B -= 1 # decrease B by 1\n",
    "    \n",
    "    # Update the time by adding the time interval\n",
    "    t += tau\n",
    "    \n",
    "    # Append the time and molecule values to the lists\n",
    "    t_list.append(t)\n",
    "    x_A_list.append(x_A)\n",
    "    x_B_list.append(x_B)\n",
    "\n",
    "# Plot or output the results\n",
    "plt.plot(t_list, x_A_list, label=\"A\")\n",
    "plt.plot(t_list, x_B_list, label=\"B\")\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Number of molecules\")\n",
    "plt.title(\"Gillespie algorithm for a simple chemical system\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Gaurateen of this algorihtm:\n",
    "\n",
    "- Theorem: Two random variable can simulate a random walk process. (?)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stochastic reaction networks with $M$ species (each with count $N$)\n",
    "\n",
    "$$\n",
    "\\sum_{j=1}^M r_{ki}X_i \\xrightarrow{c_k} \\sum_{j=1}^M p_{ki}X_i\n",
    "$$\n",
    "\n",
    "We have equation\n",
    "\n",
    "$$\n",
    "\\partial_t P_t(\\bold{n}) = \\sum_{k=1}^K [a_k(\\bold{n} - s_k)P_t(\\bold{n} - s_k) - a_k(\\bold{n})P_t(\\bold{n})]\n",
    "$$\n",
    "\n",
    "Then we can derive the final state of the system with the matrix $r,p,c$. (The matrix $r$ is the stoichiometry matrix, $p$ is the propensity matrix, and $c$ is the reaction rate matrix.)\n",
    "\n",
    "$$\n",
    "c = p - r\n",
    "$$\n",
    "\n",
    "The chemical transformation matrx $\\mathbb{T}$ is defined as\n",
    "\n",
    "A numeric example:\n",
    "\n",
    "$$\n",
    "\\emptyset \\xrightarrow{1} A\\\\\n",
    "A \\xrightarrow{1} \\emptyset\n",
    "$$\n",
    "\n",
    "We have\n",
    "\n",
    "$$\n",
    "r = \\begin{bmatrix}\n",
    "-1\\\\\n",
    "1\n",
    "\\end{bmatrix}\\\\\n",
    "p = \\begin{bmatrix}\n",
    "1\\\\\n",
    "1\n",
    "\\end{bmatrix}\\\\\n",
    "c = \\begin{bmatrix}\n",
    "2\\\\\n",
    "0\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "The joint probability $P_t(\\bold{n})$ can be difficult to handle, since $\\bold{n}$ is a vector of length $N^M$.\n",
    "\n",
    "This quickly become prohibitive for large $M$. (e.g. 10 species with 16 counts each, $N^M = 10^{16}$, here counts means the number of molecules of each species)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can propose a parameterization of the joint probability $P_t^{\\theta_{t}}(\\bold{n})$, and then minimize the loss defined as\n",
    "\n",
    "$$\n",
    "\\mathcal{L}=D_\\text{KL}(P_{t^+\\delta t}^{\\theta_{t+\\delta t}}(\\bold{n})||\\mathbb{T}P_t^{\\theta_t}(\\bold{n}))\n",
    "$$\n",
    "\n",
    "where $D_\\text{KL}$ is the Kullback-Leibler divergence, $P_{t^+\\delta t}^{\\theta_{t+\\delta t}}(\\bold{n})$ is the joint probability at time $t+\\delta t$ with parameter $\\theta_{t+\\delta t}$, $P_t^{\\theta_t}(\\bold{n})$ is the joint probability at time $t$ with parameter $\\theta_t$, and $\\mathbb{T}$ is the chemical transformation matrix.\n",
    "\n",
    "The dimension of the matrix $\\mathbb{T}$ is $N^M\\times N^M$, where $N$ is the number of molecules, and $M$ is the number of species.\n",
    "\n",
    "The loss function is the Kullback-Leibler divergence between the joint probability at time $t+\\delta t$ and the joint probability at time $t$ transformed by the chemical transformation matrix. \n",
    "\n",
    "> Why not implement the time $t$ as an embedding?\n",
    "> \n",
    "> Tried , not promising due to the long step negative probability problem.\n",
    "\n",
    "> Is this size big enough? Given $N$ and $M$, the distribuiton is of dimension $N^M$, where $\\bold{n}$ indicates the number state of he system\n",
    "> \n",
    "> e.g. $\\ket{1,0,0}$ means there is one molecule of species 1, and no molecule of species 2 and 3.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## New directions\n",
    "\n",
    "This training use a neural network to replace the transition matrix via sampling from the matrix. Such method shows promising result in simulating chemical reaction.\n",
    "\n",
    "It is essentially based on a variational method with auto-differentiation. Another promising direction is using \"stochastic auto-differentiation\", which is developed in SciML. See talk [Agent-based modelling for prediction of chaotic systems](https://iclr.cc/virtual/2023/14704) for more information."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}