{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\"Open " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Diffusion Analysis with MRtrix: Part 2 - Constrained Spherical Deconvolution and Tissue estimation\n", "\n", "\n", "This notebook is the second in a series on Diffusion Analysis using MRtrix and builds upon the preprocessed data generated in **Diffusion Analysis with MRtrix: Part 1 - Preprocessing**. The data necessary for this example will be downloaded from OSF. \n", "\n", "\n", "####\n", "Author: Monika Doerig\n", "\n", "Citation:\n", "\n", "__Andy's Brain Book:__\n", "\n", "- This MRtrix example is based on the [Diffusion Analysis with MRtrix](https://andysbrainbook.readthedocs.io/en/latest/MRtrix/MRtrix_Introduction.html#) chapter from Andy’s Brain Book (Jahn, 2022. [doi:10.5281/zenodo.5879293](https://zenodo.org/records/5879294))\n", "\n", "__Original Data from OpenNeuro:__\n", "- Hannelore Aerts and Daniele Marinazzo (2018). BTC_preop. [OpenNeuro Dataset ds001226](https://openneuro.org/datasets/ds001226/versions/00001)\n", "\n", "__Preprocessed Diffusion Data from OSF:__\n", "- Dörig, M. (2024, November 19). Diffusion MRI Analysis with MRtrix: An Interactive Three-Part Series on Neurodesk. Retrieved from [OSF](https://osf.io/y2dq4/)\n", "\n", "__MRtrix3:__ \n", "- Tournier, J.-D.; Smith, R. E.; Raffelt, D.; Tabbara, R.; Dhollander, T.; Pietsch, M.; Christiaens, D.; Jeurissen, B.; Yeh, C.-H. & Connelly, A. MRtrix3: A fast, flexible and open software framework for medical image processing and visualisation. NeuroImage, 2019, 202, 116137. https://doi.org/10.1016/j.neuroimage.2019.116137\n", "- For more details: https://www.mrtrix.org/\n", "\n", "__DIPY:__\n", "- Garyfallidis, E., Brett, M., Amirbekian, B., Rokem, A., van der Walt, S., Descoteaux, M., Nimmo-Smith, I., & Dipy Contributors (2014). Dipy, a library for the analysis of diffusion MRI data. Frontiers in neuroinformatics, 8, 8. https://doi.org/10.3389/fninf.2014.00008" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setup Neurodesk" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "vscode": { "languageId": "plaintext" } }, "outputs": [], "source": [ "%%capture\n", "import os\n", "import sys\n", "IN_COLAB = 'google.colab' in sys.modules\n", "\n", "if IN_COLAB:\n", " os.environ[\"LD_PRELOAD\"] = \"\";\n", " os.environ[\"APPTAINER_BINDPATH\"] = \"/content,/tmp,/cvmfs\"\n", " os.environ[\"MPLCONFIGDIR\"] = \"/content/matplotlib-mpldir\"\n", " os.environ[\"LMOD_CMD\"] = \"/usr/share/lmod/lmod/libexec/lmod\"\n", "\n", " !curl -J -O https://raw.githubusercontent.com/NeuroDesk/neurocommand/main/googlecolab_setup.sh\n", " !chmod +x googlecolab_setup.sh\n", " !./googlecolab_setup.sh\n", "\n", " os.environ[\"MODULEPATH\"] = ':'.join(map(str, list(map(lambda x: os.path.join(os.path.abspath('/cvmfs/neurodesk.ardc.edu.au/neurodesk-modules/'), x),os.listdir('/cvmfs/neurodesk.ardc.edu.au/neurodesk-modules/')))))" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "vendor_id\t: GenuineIntel\n", "model name\t: Intel(R) Xeon(R) Gold 6126 CPU @ 2.60GHz\n" ] } ], "source": [ "# Output CPU information:\n", "!cat /proc/cpuinfo | grep 'vendor' | uniq\n", "!cat /proc/cpuinfo | grep 'model name' | uniq" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Import Python Modules" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "%%capture\n", "! pip install nibabel matplotlib dipy fury" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "import os\n", "import subprocess\n", "import nibabel as nib\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import tempfile\n", "from ipyniivue import AnyNiivue\n", "from IPython.display import display, Markdown\n", "from ipywidgets import RadioButtons, VBox, HBox, Dropdown, Output\n", "from dipy.data import default_sphere\n", "from dipy.core.gradients import gradient_table\n", "from dipy.io.image import load_nifti\n", "from dipy.reconst.csdeconv import auto_response_ssst\n", "from dipy.reconst.csdeconv import ConstrainedSphericalDeconvModel\n", "from dipy.viz import window, actor\n", "from dipy.sims.voxel import single_tensor_odf\n", "from dipy.direction import peaks_from_model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction" ] }, { "cell_type": "markdown", "metadata": { "vscode": { "languageId": "plaintext" } }, "source": [ "In this notebook, we focus on Constrained Spherical Deconvolution (CSD) to estimate Fiber Orientation Distributions (FODs) for the different tissue types. Additionally, FODs are normalized and tissue boundaries are derived from the anatomical image. The boundaries will be the starting point for generating streamlines in subsequent analyses. \n", "\n", "In the next part of this series, we will address registration and streamline fitting, advancing toward tractography." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Load packages" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['mrtrix3/3.0.4', 'fsl/6.0.7.4']" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import lmod\n", "await lmod.load('mrtrix3/3.0.4')\n", "await lmod.load('fsl/6.0.7.4')\n", "await lmod.list()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Downloading the preprocessed data from MRtrix Part 1" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "100%|██████████████████████████████████| 80.0k/80.0k [00:00<00:00, 73.8Mbytes/s]\n", "100%|█████████████████████████████████████| 226M/226M [00:02<00:00, 108Mbytes/s]\n", "100%|██████████████████████████████████| 8.31M/8.31M [00:00<00:00, 56.3Mbytes/s]\n" ] } ], "source": [ "! osf -p y2dq4 fetch preprocessed/mask.mif MRtrix_preprocessed/mask.mif\n", "! osf -p y2dq4 fetch preprocessed/sub-02_den_preproc_unbiased.mif MRtrix_preprocessed/sub-02_den_preproc_unbiased.mif\n", "! osf -p y2dq4 fetch preprocessed/sub-CON02_ses-preop_T1w.nii.gz MRtrix_preprocessed/sub-CON02_ses-preop_T1w.nii.gz\n", "! osf -p y2dq4 fetch preprocessed/sub-02_AP.bval MRtrix_preprocessed/sub-02_AP.bval\n", "! osf -p y2dq4 fetch preprocessed/sub-02_AP.bvec MRtrix_preprocessed/sub-02_AP.bvec" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'/home/jovyan/Git_repositories/MRtrix_preprocessed'" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Save the original directory\n", "original_dir = os.getcwd()\n", "# Navigate to the dowloaded data\n", "os.chdir(\"MRtrix_preprocessed\")\n", "os.getcwd()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Try typing one of the commands from the library, such as ```mrconvert```. If MRtrix has been installed correctly, you should see the help page printed by default when no arguments are passed to the command:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "MRtrix 3.0.4 mrconvert Mar 20 2024\n", "\n", " mrconvert: part of the MRtrix3 package\n", "\n", "SYNOPSIS\n", "\n", " Perform conversion between different file types and optionally extract a\n", " subset of the input image\n", "\n", "USAGE\n", "\n", " mrconvert [ options ] input output\n", "\n", " input the input image.\n", "\n", " output the output image.\n", "\n", "\n" ] } ], "source": [ "! mrconvert | head -n 18" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Constrained Spherical Deconvolution\n", "In order to determine the orientation of diffusion within each voxel, we will create a basis function from the subject’s own data. By extracting the diffusion signal from representative grey matter, white matter, and cerebrospinal fluid voxels, we will build a model to estimate what the signal should look like in different orientations and when we apply different b-values. The concept is similar to using a hemodynamic response function (HRF) as a basis function for fMRI data: We have a canonical shape of what we believe the fMRI signal should look like in response to a single event, and then we modulate it to fit the observed data.\n", "\n", "The response function is similar to the canonical HRF we use in fMRI studies. In this case, however, we’re estimating the response function for each tissue type. If you happened to collect your diffusion data with multiple b-values, then this approach in MRtrix is called multi-shell multi-tissue (MSMT).\n", "\n", "### dwi2response\n", "\n", "Unlike most fMRI studies which use a basis function that has been created beforehand, MRtrix will derive a basis function from the diffusion data; using an individual subject’s data is more precise and specific to that subject. The command ```dwi2response``` has several different algorithms that you can choose from, but for this tutorial we will use the “dhollander” algorithm:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "dwi2response: \u001b[03;32m\u001b[0m\n", "dwi2response: \u001b[03;32mNote that this script makes use of commands / algorithms that have relevant articles for citation. Please consult the help page (-help option) for more information.\u001b[0m\n", "dwi2response: \u001b[03;32m\u001b[0m\n", "dwi2response: \u001b[03;32mGenerated scratch directory: /home/jovyan/Git_repositories/MRtrix_preprocessed/dwi2response-tmp-8O8KLK/\u001b[0m\n", "dwi2response: \u001b[03;32mImporting DWI data (/home/jovyan/Git_repositories/MRtrix_preprocessed/sub-02_den_preproc_unbiased.mif)...\u001b[0m\n", "dwi2response: \u001b[03;32mChanging to scratch directory (/home/jovyan/Git_repositories/MRtrix_preprocessed/dwi2response-tmp-8O8KLK/)\u001b[0m\n", "dwi2response: \u001b[03;32mComputing brain mask (dwi2mask)...\u001b[0m\n", "dwi2response: \u001b[03;32m-------\u001b[0m\n", "dwi2response: \u001b[03;32m4 unique b-value(s) detected: 0,700,1200,2800 with 6,16,30,50 volumes\u001b[0m\n", "dwi2response: \u001b[03;32m-------\u001b[0m\n", "dwi2response: \u001b[03;32mPreparation:\u001b[0m\n", "dwi2response: \u001b[03;32m* Eroding brain mask by 3 pass(es)...\u001b[0m\n", "dwi2response: \u001b[03;32m [ mask: 94895 -> 66960 ]\u001b[0m\n", "dwi2response: \u001b[03;32m* Computing signal decay metric (SDM):\u001b[0m\n", "dwi2response: \u001b[03;32m * b=0...\u001b[0m\n", "dwi2response: \u001b[03;32m * b=700...\u001b[0m\n", "dwi2response: \u001b[03;32m * b=1200...\u001b[0m\n", "dwi2response: \u001b[03;32m * b=2800...\u001b[0m\n", "dwi2response: \u001b[03;32m* Removing erroneous voxels from mask and correcting SDM...\u001b[0m\n", "dwi2response: \u001b[03;32m [ mask: 66960 -> 66923 ]\u001b[0m\n", "dwi2response: \u001b[03;32m-------\u001b[0m\n", "dwi2response: \u001b[03;32mCrude segmentation:\u001b[0m\n", "dwi2response: \u001b[03;32m* Crude WM versus GM-CSF separation (at FA=0.2)...\u001b[0m\n", "dwi2response: \u001b[03;32m [ 66923 -> 36808 (WM) & 30115 (GM-CSF) ]\u001b[0m\n", "dwi2response: \u001b[03;32m* Crude GM versus CSF separation...\u001b[0m\n", "dwi2response: \u001b[03;32m [ 30115 -> 19977 (GM) & 10138 (CSF) ]\u001b[0m\n", "dwi2response: \u001b[03;32m-------\u001b[0m\n", "dwi2response: \u001b[03;32mRefined segmentation:\u001b[0m\n", "dwi2response: \u001b[03;32m* Refining WM...\u001b[0m\n", "dwi2response: \u001b[03;32m [ WM: 36808 -> 33730 ]\u001b[0m\n", "dwi2response: \u001b[03;32m* Refining GM...\u001b[0m\n", "dwi2response: \u001b[03;32m [ GM: 19977 -> 11830 ]\u001b[0m\n", "dwi2response: \u001b[03;32m* Refining CSF...\u001b[0m\n", "dwi2response: \u001b[03;32m [ CSF: 10138 -> 5102 ]\u001b[0m\n", "dwi2response: \u001b[03;32m-------\u001b[0m\n", "dwi2response: \u001b[03;32mFinal voxel selection and response function estimation:\u001b[0m\n", "dwi2response: \u001b[03;32m* CSF:\u001b[0m\n", "dwi2response: \u001b[03;32m * Selecting final voxels (10.0% of refined CSF)...\u001b[0m\n", "dwi2response: \u001b[03;32m [ CSF: 5102 -> 510 ]\u001b[0m\n", "dwi2response: \u001b[03;32m * Estimating response function...\u001b[0m\n", "dwi2response: \u001b[03;32m* GM:\u001b[0m\n", "dwi2response: \u001b[03;32m * Selecting final voxels (2.0% of refined GM)...\u001b[0m\n", "dwi2response: \u001b[03;32m [ GM: 11830 -> 237 ]\u001b[0m\n", "dwi2response: \u001b[03;32m * Estimating response function...\u001b[0m\n", "dwi2response: \u001b[03;32m* Single-fibre WM:\u001b[0m\n", "dwi2response: \u001b[03;32m * Selecting final voxels (0.5% of refined WM)...\u001b[0m\n", "dwi2response: \u001b[03;32m Selecting WM single-fibre voxels using built-in (Dhollander et al., 2019) algorithm\u001b[0m\n", "dwi2response: \u001b[03;32m [ WM: 33730 -> 169 (single-fibre) ]\u001b[0m\n", "dwi2response: \u001b[03;32m * Estimating response function...\u001b[0m\n", "dwi2response: \u001b[03;32m-------\u001b[0m\n", "dwi2response: \u001b[03;32mGenerating outputs...\u001b[0m\n", "dwi2response: \u001b[03;32m-------\u001b[0m\n", "dwi2response: \u001b[03;32mChanging back to original directory (/home/jovyan/Git_repositories/MRtrix_preprocessed)\u001b[0m\n", "dwi2response: \u001b[03;32mDeleting scratch directory (/home/jovyan/Git_repositories/MRtrix_preprocessed/dwi2response-tmp-8O8KLK/)\u001b[0m\n" ] } ], "source": [ "! dwi2response dhollander sub-02_den_preproc_unbiased.mif wm.txt gm.txt csf.txt -voxels voxels.mif" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### \n", "This command uses an algorithm to deconvolve the fiber orientation distributions (FODs) - in other words, it tries to decompose the diffusion signal into a set of smaller individual fiber orientations. You have several algorithms to choose from, but the most common are tournier and dhollander. The tournier algorithm is used for single-shell data and for a single tissue type (e.g., white matter). The dhollander algorithm can be used for either single- or multi-shell data, and for multiple tissue types. Estimating the FOD for each tissue type will later help us do anatomically constrained tractography.\n", "\n", "The next argument specifies your input data, and the resulting response functions for the different tissue types. The order matters; you can call the output files whatever you want, but it makes the most sense to label them as some kind of variation on the phrases “white matter”, “grey matter”, and “cerebrospinal fluid” (here, labeled as “wm.txt”, “gm.txt”, and “csf.txt”). The last option, “-voxels”, specifies an output dataset that shows which voxels from the image were used to construct the basis functions for each tissue type. This dataset can be viewed by typing the following:\n", "\n", "```javascript\n", "mrview sub-02_den_preproc_unbiased.mif -overlay.load voxels.mif\n", "```\n", "\n", "We will visualize the voxels used to construct a basis function for each tissue type using **Matplotlib**. CSF voxels should be colored red, gray matter voxels green, and white matter voxels blue:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAABcwAAAH/CAYAAABq7J4HAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAACfjElEQVR4nO3deZQeZ3Xg/yvZ2qWWWlJr3yXbknfjFa844AWMHWDMOkMwAUIgAySEhBAYbEOCCSQh/JwMGWAO9mQMBCeQ2D8wmHXA2NgGvK/a97UltVZbttW/P/Kzhq7nK/uq61X3293fzzk+B13qfavqqaeq3q6nnnsHdXZ2doYkSZIkSZIkSQPc4N7eAEmSJEmSJEmSmoEPzCVJkiRJkiRJCh+YS5IkSZIkSZIUET4wlyRJkiRJkiQpInxgLkmSJEmSJElSRPjAXJIkSZIkSZKkiPCBuSRJkiRJkiRJEeEDc0mSJEmSJEmSIsIH5pIkSZIkSZIkRYQPzKVetWLFihg0aFDccMMNh/zZG264IQYNGhQrVqxo2PYMGjQorrnmmsO6DkmSJEmSJKlZ+cBc6ob//t//ewwaNCjOPPPM3t6UlDvuuCNe+cpXxvTp02P48OExa9asuPzyy+OrX/1qb2+aJEmSJEmS1DR8YC51w0033RRz5syJe+65J5YsWdLt75k9e3bs3bs33vrWtzZw67q6+eab4/zzz4+NGzfGBz7wgbj++uvjv/yX/xLbtm2LL33pSy/42be+9a2xd+/emD179mHbPkmSJEmSJKlZHNnbGyD1NcuXL48777wzvvnNb8a73/3uuOmmm+Lqq6/u1ncNGjQohg8f3uAt7Oqaa66JY489Nn7xi1/E0KFDu/x/mzZtesHPHnHEEXHEEUcczs2TJEmSJEmSmoZvmEuH6KabborW1ta47LLL4sorr4ybbrqpWObqq6+OwYMHxw9/+MMu8d/7vd+LoUOHxgMPPBARnMP8wQcfjKuuuirmzZsXw4cPjylTpsTv/u7vRnt7e7e2d+nSpXH66acXD8sjIiZNmvSCnz1YDvPbbrstLrjgghgzZky0tLTE6aefXqR3ufvuu+PSSy+NsWPHxsiRI+OCCy6In//8593aB0mSJEmSJKkn+MBcOkQ33XRTvO51r4uhQ4fGm9/85li8eHHce++9XZb52Mc+FieffHK84x3viJ07d0ZExPe+97340pe+FB//+MfjpJNOOuj3f//7349ly5bF29/+9rj++uvjTW96U3z961+PV73qVdHZ2XnI2zt79uz44Q9/GGvWrDnkz5IbbrghLrvssti6dWt85CMfiU9/+tNx8sknx3e/+90Dy/zoRz+K888/P3bs2BFXX311fOpTn4rt27fHb/3Wb8U999zTkO2QJEmSJEmSGs2ULNIh+NWvfhWPP/54XH/99RERce6558aMGTPipptuitNPP/3AckOGDIn/9b/+V5x66qnxwQ9+MD772c/GO97xjjjttNPiz/7sz15wHe9973vjj//4j7vEzjrrrHjzm98cd9xxR5x33nmHtM0f/vCH4x3veEfMnz8/zjnnnDj33HPj4osvjrPPPjsGDz60MbOOjo54//vfH2eccUb85Cc/6ZJO5vmH+Z2dnfH7v//7ceGFF8Ztt90WgwYNioiId7/73XHcccfFxz72sbj99tsPab2SJEmSJElST/ANc+kQ3HTTTTF58uS48MILI+I/cpC/8Y1vjK9//evx3HPPdVn2+OOPj2uvvTa+/OUvxyWXXBJbtmyJG2+8MY488oXHqUaMGHHgfz/11FOxZcuWOOussyIi4te//vUhb/Pv/u7vxne/+9142cteFnfccUd88pOfjPPOOy+OOuqouPPOOw/pu77//e/Hzp0748/+7M+K3OvPPxi///77Y/HixfGWt7wl2tvbY8uWLbFly5bYvXt3vPzlL4+f/vSnsX///kPeD0mSJEmSJOlw84G5lPTcc8/F17/+9bjwwgtj+fLlsWTJkliyZEmceeaZsXHjxiJfeUTEn/zJn8RJJ50U99xzT1x99dVx7LHHvuh6tm7dGh/4wAdi8uTJMWLEiGhra4u5c+dGxH+84d0dl1xySXzve9+L7du3x09/+tP4gz/4g1i5cmW8+tWvftHCn79p6dKlEfEfgwEHs3jx4oiIeNvb3hZtbW1d/vvyl78cTz/9dLf3Q5IkSZIkSTqcTMkiJf3oRz+K9evXx9e//vX4+te/Xvz/N910U1x88cVdYsuWLTvwAPmhhx5KrecNb3hD3HnnnfEnf/IncfLJJ8fo0aNj//79cemll9Z+M3vkyJFx3nnnxXnnnRcTJ06Ma6+9Nm677bZ429veVut7f9Pz2/jZz342Tj75ZFxm9OjRDVufJEmSJEmS1Cg+MJeSbrrpppg0aVL8wz/8Q/H/ffOb34xvfetb8Y//+I8HUqrs378/rrrqqmhpaYk//MM/jE996lNx5ZVXxute97qDrmPbtm3xwx/+MK699tr4+Mc/fiD+/EP3RjrttNMiImL9+vXpz8yfPz8iIh5++OFYsGDBCy7T0tISr3jFK2pupSRJkiRJktRzfGAuJezduze++c1vxutf//q48sori/9/2rRp8bWvfS1uueWWeOMb3xgREX/7t38bd955Z9xyyy1x2WWXxU9+8pN4z3veE+eff35MnDgR13PEEUdExP8toPm8v/u7v+v2tv/whz+Ml7/85UX8O9/5TkREHHPMMenvuvjii2PMmDFx3XXXxaWXXloU/Rw0aFCceuqpMX/+/Pjrv/7reMtb3lK8Tb558+Zoa2vr5t5IkiRJkiRJh48PzKWEW265JXbu3BlXXHEF/v9nnXVWtLW1xU033RRvfOMb47HHHov/9t/+W1x11VVx+eWXR0TEDTfcECeffHK8973vjW984xv4PS0tLXH++efHZz7zmXjmmWdi+vTpcfvtt8fy5cu7ve2//du/HXPnzo3LL7885s+fH7t3744f/OAHceutt8bpp59+YPsyWlpa4nOf+1y8853vjNNPPz3e8pa3RGtrazzwwAOxZ8+euPHGG2Pw4MHx5S9/OV75ylfGcccdF29/+9tj+vTpsXbt2vjxj38cLS0tceutt3Z7fyRJkiRJkqTDxQfmUsJNN90Uw4cPj4suugj//8GDB8dll10WN910U7S3t8fb3va2mDhxYpc3w4866qi47rrr4gMf+EB84xvfiDe84Q34XV/96lfjfe97X/zDP/xDdHZ2xsUXXxy33XZbTJs2rVvb/uUvfzn+/d//Pb7xjW/EunXrorOzM+bNmxcf/ehH48Mf/nAceeShXQbe8Y53xKRJk+LTn/50fPKTn4whQ4bEwoUL44/+6I8OLPOyl70s7rrrrvjkJz8Zf//3fx+7du2KKVOmxJlnnhnvfve7u7UfkiRJkiRJ0uE2qLOa+0GSJEmSJEmSpAFocG9vgCRJkiRJkiRJzcAH5pIkSZIkSZIkhQ/MJUmSJEmSJEmKCB+YS5IkSZIkSZIUET4wlyRJkiRJkiQpInxgLkmSJEmSJElSRPjAXJIkSZIkSZKkiIg4MrvgoEGDDud2qEkMHlyOoXR2dqZikqQ8r6OSJEmSJDUf3zCXJEmSJEmSJCl8YC5JkiRJkiRJUkT4wFySJEmSJEmSpIg4hBzmGhj279/f25ugJkf1DMzFLEmSJEmSpP7AN8wlSZIkSZIkSQofmEuSJEmSJEmSFBE+MJckSZIkSZIkKSJ8YC5JkiRJkiRJUkRY9BNZ1FBSf+H1TJIkSZIkKc83zCVJkiRJkiRJCh+YS5IkSZIkSZIUET4wlyRJkiRJkiQpInxgLkmSJEmSJElSRFj0E1kQTzo4z4/+yeKgkiRJkiRJvmEuSZIkSZIkSVJE+MBckiRJkiRJkqSI8IG5JEmSJEmSJEkR4QNzSZIkSZIkSZIiwqKfad0tiNfMhfTqbBt9ljTLvkoDFZ2D2fNXkiRJkiRpoPENc0mSJEmSJEmSwgfmkiRJkiRJkiRFhA/MJUmSJEmSJEmKiD6Yw7zRebezn+1uLu5mzuGd3bbBg8txlZ5ocx1efe3Y9LXtlSRJkiRJUt/jG+aSJEmSJEmSJIUPzCVJkiRJkiRJiggfmEuSJEmSJEmSFBE+MJckSZIkSZIkKSL6YNFPKvJHxQDrfLa/FhKs7ivt+/79+1/0cwdbjoqD0nKZbYvov8ehN/SH9j3c25tto/7QlpIkSZIkSWK+YS5JkiRJkiRJUvjAXJIkSZIkSZKkiPCBuSRJkiRJkiRJEeEDc0mSJEmSJEmSIqIPFv0k2cJ82eKgfU22CGE1Rp+jwp11Ch9m29xCiodXncK43f2urDrnb6afN3o7SKMLD2c+J0mSJEmSpMbzDXNJkiRJkiRJksIH5pIkSZIkSZIkRYQPzCVJkiRJkiRJiggfmEuSJEmSJEmSFBH9pOhnnQKUdZZrFt0tJHjEEUcUsWeffbah20HrIM8991y311vV145fb8m2yeEulpstPrt///4idriPa/b7qZ/T9nZ3HXX6dH8tdixJkiRJknQ4+Ia5JEmSJEmSJEnhA3NJkiRJkiRJkiLCB+aSJEmSJEmSJEWED8wlSZIkSZIkSYqIflL0M1v8bqAXfqzuPxVWJHWWq9Pm3S102B8KuTZatk0y7VSniCRtR7bf1Olf2QKcGXX2v5H9sM53Nbp9JUmSJEmS+gvfMJckSZIkSZIkKXxgLkmSJEmSJElSRPjAXJIkSZIkSZKkiPCBuSRJkiRJkiRJEdGLRT+rBesaXVwuWxCvTgG/vibT5kOGDCli+/btK2JHHHHEi35/RMSzzz5bxA53Uc7eKvDZzH2JjhcdmwzaTyoYWafQJn0frZfW0cgCn9lCmHWuLZm2a3SBzzrLSZIkSZIk9We+YS5JkiRJkiRJUvjAXJIkSZIkSZKkiPCBuSRJkiRJkiRJEeEDc0mSJEmSJEmSIuIwFP3sbnG6bIG8RheMrFPo7nAXLj3caHuzBROpiORzzz1XxI48suxizzzzTGod3T2uPVFAlmQLVTayOGi2ECZtG8XouFbbM9tvsoVAs99XZ18zy1D/bfR1JLNttBx9f51Cq9nCpZIkSZIkSQONb5hLkiRJkiRJkhQ+MJckSZIkSZIkKSJ8YC5JkiRJkiRJUkT4wFySJEmSJEmSpIiIGNSZrPTWyEKF9H2NLL5Z9/v6mu4WTO1u0ceDrfPZZ58tYtkih9lihd1Vpz9Qm2TbjgpJDhkypMu/qTDqvn37ili28CNtBx0bWi99X7VIK20HfX/2vKTlhg4dWsSoTQhtS2Y7SPbcyi6XLYJb7SPZcytbVDVb9FSH10C6T0mSJEmS1Ff4hrkkSZIkSZIkSeEDc0mSJEmSJEmSIsIH5pIkSZIkSZIkRYQPzCVJkiRJkiRJiogmemA+aNCg1H91vm+gy7Tv4MGDi/+yx4b+279/f+o/0tnZWfyX0d3PRfxHwczqf0ceeWTxX519feaZZ7r8R8vs27ev+C+73N69e4v/aL+y7fTcc891+Y+2Y/To0cV/1JeoLUeNGlX8V13nc889Fy0tLcV/w4YNK/6jdVT3k44ftVEWtSW1E+0XtVP1c7S9je7nkqTmc8011zT0N+yKFSti0KBBccMNNxy2dUiS1Ex+8pOfxKBBg+InP/lJj6zPe7fUPzTNA3NJkiQpa+nSpfHud7875s2bF8OHD4+WlpY455xz4vOf/3zs3bu3tzevx916661xwQUXxKRJk2LkyJExb968eMMb3hDf/e53e3vTJEmKhx56KK688sqYPXt2DB8+PKZPnx4XXXRRXH/99T2+LV/96lfj7/7u74r4unXr4pprron777+/R7bDe7fUvHxgLkmSpD7l29/+dpxwwgnxjW98Iy6//PK4/vrr47rrrotZs2bFn/zJn8QHPvCB3t7EHvXXf/3XccUVV8SgQYPiIx/5SHzuc5+L//Sf/lMsXrw4vv71r7/gZz/2sY8NyAEGSVLPufPOO+O0006LBx54IN71rnfF3//938c73/nOGDx4cHz+858/rOs+//zzY+/evXH++ecfiL3QA/Nrr722Rx6Ye++WmtuRvb0BkiRJUtby5cvjTW96U8yePTt+9KMfxdSpUw/8f3/wB38QS5YsiW9/+9u119PZ2RlPPfVUjBgxovZ3HU7PPvtsfPKTn4yLLroobr/99uL/37Rp0wt+/vn0ZZIkHS5/+Zd/GWPHjo177703xo0b1+X/e7H7VF2DBw+O4cOHH9Z1HCrv3VLz8w1zSZIk9Rmf+cxnYteuXfE//+f/7PKw/HkLFizo8ob583+Uzp8/P4YNGxZz5syJP//zP4+nn366y+fmzJkTr371q+N73/tenHbaaTFixIj4H//jf0RExLJly+L1r399jB8/PkaOHBlnnXVW8VD++Ryp3/jGN+Iv//IvY8aMGTF8+PB4+ctfHkuWLOmy7M9+9rN4/etfH7NmzYphw4bFzJkz44/+6I+69bbYli1bYseOHXHOOefg/z9p0qQX/PzB8qD+7//9v+OMM86IkSNHRmtra5x//vnFH/W33XZbnHfeeTFq1KgYM2ZMXHbZZfHII48c8j5Ikvq3pUuXxnHHHVc8LI8o71Nf+cpX4rd+67di0qRJMWzYsDj22GPjC1/4QvG5/fv3xzXXXBPTpk2LkSNHxoUXXhiPPvpozJkzJ6666qoDy1VzmL/sZS+Lb3/727Fy5coD9aHmzJkTP/nJT+L000+PiIi3v/3tB/6/53OHe++WBpZeG5LKFKjLFjGg5Q6lAF5fkt3XwYPLsZDqcvQ5KsBJ30Voueeee+5FtyOC94tGTJ999tnUtmS+n7Z36NChRYy295lnnkmtY8iQIS/6WWrzkSNHFjFqS9oOWo72a/v27UWMij9mzqU9e/YUsTFjxhSxffv2pbajpaWliG3btq2IDRs2rIhlClhSu2WvGdRvsufNwQrcvpjscWl04ZeBdG2V1HfceuutMW/evDj77LNTy7/zne+MG2+8Ma688sr44z/+47j77rvjuuuui8ceeyy+9a1vdVn2iSeeiDe/+c3x7ne/O971rnfFMcccExs3boyzzz479uzZE+9///tjwoQJceONN8YVV1wR//Iv/xKvfe1ru3zHpz/96Rg8eHB86EMfio6OjvjMZz4T//k//+e4++67Dyxz8803x549e+I973lPTJgwIe655564/vrrY82aNXHzzTcfUntMmjQpRowYEbfeemu8733vi/Hjxx/S58m1114b11xzTZx99tnxiU98IoYOHRp33313/OhHP4qLL744IiL+6Z/+Kd72trfFJZdcEn/1V38Ve/bsiS984Qtx7rnnxn333Rdz5sypvR2SpP5h9uzZcdddd8XDDz8cxx9//Asu+4UvfCGOO+64uOKKK+LII4+MW2+9Nd773vfG/v374w/+4A8OLPeRj3wkPvOZz8Tll18el1xySTzwwANxySWXxFNPPfWC3//Rj340Ojo6Ys2aNfG5z30uIiJGjx4dixYtik984hPx8Y9/PH7v934vzjvvvIiIA783vHdLA0xnUkT0+H+DBg0q/qPlBg8eXPyX/Wxf+6+RbULLDBkypPhv2LBhqf+OPPLI4r/s9tJy9H2NbLcjjjii+G/EiBHFf8OHDy/+o7aj76PPVpcZOXJk8d/o0aOL/2jbsm3e2tpa/Jdtc9qv6n9Dhw4t/pswYULx35gxY4r/6Hi1tLQU/1GbZ49XtU/XuY7Q/tNy2e/L/Nfd43LEEUfgdtDxqnO98b/u/Sfp0HV0dHRGROdv//Zvp5a///77OyOi853vfGeX+Ic+9KHOiOj80Y9+dCA2e/bszojo/O53v9tl2T/8wz/sjIjOn/3sZwdiO3fu7Jw7d27nnDlzOp977rnOzs7Ozh//+MedEdG5aNGizqeffvrAsp///Oc7I6LzoYceOhDbs2dPsa3XXXdd56BBgzpXrlx5IHb11Venrhcf//jHOyOic9SoUZ2vfOUrO//yL/+y81e/+lWx3PLlyzsjovMrX/nKQdexePHizsGDB3e+9rWvPbBvz9u/f/+B/R83blznu971ri7//4YNGzrHjh1bxCVJA9vtt99+4O+Tl770pZ1/+qd/2vm9732vc9++fcWydI+85JJLOufNm3fg3xs2bOg88sgjO1/zmtd0We6aa67pjIjOt73tbQdiz9+ff/zjHx+IXXbZZZ2zZ88u1nPvvfcW98kX2i7v3VL/ZUoWSZIk9Qk7duyICJ7FRL7zne9ERMQHP/jBLvE//uM/jogo0qrMnTs3LrnkkuI7zjjjjDj33HMPxEaPHh2/93u/FytWrIhHH320y/Jvf/vbu8zsev4NtWXLlh2I/WZe9N27d8eWLVvi7LPPjs7OzrjvvvtS+/abrr322vjqV78ap5xySnzve9+Lj370o3HqqafGS17yknjssccO6bv+7d/+Lfbv3x8f//jHi9lSz888+v73vx/bt2+PN7/5zbFly5YD/x1xxBFx5plnxo9//OND3gdJUv910UUXxV133RVXXHFFPPDAA/GZz3wmLrnkkpg+fXrccsstXZb9zXtkR0dHbNmyJS644IJYtmxZdHR0RETED3/4w3j22Wfjve99b5fPvu997zts++C9WxpYfGAuSZKkPuH5lF07d+5MLb9y5coYPHhwLFiwoEt8ypQpMW7cuFi5cmWX+Ny5c/E7jjnmmCK+aNGiA///b5o1a1aXf7e2tkZE19Riq1atiquuuirGjx8fo0ePjra2trjgggsiIg48DDhUb37zm+NnP/tZbNu2LW6//fZ4y1veEvfdd19cfvnlLzo9/TctXbo0Bg8eHMcee+xBl1m8eHFERPzWb/1WtLW1dfnv9ttvP+wF3CRJfc/pp58e3/zmN2Pbtm1xzz33xEc+8pHYuXNnXHnllV0Gn3/+85/HK17xihg1alSMGzcu2tra4s///M8j4v/eI5+/91bv7+PHjz9w3200793SwNLnyurWyanb6Hy81e+r811ZtI7sfmW2L5tjuSfyy1O+8sx6KXc05UOvk2OacmdTfm5SfSuO8pzTvtP3U3tQvnL6LLVJJhc37TvdzOm76Ni3tbUVMSqcMnbs2NQ6snnzM8vU6SN1zsvqemmZbN707GdJT1zTqsybLumFtLS0xLRp0+Lhhx8+pM9lf7f85ttj3XWwWhrPX8uee+65uOiii2Lr1q3x4Q9/OBYuXBijRo2KtWvXxlVXXdXtehfPa2lpiYsuuiguuuiiGDJkSNx4441x9913H/ijvhGe38Z/+qd/iilTphT/P/3GkCQp4j/+Xj399NPj9NNPj6OPPjre/va3x8033xxXX311LF26NF7+8pfHwoUL42//9m9j5syZMXTo0PjOd74Tn/vc52rfI7vLe7c08HhGSJIkqc949atfHV/84hfjrrvuipe+9KUvuOzs2bNj//79sXjx4gNvhEdEbNy4MbZv3x6zZ89+0fXNnj07nnjiiSL++OOPH/j/D8VDDz0UTz75ZNx4443xO7/zOwfi3//+9w/pezJOO+20uPHGG2P9+vXpz8yfPz/2798fjz76aJx88skHXSbiP4qWveIVr2jEpkqSBqDTTjstIuLAferWW2+Np59+Om655ZYuM7aq6UKev/cuWbKky+yw9vb2LjO6DuZgA+kHi3vvlgYeU7JIkiSpz/jTP/3TGDVqVLzzne+MjRs3Fv//0qVL4/Of/3xERLzqVa+KiIi/+7u/67LM3/7t30ZExGWXXfai63vVq14V99xzT9x1110HYrt3744vfvGLMWfOnBec/kyefwP9N2fPdHZ2HtjmQ7Vnz54u2/abbrvttogITClzMK95zWti8ODB8YlPfKJ4Y+75bb7kkkuipaUlPvWpT8UzzzxTfMfmzZvT65Mk9X8//vGPcdbo87VGnr9P0T2yo6MjvvKVr3T53Mtf/vI48sgj4wtf+EKX+N///d+ntmfUqFGYRmXUqFEREbF9+/Yuce/d0sDjG+aSJEnqM+bPnx9f/epX441vfGMsWrQofud3fieOP/742LdvX9x5551x8803x1VXXRURESeddFK87W1viy9+8Yuxffv2uOCCC+Kee+6JG2+8MV7zmtfEhRde+KLr+7M/+7P42te+Fq985Svj/e9/f4wfPz5uvPHGWL58efzrv/4rpsN6IQsXLoz58+fHhz70oVi7dm20tLTEv/7rv6beiCN79uyJs88+O84666y49NJLY+bMmbF9+/b4t3/7t/jZz34Wr3nNa+KUU05Jf9+CBQviox/9aHzyk5+M8847L173utfFsGHD4t57741p06bFddddFy0tLfGFL3wh3vrWt8ZLXvKSeNOb3hRtbW2xatWq+Pa3vx3nnHNO+qGFJKn/e9/73hd79uyJ1772tbFw4cID9+x//ud/jjlz5sTb3/72iIi4+OKLY+jQoXH55ZfHu9/97ti1a1d86UtfikmTJnV543ry5MnxgQ98IP7mb/4mrrjiirj00kvjgQceiNtuuy0mTpz4oqnYTj311Pjnf/7n+OAHPxinn356jB49Oi6//PKYP39+jBs3Lv7xH/8xxowZE6NGjYozzzzTe7c0APnAXJIkSX3KFVdcEQ8++GB89rOfjX//93+PL3zhCzFs2LA48cQT42/+5m/iXe9614Flv/zlL8e8efPihhtuiG9961sxZcqU+MhHPhJXX311al2TJ0+OO++8Mz784Q/H9ddfH0899VSceOKJceutt6beUK8aMmRI3HrrrfH+978/rrvuuhg+fHi89rWvjf/6X/9rnHTSSYf8fePGjYsvfelL8e1vfzu+8pWvxIYNG+KII46IY445Jj772c/G+9///kP+zk984hMxd+7cuP766+OjH/1ojBw5Mk488cR461vfemCZt7zlLTFt2rT49Kc/HZ/97Gfj6aefjunTp8d555134MGHJEkREX/9138dN998c3znO9+JL37xi7Fv376YNWtWvPe9742PfexjMW7cuIj4j7eq/+Vf/iU+9rGPxYc+9KGYMmVKvOc974m2trb43d/93S7f+Vd/9VcxcuTI+NKXvhQ/+MEP4qUvfWncfvvtce6558bw4cNfcHve+973xv333x9f+cpX4nOf+1zMnj07Lr/88gP5wz/ykY/E7//+78ezzz4bX/nKV+Kqq67y3i0NMIM6k9XUssWSGqkn1tlfi9rROjJFSrNFA2k5KrZYR6adaDuyb3pRYc1sIUy6AT8/fes3UfGwmTNndvl3tQjowbZt9+7dRYyKbe7atauIbdmypdvrqBb5pAIhtA90HCZPnlzEaCrc81XPfxMVR21vby9ie/bsKWLVgqHUbtQedfp0tvBL5pzLXgtpndSns8VHu1s8uI46hVb7GouZSpIkSY2zffv2aG1tjb/4i7+Ij370o729OZL6MHOYS5IkSZIkqc/Yu3dvEXu+ZsnLXvaynt0YSf2OKVkkSZIkSZLUZ/zzP/9z3HDDDfGqV70qRo8eHXfccUd87Wtfi4svvjjOOeec3t48SX2cD8wlSZIkSZLUZ5x44olx5JFHxmc+85nYsWPHgUKgf/EXf9HbmyapH/CBuSRJkiRJkvqMl7zkJfGDH/ygtzdDUj/VNEU/+0OBz+w+9ESBz+5+jmLZAp/Z/aLlsuutooKG9F379u1LLUfFPGl7Z8+enYpRocpmsXz58iL2yCOPFLETTzyxy7+POOKIYpnp06c3bsMOwR133FHEnn766SJWPf5UpJRQIVDqI1RolfocoX6YKTaaPe+zRTR7o8BnHT1RPPlw62vbK0mSJEnSQGDRT0mSJEmSJEmSwpQskiT1Cz0xU0uSJHVVZ8aY925Jknpe5t7tG+aSJEmSJEmSJIUPzCVJkiRJkiRJioiaRT8bWXQt+111pq31RtHP3irq1t1toYKOpNEFArP9a+jQoV3+TcULn3nmmSJG+1X9roiICRMmFLGTTjqpiNVRbadly5YVy4wfP76IUXuMGzeuYdsVEXH33XcXsTPPPLNh3793794ilj02ddx3331d/r1z585imT179hQxOg60D1Tgk2LUX6koZ7WP1Lm2UPvSZ6nQaDNd0/qjZikALUmSus+ULJIk9S2Ze7c5zCVJktRU6CESDTAeeWT5U3bkyJFFjAYFaUDx2WefLWJPP/10EaNBUdo+GihtJrTNFKOB7GHDhhUxemmB0B8p1PbUzg7c9k0OwkuSpL7EB+aSJEmSJGnAyA7MUozQABB9NrveLBqYpUFi2j76bDPNmBhog2o0kD9kyJAilh2AzPaNRvcDqb8wh7kkSZIkSZIkSeEDc0mSJEmSJEmSIqJmSpbemKbR6EKgh1udYqZ12re7n83m+MwW/cweG/os5casThei3JbDhw9/0c9FcDFPKvrZaNX8npSzc9OmTUWM9pXyqo4YMSL12YkTJxaxefPmFbHueuihh4rY6NGjixhtW2traxGbNGlSt7fllFNO6fLvjo6OYpmNGzcWsW3bthWx3bt3FzHqq9njQMe/qrtFcSM4n2w2p26zTMUz76kkSZIkSeop5jCXJElSr6HBPhr8pjyeVPSTBiL37NlTxGhAMftCAGn2Ap+Etpli1KZ79+5NrYOOJQ00Z194oIFrGpBWc3GgW5Ik9SU+MJckSZIkSQNGdhCHBmZpBna20GZ25iR9X2ZmagQXEe2Jwo4OjNVDx4hmlNPxpUHn7ED0U089VcSyLxVI/Zk5zCVJkiRJkiRJCh+YS5IkSZIkSZIUEU2ekiU7XalOsck6xeSqn+3u55odTQcjddqc8pfSdKFqjkoqrEjT5k499dQiRgUoyQMPPFDEWlpailhbW1tqHdV9PeaYY4plKC8o5V+lKVrr1q0rYnQMqegn7UMGFcdcuHBharkdO3YUsQ0bNhQxKspJRUpHjRp10O183tixY1OxpUuXFjGankb5VKkf0vS57l6DslNBs9NDJUmSJEmS1OQPzCVJktR/0EDymDFjihgV+KRBYhqwpOVosLORL00cyvdlB09pwDObl5b0Vq5ayo1KMXoJYuTIkUWMXlrI5l+ltqLlJEmSNLD5wFySJEmSJPVLNChJs0FpsJY+S4OXNIhIRToplp3RnZ1Znx1IJdnB3zqzy2mmLLV9tp0pRutodMHVbCHV7PEl2cHp6kz8CG4Dmo2d3ebeGqCXeos5zCVJkiRJkiRJCh+YS5IkSZIkSZIUEU2UkqUn8ioe7kKgjS7m2ejv624bZ7cj2+Y0/Y3QFKLqZ2mZCy64IPX9WZRvlQpVDhs2LPXZDMrjSUU/d+7cWcQon+uUKVNS670FYlckPtfa2lrEaMpWdiobFeBsb28vYjSlrLv7QObPn1/EqJ9v2rSpiFHh1uHDhxexzBS97DS+OtMvm5nT+iRJkiRJUk9pmgfmkiRJvaVOAUjxYDgNJNPAIX22o6OjiFFhx2wRxzovIQwdOrSITZ8+vYjRYO/MmTNT61i1alUR27JlSxGjdqE2mDp1ahGjApo0ILt169Yitnnz5tR6s7ll6XjQ8aXBZ3qpgArHEsqHS21KLymoHq+xkiSpL/GBuSRJkiRJ6vNocIYGPimW/T5Cg7UUyw4U0XqzBUhpFjYNctJAKs0IzrZBdlsoRgPqNJt43LhxRYwKhtKA69NPP13EaNC+paWliNEAM83spn2j2emrV68uYtRfCG0L9atdu3YVMTq+1A/oWNLANm1LdvukZmcOc0mSJEmSJEmSwgfmkiRJkiRJkiRFxABLyVKnEOjh1uhto++rTqvJToupM32GpvLQtCeatkNTpqrTss4999xubxuhqVuUB3TJkiVFjPJ7Uk7NbJ7NKip6Se1G7Uu5Vg836r90TGlqF+3Dcccdl1pvpsAnFUulYz9p0qQitmDBgiK2cePGIkbTHiknamYKYaML70qSJEmSJIkNqAfmkiRJxNyKeZSfk3KK0oAdDYhmB4QptycNTmbzvtKA5dlnn13ERo8eXcQabdasWUWM2nnDhg1F7LTTTitiNOCbNX/+/NRy27dvL2IPPfRQarlscVBCg+8Uo5ys1E+zLzJYCLQer7GSJKkv8YG5JEmSJEk9jAbzHFzIyxb4pIKSNBOaBt8IfZYGemk5Or40+EuDpjQYetRRRxUxGuilbcmifaNtzs4apwFrGuSsI7t9FKM+RO03YcKE1LbQbO+2trYitmLFiiK2devWIkZtmn2BgGZV04xv6n903OjcosKitH30IgR9tj/ridngdWah17kfNXod2e/LvjiTYQ5zSZIkSZIkSZLCB+aSJEmSJEmSJEWEKVlqTTGofrbRU+qyUycOd6G/7FQyQtOZaNoTTb2hPJPHH398l3/v27evWIZiY8eOTW3H3XffXcROOumkIrZo0aIi9uijjxaxxx9/vIiNGDGia+CWYpHYd2m5DzR9itpo2rRp5Rf2gi1bthSxHTt2FLGiPYL7F03Toylq48ePf9Fte+KJJ4oYTXGk7W1paSlilPP29ttvL2LDhg0rYpkiuHR+UBvR1Dnq55IkSZIkSWID/oG5JEmSGA3q0UBnNp/rpk2bihgNrmfzyBIaYD3rrLNSyxHKsbl69eoiRnk8KY8sLUdt0NraWsSyx6OObB5UGrQ/77zzUuv46U9/WsQ6OjpS66V+Rcvt3r27iNEA9MSJE4sYDY5T7lYaWDf/tCRJUt/nA3NJkiRJknqYAyyMBsFoYDYbo4FAGpjNFoXMHjdajraPBnBPP/30IkaDeT2BZrHWKQqZRYOSNPBOBTRphjm1PW0zafSsfhqspYHy9vb2IrZ8+fLUctkBZjo/aH/p+2h2P+0bLUdtny0cW+e8bCa9tc09sd5Gr6POdbe7zGEuSZIkSZIkSVL4hnlD9dYoTXdzp2e/i0YcSZ3RWfrsCSecUMSqI5Pr1q0rlqGpyTTCTKPio0aNKmK0jtmzZxcxmpq8bNmyIladdr3jvHLkPDtNnGzevLmI0ah7I+3atSu1HXTsaao7vT2xbdu2IkZT4inX+549e7r8m6a00zqz5zSdNzT9/9577y1iNNpefRNj586dqe1odN2DvjhKL0mSJEmSVIdvmEuSJEmSJEmSFL5hLkmSVEt/maFBs5uoGCXlAKXZQjT7iNDMoKzjjz++iM2dOzf12RUrVqS2hfaDikdm8+ZSPtdTTjmliNHxmDlzZhGjwqJZtB80q2njxo1FjNqKvo/y8J5//vlFjHLV3nnnnal1ZGdI0qyuDRs2FDHKI0vHl45lnf4sSZKk5uADc0mSJEmSelh/GXDtLTSQRe1XTc8YwUUD6xQRzaJBzlNPPTW1XLOjwXNKvUrH6MEHHyxiHR0dRYwKQFIhUEqhOn369CJGg7APP/xwEZsyZUoqVgeljJ00aVIRo6KaTz75ZBFbtWpVEaP2q4POo61btxYxOubZFLn08gad09T/6pyrPaGZ7gHNtC3NwpQskiRJkiRJkiSFb5g3NRrhOdzfnx35I8OGDSti2VG+OXPmFLExY8YUseooKa0zW6ST0PbSqPOaNWuKGLVTZto6jfLSftHIKhWv3L59exHLugJi1T5B308xGvWlfaCpyzQVev369UWM+hJtS3W5qVOnFsvQcdiyZUsRo+NMb4BQoVmack9vU1TRufDUU08VMWoPOqezI8UDfURZkiRJkiQNPL5hLkmSJEmSJElS+Ia5JElSLX1xNgbNZKIChiNGjChiu3fvLmKUZ5RmQWVnsh199NFFbMGCBUUsa+nSpUVs9erVRYxmB9FMM8rJSrONqMgkzayidqGZgFTQtM6MRFpvthAofZb2d9u2bUWMZpK1tLQUsUsvvbSI0UzChx56qIjRrEHKtUpom6kIKx3zbLHRvnjdkCRJGih8YC5JkiRJUh+SHSyj5bKfzab76wk04EUDrjSgSZ+lfWt00TsadKaBTxrgo8FpSnlJg9M0KEloMJmKTFJ6SCoySe1H20LpTbNFWKlN6RgtW7asiC1evLiIZfsBDepmB7sbjdpl4cKFqRj1oZUrVxYxOkbUT7OD54S+j2KU1pXO8+yLGtn19oTs9SXbr+pcrxzIL5mSRZIkSZIkSZKk8A3zhsqOSGeXq/PZjDrfT6PM2SKPNMp8zDHHHHQ7f9OUKVO6/Hvt2rWpz2WNGjWqiGULS9JUXRqJro580vRval+aPk+jqCNuL0fd431liOzatauIbd68ucu/d+zYkdoOGvkmNJ2ZCnDSFHbqr7R91Tcl9t1cjnrve2UZozbP9hHa/wkTJhQxOv7VNqfzaOLEiUVs06ZNRYxG32m6Pe0DafSbN5IkSZIkSc3EN8wlSZIkSZIkSQrfMJckSerXKMcrzWSi5Z566qkiRrORaKYRzWahmTYvfelLi1ij0awfmr1D+0GzjbKz9CgnK820y+ampGOUlc3ZmZ1xRG1F/SU74yxr2rRpqRgVAl2xYkURy+ZopsKnNEuOZqLRbC+aOZctSipJkqTDywfmkiRJkiT1sN4q0FYnHWhvFZ/LDnzSABUN5mX3g9ZLg42UinPy5MlFjNJ40n48+OCDRSybonLOnDlFrK2trYh1dHQUMRpQpxil0KSB8o0bNxYxQm1A/SBb6JXahbaZBn+zx5fWQQUvs+08fPjwIkZ9lwrCUppP2jfqf4sWLUrFHn/88SJGxVVJtqgroeNGbU99g/pVNtUt9QPq4z2h0SlQTbOaY0oWSZIkSZIkSZJigL1h3uiinIdbnRH9DBrRo++n0TaavkqfpdG7E044IbV91cKHtC006kejiFQMsVoIMoKn1tLUc9p/Gm2kkdRqX6KRUGpzLPA5oizw+fTryhHYZyaV7UQj8bT/1WNIU5KpzVevXl3EqD3oTYns1H7af+qH1ZH09WevL5ahKdTTp08vYtu2bStiTz75ZBE766yzihg56aSTitjPf/7zLv+mNxUIvZWQffOBNMu1UJIkSZIkqaf4hrkkSZIkSZIkSTHA3jCXJEnqz2jWEuW1pJlSlCOScmfSbCGakXLaaacVsSlTphSx5cuXFzGaCUMzjY466qgi9l2Y9XTZrFlFjGZCEZqRRbK5ZefOnZv6PpLN90nHd8eOHd2OEZpxRG1FeVqztmzZUsQoH+727duL2DnnnFPEaJbj97///SJG5wLZs2dPEaM2yM4Ao+/ri4VAnaEmSZL6Oh+YS5IkSZLUw2gQjAa86sgOVtQp5kmfzQ7w1VkHLUcDvdSm2XbOFhKcPXt2EaMBw3HjxhUxSkVKg2+0XtqPtWvXFrENGzYUMUIDd2TixIlFbM2aNUWMUqjSvtF+UKFNWo4G2amtaACSlqMB+uzgL6UvpbbP9isq8LlixYrU9tEAPa03e64uXLgwFVuyZEkRo+KgNDhN25ItUkwxOpbU/2i9lF6YjhH1v0xq3t7UTNvSzEzJIkmSJEmSJElSDLA3zLOjKI0ehc98f3bEPDvi1sipkLROilEBxpkzZxYxmpZKo65UbLOKRsBpP2nUvq2trYhRIVAa1aVRSRpFpHaqLkejlDStm0Zgaf9pdHnp0qVFjKax04hztT3pOFN7TJ48uYjR9HfqD/TmQDYFAG1LdaSfvovekth2Y1ngc88ryjZvbW0tYlnU5scdd1yXf//iF78otw2Kj2aL9tJy1OaOPEuSJEmSpIHGN8wlSZIkSZIkSYoB9oa5JElSs6szSyxbXJByOtLMLvos5XSkGVo0++aRRx4pYpRvltY7evToIvbAAw8UsStOP72IfRtm8yyaOrWI0ewdWi+hHJbjx48vYjSzKCv7WTq+NJuPZq/RZ7N5c0eNGlXE6uQxpll+NAuPzo/vfOc7Rez4448vYhdddFER+/GPf1zEqJ8SytdL6Dyi85faoNk5Q02SJPV1PjCXJEmSJKkJNDK15sG+Lyub+pNki2pmB9UyqSMjeNCKtpm+j2KUyvElL3lJEZsyZUoRy8oWCKQijjSYR8U3aR2UajM7MEvpN7OFRWnfaL0zZswoYjSwTelXaUCYYiNHjixiNOBPg5yEto8G1Gkwns4FepmBBo5pcJXSoU74+YQituCDC1Lbl7VgQfl98+fPL2L33ntvEcu2H/WXbNpR+j5qe3pZgNK/0nrpOkQxB5ibmylZJEmSJEmSJEkK3zBH2VH9aiw72k6jV9kR+OxytC3V9dJ30SgzofagEeBZs2alvo8KE1KsWuSStoNGz2k0edmyZUWM9p+mJmeLfmb6DX2ORj1pmi6hUWhqSxqVrRbHjIjYtGnTi64z+5YALUdvRdBIP41yZ49NdR10ftB27L0oN+2bjuHKlSuL2PLly4vYy172siI2YULXkf+xY8cWy9BbBHScab9ouexIe6OLJ0uSJEmSJDUT3zCXJEmSJEmSJCl8w1ySJPWSRudp7S+ybUCzcSjfYjaPJ80+oVk7Z555ZhGjWUBU4JPWS7O7KM8jzXqhbab8lwH5SI866qgi9rOf/ayI0QwvmuGUnWlYnS0XwbOysrP+boHYFdAPsqhNe+KzNCOK2jk7K3EqFHV98sknixgd3wsvvLCI/frXvy5i2Xy9NIuL9m3MmDGp7+uLhUAlSZL6Eh+YS5IkSZLUw7LpNuto9EB0thBodhAx+300IEyDUdSm2aKfNIBGqRezBT5p3yh95pYtW1LrbW9vL2JHH310EaP9oHSRtH3bt28vYjNnzixiWCTxtvIYPXt6OThIx4gKmlbTVB7M1q1bi1hHR0dqvdTOtC2Eim9S+2VfUiDZl0voeFAfH3pleYyOmFsuRy8f0EsU2eKgtB9nnHFGEaOB6IceeqiIUZtm10uob9A6qO1HjRpVxOilB0p/S9cw4gtF7HC/fGVKFkmSJEmSJEmSwjfMUXdHJLKj49lRb0IjhzQqSrHqCGOdgn7jxo0rYvPnz099ltDoMbVJdfoqjerSFHVajkadaYprtp1oBDdTlJKWyb5BQSPT9H1UkJXQSHK1L1FbEuqD1Pep3eiz2en5NBJcXS+9dUHfn50eTn2ORpDnzJlTxDIo/cEPf/jDIpYpMnsost+XvX5JkiRJkiQ1O98wlyRJkiRJkiQpfMNckiT1EvPx5VHeSJrdQm1Ks6ey+WZf9rKXFTGaGZSVnbVFM42yM62oqOYVUPST0Hpp+2hbaHbYjh07UuulXJfZop9XpJbi2YGbN28uYjQbjGKNntVEJk+enNoWOm7ZXKtLliwpYpMmTSpiL3nJS4rY448/XsRWrVqV2j6asUb9ivaD+j2tI8trsSRJUlc+MJckSZIkqYcd7oJlh4IGNLPpRUmdYnuEBpmybZUtiEjbcuqpp6bWQWhQkgpUUqFIGhSngqHULpTuNFv085l/KdNUxrVliHReUX4ftemEO8rtm3FtLoUoDURnizPSMae0stlio5TilVKcDhs2rIhRnySUNpRitG9UjHLBggVFjAZhafCXzmkaYKZ9GzNmTOr7qKAuvUBw//33FzF6QSR7fcleS6jtKRUrtT2d0yRbCJQ00z2lJxzufTMliyRJkiRJkiRJ4Rvmh122QF6dQnrdnd5MUzwpRiPZxx13XGrbaJorjWrSvra0tBSx9vb2F10mO+2VlqPRaRoJpO/LFrSs7mudqeg0Wk1tSVOQSebtBiqESaPS2ana1G7Up2mqe3ZEsdqe9P30NgDtF20v9QeKzZo16wW382DovKRp6uvXry9i1EfoGFIhVDofuvum0UAb7ZYkSZIkSX2Tb5hLkiRJkiRJkhS+YS5JktQjsjMtaHYS5W+kWTCUG5Vm1dDMkFe+8pWpz9JsJJrxRagNKB9ktpgnzQLbu3dvaltoP7IzX6j9aN9oHTTLh2YSUR7P7LYsXry4iFH70WwoantajvpGnbybNCOKZj/R9mVnZVIOVWqXrIULFxYxKq764IMPFjFqP5rVSecW9ZfeKswqSZLUH/nAXJIkSZKkHpZNrZkd6KgzIJJNL5kdhM2ugwa8sjHalmzqTDJ27NhULIsKItKgaXbwlwbLli1bVsSokCUdS0y/+GoY3D8tlz520aJFRYwGjjteVhbL/OWYXxYx2l/6vmyxRyrgOnHixCKWtW7duiJG7VwHpSWlfk8vH7S2tqa+j84Z6hvUztQG2YKXc+fOTS1H5+AFF1xQxOicueeee4oYDU6TbCpSuu5ScVraN0qdSwP5dN1wUPzwMyWLJEmSJEmSJEnhG+ZNIztNO1sclEYdq8vR6BhNc6XvouVoGm61SGcEj3TSFFkqwlgdNaQp6jQamp26SsvRaG322GTewMhORafRYDr29Nlvw0j3a4pIrtgmFSSlfafjQKPNNN2aZIuD0vZVl8uODNMoOrUvjfhm0xN01/HHH1/EskU/afQ9+zZEdwsZO9otSZIkSZL6At8wlyRJkiRJkiQpfMNckiSp12RmhEXwDB2axUVopsnJJ5+c2pYHHnigiGVn6dCsLZoFRNtHM3eyM3yy37dly5YiRrO5aAZOdpZednYUzRCiHJuU6/IXv/hFEcvmKKWZetkZQbS/lFs2KzuDifoQodlUNKOR2plmOdIMSULFWqkY74oVK4oY9Q2a7ZYtAkyzGCVJkvTifGAuSZIkSVIPowGgbKrOnpAdoCLZbc6khIzgAUgaCKTBI4oRSnlIaSZXr15dxGiAatasWakYDaoRagNaL8VowJUGFqdPn17EqE8SGmA+7bTTitiGDRtS35dNj0mpMGnwcurUqan1PvLII0WMtpkG48eNG1fE6hSOnTBhQhGj/kfnaltbW2od1MepD9H5Ri90zJgxo4jRsdy5c2cRo8KYWTSYTMVBH3vssSJGbUrXMEpDnD0/6Hyjczpb6JX6X51rtkqmZJEkSZIkSZIkKXzDvFdkR9tppIpGoLJFP2mkKvO5bPFCGjHbsWNHEaMR9OXLlxcxmjZbHemkUezs9NNsMVP6vuzbIHS8qrJT1mlfKUajvFfAejuT0+L37NkDn+6K9iH7xge1G20HfV/2fKjGsm8rZGPZY7906dIidtRRRxWxDGrz7FRuKrSa7XN0HDLT6LtbLFSSJEmSJKkn+Ya5JEmSJEmSJEnhG+aSJEk9gmZQZHNx0myfbJ7CY445pohRTsxVq1YVMcoHSTOvaDZLdhZJdtYWzaKhWS80423t2rVFjHJnZotlZmeH0XLUppSzM5sTk3KPZveDYiQ7o49y5NKMKypySm1AeV/puFFbZbeZ2qqjo6OIZYt+kmOPPbaIZWdh0jlDeUup/agN6PzIcmaYJEkaKHxgLkmSJElSD8sOjPXWYAUNVNK2ZAdws4Oh2XSHFKN2ocHBOXPmFDEakKMBaxp4onUsXry4iC1YsKCIzZw5s4hRylIaGKOBNipuSQUgt2/fXsQ2b95cxDZt2lTEaHCfUk1SKkhq+yw6vu3t7UWMXgwgtL90zGnAlY4H9YNp06altoXQOiZPnlzEqA0y6WEjInbt2lXEaHCfCprOnz8/tQ7qpytWrChidDwo5eykSZNSy5FFixYVMTov77jjjiKWffGDlqNrHb30kE0dTGls6fuavRBo9j5IDvd90AfmfUydzlTNCU43enqjJptjmSpH082BTnY6iWlbqlWJ6QShG0P2rTRaLnuxy663etPLvqlH6EcOvU20bdu2IkZvJ9HNsvqDgdqDtjeTM/9gy1EfoYrU2Tfjqscmm5uc+mX2hkPLUX+47777itj48eO7/Jt+4NKPJ/rhSseUfkxQH6Hl6Adk5kZFy/hWmiRJkiRJajbmMJckSZIkSZIkKXxgLkmSJEmSJElSRJiSRZIkqUfyw9I6qunSDiabMoxSKWVzhVKO0rlz5xaxZcuWFTEqbkmpvyiXKaWIo7RTlPqM1kupzihNVnb7KCVYNo0X5QClz9bJxUnHiFKl0f5SH6L2o/2gdqHjRjl8KV0Y5UYltH2UxzObQoz6H+1vo5111llF7Pbbby9idN2gY0n9itQp+tkTqdQsLCpJkpqBD8wlSZIkSeph2cEAGkjoicGFbNFAkt1milGdomxhPfo+GlylQb81a9YUMRqkowFNGmijelNPPPFEEaMBr+y+0UAlFf4bO3Zsar3Z4qo0cPfoo48WsdbW1iJGNaxo0JRqc1EbTJw4sYhlUbscc8wxRYzaigqzUuFYetGAYlSolAbF6UWIzn8v+8umy8oXIWi9dL4dffTRRYyOZRZ9ll7UoDp2NABO9bdoHbNnzy5i9BIKXSOoDaiPZ88ZimUHxemFDnpZgK5X2eLIvaWZtqXKB+aHWZ0fMrRc9qaReYOLbkobN24sYnSRJtkbOL3ZQjdN+hFT/SGSLdyZfTuMLk60bXWO4fTp07v8m35w0Y1i7dq1RWzbthuK2Lp15QVxy5Yrixjt/549XytiEV2PK/XBUaPeWMToxwe9SUk/6KiPZN/mo+Wq66DP1Xmjjm6wdGOn85LWWz2XqDDoOeecU8ROOOGEInb//fcXMZItbtvd4p0W/ZQkSZIkSX2BOcwlSZIkSZIkSQofmEuSJEmSJEmSFBGmZJEkSeqRAp+UlolSP1FeRkpfRWmTzj333NT2UT5D2hbKxUmpryhFGsUorVU2TRalv8oeN2o/yv2Yza9LaN9aWlqKGO3H+PHjixi1C6VFo1ymlNps5cqVRSxbGJP6SyYFWgQfy2x+YkpRR8uRbB+i9IF0ru7YsaOIrV69uoiNHj26iFEOVXLGGWcUsV/84hdFjPoBFVKlYrKUfpD6bm8xXZskSWoGPjCXJEmSJKlJ1SkOWuf7Go0GrWjwlwbV6qxj8uTJRYxq/1ARQorRoBXVaqJBv927dxcx2l+qB0UDhnTMN2zYUMRocHDBggVFjAYCaZtpkI7angZSyapVq4pY9uUDOr40eEnfR/XKKEbFMtva2ooYDZRTnThqP2orGtg+6qijitjot5T7u273utT30aD40qVLixj1P3r5gF4CmDVrVhFbuHBhEaN+QG2fvYbRSwVZ1dpzEXxePvDAA0WMCv5S+9H1j/opDWx3t5bhwbbFAeuSD8x7QZ1CoNlK5XSzqr4lRtWk6/w4mTJlShGjGwbdcKmgZaYyOF1gKJatDJytRk7os/R2E928q6ggK8VWrXpfEdu8eXMR27273A7uc1S49U1d/k03ialTywrS9KOR0MV6yZIlRSx7bOhYV28wdCOhm1D27chshXn6QUHnQ3Ud9KONUJvTjyLa3uwPXFqO3lbLaPY/6iRJkiRJ0sBjDnNJkiRJkiRJksIH5pIkSZIkSZIkRYQpWSRJkmrJplqjopqUNonSQVG+1KlTpxYxSrlEaB2URiyb97BOIcFs/sZsEUzat2wKt2zBS/osbQuh/aAcuZRijXJxrl+/PrUttM2UZi3b9tkUeLQcpffKpCiL4Dy81O8p9SBtH30fpQqkc4HaKtsPCKXfq5NWkM5LSg/XTEU/JUmSmoEPzCVJkiRJ6uN6q+4LDR5lByppkJNitG80QDVhwoQiNm/evCL2xBNPFDEaHKTBZKqHlB2opEE6Wo72g4pMUnFGqkP02GOPFTEaeKd9o8F9Gviktsoey+yAP/W1FStWpL6P2rm1tbWI0QApFXCldqF6ZXPnzi1i27ZtS62DBj5X/j8ry3X8YbkOOpbUVlSAlAaxqa7Wjh07ihi1C61j0aJFRYwG8keNGlXEqK/R92VfAqDCrNOmTStidHypoCkdN9pm2haSfRkkex2iczC7LQOJD8x7QZ1Cd9kfCvQ2SnU5eiNp5syZqW0jdCISuliSTNVfupDTzZEuTtRu2Zs3fR+hN3aq+58pAnowdHGmHz5UbZpubpMmvaeItbV1LehJN7GsWyB2GcSoj1B/yMaqsjcD6iPUHyhGN3baNlqu2s/pOGdNmjSpiK1bV1ZMp+2g40A/KulHuSRJkiRJUl9kDnNJkiRJkiRJksIH5pIkSZIkSZIkRYQpWSRJkhouW+iQclhmU62deOKJ3dw6zkOZzYVI+0EppyiHJaV/yua5zRb9zOZvpDbtbqqvg6HUVrQf2VyrlAKL0tNlU49RP8gWmcyk/zuYbDtTzk5C7UJtQNtMMerj2eWmT59+0O3sjosuuqiI3XnnnUWM2oCuL2PGjCli2eORzY/dW3m0JUmSGsUH5pIkSZIk9UN1Bj+yaICPBuSyA1k0+Eb7QYON48ePL2I0UESFGKnuExXvyxbBJOPGjStiLS0tRYwGDKkOEdXMyqLjQYVFR4wYUcToeKxevbqIUV+jtsoOElPbU1tRwUtC66A+SWhb6Pu2bt1axKgfUKHXnTt3FrFnXlXWU3synixiVLTyuOOOS60je05T22/atKmIbd68uYh9//vfL2LDhg0rYvSyxamnnpr6LPW/Bx98sIhRfTdCbUp1AOk8+ulPf5ravjo1D+naSe1ChW174l7R1/SLB+bNfGCzhTuzb3bQhYtuBnTzq54odDGnm2HWhg0bili2kjL92KEbaXW/shV/M991sO3IViSmHyu03pUru1a1pptWHbTO+fPnFzHaXjpejXQFxOiHafY4ECrImvkBSzcSWme2IjwVX6WbPa232r9omSw69vQjJvvDkNqkTmFcSZIkSZKkZmIOc0mSJEmSJEmSwgfmkiRJkiRJkiRFRD9JySJJktRodYo9ZospZgssUuqjRqP1ZguB0mdHjhxZxLI5Xkm2EGijCxNmc+kSSn+W3Y9szlj6LOWmzKaTo3bJpmzLph2kfaM8wdl0YbRvlEIsG8umX5s9e3YRo1zJhApyktbW1tS2UKrFbL7ZbJ9stGZOqylJkgY2H5hLkiRJktQP9cQgRJ1B3ez2ZWsc0YDX+vXrixjVG5o1a1YRe+SRR4oYbfPu3buLGA36dXR0FLH29vbUOmhgjAanaRBswYIFRYwGrZYsWZL6PkK11bJFSbP1j7KD+zQISwPg1KbU1+hYUvvRgDXVo6J+kC2GSm1Ag5xUhJXWQQOzVJQ0W1eLithSO1OxUTpGJ5xwQmpb6Jg/+uijB93O30QDzFRkl843ais6F84777widscddxQx6lfZwW5aLxWnzb6QMND1iwfmA+lNBDop6KSlwpqjR4/u8m9qt+nTp3d72+gmTDc+upDRBZ4uFNUTO3vzrlNUNXtzJHRDqh6bNWvWFMvMmDEj9f11HO4Cn4T6A/1AzB4HulHSOqrfRzfsbF+iddJNkipjZ3/wVIty0r5Tu5100klFjH6wUP/NtFsE/+jNHgdJkiRJkqRmZw5zSZIkSZIkSZLCB+aSJEmSJEmSJEVEP0nJIkmS1GjZAoaUvihbbI+Wo5RGJ5988sE2s2HmzJlTxNauXVvEKI0apYejdFLZwoSNLjiYzdVIxzyb+o/2g3Jx0r5l85ZSLJvrklCeVtoPWkd236j9KJ9mNjVatkAltSm1C6VVo/SBixYtKmLjx48/2GZ2QQU+Kc/tsmXLiti5555bxE477bQi9stf/jK1LbTeUaNGFbFdu3alvm8gpcaUJEkDhw/MJUmSJElSF9k6UDRARYO/NGhVp/jchAkTili2uOCkSZOKGNXloUEmGiSmgTEagKR6QNmBT/q+ap2yCC78RwUgaR00YJgdwM0WZqVjTjEamCXU/7L9KvsiRLbgZfZFCGpnitF6aZtp4Jj6GvUNKlhLddZogHTatGlFjAaTqVYeHaP58+cXMWoDOt+WL19exOh6QN9H5zQVZr3vvvuKGL1wMnfu3CI2YsSIInbGGWcUsXvuuaeIUVtlr7HZ6xC1FaH+V2fwvNHf10g+MG8S1CHohkMnN130tm/f/qLroAsjncR10M3r6KOPTn2WihquWrWqy7+zla5pObro0HGg5bJveNExrF7YqvsUwQUjszfpw+2W5HJXQCxbpZou4NkfYJk33rL9IVvc9blvlT9QVly6oohlC2FW+zDdwKk96EfCvHnzihhdM+hNsmyb0024eu5bBFSSJEmSJPUF5jCXJEmSJEmSJCl8YC5JkiRJkiRJUkSYkkWSJCmNUndlCw5SaiJKfUSyxQWpSCfljSRjx44tYuvXry9ilMKK8m5SW2ULO2ZzTmZl027RsWxpaSlikydPLmIrVqwoYtxfvlXERo78zxDLpeMaM2ZMEaO0e1xY9FaIvbaIUbtkj2U2BV4d2XOL+hC1H8Wy5yCh/kzHiIpvUopCylmczYucVSdVHWmWfKQRzZ0vVZIkNQcfmEuSJEmSpC6yhR1p0IoGcbLFI2nAhgYvaRBsx44dRSw7QEp1jqg46J49e4pYdtCFBsayA17ZIn90jLLbR+ul40HtR9uXHTzPFt+kWkxUhy27vzSgmX2ZIbtvdMzpuNFAKr2QQP2UYjTwnh0wpD5J59vmzZuLGJ2rtC0LFy4sYjS4TwVrN2zYUMTIjBkzihgNilNxX+rjGzduLGJUgHTixImp9VJbZa9NdCypD1F/ps/WqftXRzMPWDf8gfnhHrEf6G8E0IWGLsrUJtVY9kZAF/OVK1cWMXqDjS54jax2TcvQBYZ0t90OFst+X7Xdqfoy/QhrlqKfWXQcqEI9vSVFPwyyP/zoGlHtJ/Q5+oFHxSzpptFxQXkMaf9pH+h8qB5r6g/0Y5Fu9ITeIqWin9k3ZEm17Rpd9HOg3wskSZIkSdLhYQ5zSZIkSZIkSZLCB+aSJEmSJEmSJEWEOcwlSZLSsjn/smmTKJ0XfV+2wCIV6aTcmZQ2asKECUWMiltSLknat2wKNpJtA0qTlc2vS+ug5SjHJmltbU0t98gj/6WI7d1bpiLbufOmIrZly+uK2L59/1LEhg9/fRGjYz5r1vuK2OTJpxaxrF/+8pdFrE4/yKbFo2NO5yC1AeVpnT9/fnYTU6ZMmVLE6Jyh3KOUMo3OVTo/KK1iNncwpaWjtuqLTOEmSZJejA/MJUmSJElSF9lBXao/lC20STVzKEYDTzTISTWTaABozZo1RWzevHlFbMGCBUVsyZIlRYzqTdHAHbVfe3t7EaM2JVTklLZlzpw5RYzagAaiaVuyBWGpDWg56hvZulM0OEjfR30yO8ierZ/ELzj8O8TK7xs58j8VMTqWtH3UD+i40fdRnTyKZYvTbt26tYjNmjWriJHHH3+8iFGBTzq+c+fOLWIzZ85MrXfcuHFFjAr+/vrXvy5idM48/PDDReylL31pEcvUX4vgPp4t0EuxbK227IsL/VXDH5gf7sYbSAcnewGmt1FItXBgW1tb6nN0E6GLIF1k6OZPN026mNMJW70oZH+cZd9Koxih5WgddDOrtme2cGWzuCK53Pbt24sYHS8q+kk3hDoFWTPHNVs9OvujNFtNPlPxmvoWvfGYfQOVin7SW6nUlrSv2R+kjTSQ7gWSJEmSJKnnmMNckiRJkiRJkqTwgbkkSZIkSZIkSRFhDnNJkqRasimdsqmEnn766dRylA4pmxeUigaSJ598MvVZymtJ+S87OjqKGKWmyubJzLYprYPSf1FOx2w6razjjjsutdz27R8uYqtXl2ntjjnmo0WspeW0Q9+wBqiTjit7jLJp9qgPZXMvP/roo+Vn/9/ys8d/9PgiRpYtW1bE6Fygc2b27NmpdVC70PlBbUDtl82bWydFWp1zWpIk6XDygbkkSZIkSQNEtuYTDSLWqb9Eg0LVOlsREWeddVbqs4sXLy5iNHBMA2NUWG/58uVFjGr/HH300UWsDhpAyw6C0cDi3r17ixjVl6LPUowG2rIDx9T2XHPq3yCWqxs2duxbi9iiRYuKGPWNnvD0039exFauXFnEtm0rz4Xdu79WxIYNe30RmzhxYhGjgWM6Hrt27Spi2QKpY8aMKWJ0Tm/ZsqWIUW27ESNGFLHJkyento9q+z3wwANFbPQPyxpq8z84v4hRH6frULaYcRYVSKWB92zBWjrm1PbUptlYfx3s7rcPzLM37WY+sHQzoBsYVQKnt0KqN85p06Z1e9tOPvnkIrZixYoiRhet7FtAdAyrF63szTt7wae3+rI/VrIXjur20cV01KhRRayvoZsn7Ve2wjrJFm6txqjNs8eUbpK0TuqHtF+0XPWPE2qjbIFPQj9EaB/oxw69WUqFULv7IyF7Hg2kG7UkSZIkSeo55jCXJEmSJEmSJCl8YC5JkiRJkiRJUkT045QskiRJPYFSEFEOQUrBRHkPs2nlVq9eXcQo32y2wCflbm1tbS1ilAuR9o3agLaP0rfRcpS+jVJdUYopStmUjT388MNF7KSTTipijTZu3LhUrLds2rSpiGX7M8Wy6c2oX9H3UT/N9itKbfjUK8rve7K1LIpL20f7RvmEFy5cWMRWrVpVxCi/6fz5Zf7VJUuWFLFsSsJsGr066dCoXSRJkppBv31g3h9y2WZzAFNec/qhX421tbWlvp/+kKA/OGbMmFHE6I/vbJ7wzB9O9KOfYrQP2T+4sn/U0T7Qeqv7kM2v3ltuqfz7iuTnsg8xqP8S+qM5+/CkeqyzxTrG/p+yyM++S8t1Uq7v7IOYbP+qWrp0aRGjBylU9Gf06LLQyfTp04vY1q1bX3Q7IvhcyvTh7ENB0h+u8ZIkSb0hUyvqYOhvOUK/1eg38/jx44vYtm3bitiGDRtS66X9oL+N6e8I+t1Mf89SwbxRPygHmZ57dfn3BRUwpNpPtH0tLS1FjGoJ0WepsCgtR3+H0O96ait6RkA6O/+tiB1xBBWdvbKIUcHVceNOSK23t1D/yxeO/VQRWbas/Bu1vb29iNE5SEVs6ZjTywz09zL1A7q+bNy4MbV9kyZNKmKzZ88uYoS+79FHHy1ia05dU8S2jiz/7qW2yj5DomtE1lFHHVXE6HhQIVV6tkLPFWm5bH24/lAvMsuULJIkSZIkSZIkhQ/MJUmSJEmSJEmKCB+YS5IkSZIkSZIUEf04h7kkSQMd5ZjrD/nkehO1KeWmzNbmIFRngaxfv76InXnmmanPUiFL6huUCzZbe4PyPNI6qP2ydTCy6HhQ/kZaL9VzoDye2fy//QXV0qAaGZSzmPp4tmZO9rhla4TQsczm9ty1a1cRoz5OuUxp+ygfKeUdpjpElAM5ew+o0wZZ3nskSVJfUuuBefZHWDP/wd7M20Z/eNEfn/RHIKn+8KXPrV27NrXO1tbWF/3+CN4H+lFOaPuqBQvmzZuX2o4VK1akvp9kC3zSclRgobocPSig4i1UdKInZIt8VlEBENpXkn3oQH9IZ/4QpGWwwM3rykvkrAWzihgVIaLCJnQuZf4gpXajBxC03LJly4rY+eefX8SeeuqpIkbFY2h7qT2z51fmu5rFQCpwIkmSBhYaeMr+bU+/D+m3IP1WpUKWVFiUBtVGjBiR2j76nUu/y+sUNI3fLkPZwTxqA/pbimK0LbS/1PZ03Gj7qK2yA9sUGznyDUWM2mDhwjOKmPgZCMUef/zxIrZ58/8sYrNnv6+IzZgxo4jRsyIqRkl/y9P1gPoaPQfZuXNnEZs/f37q+2ggn4oKU3+m/aDl6LpRZ4CZ0L7RczVq52zb037Q9aXO37rN/PyVmJJFkiRJkiRJkqTwgbkkSZIkSZIkSRHhA3NJkiRJkiRJkiLCop+SJPVbPVFXpK/loqsrm5Mwm/OP8pFS8T7KK0i5FSkv45o1aw66nb8JazlArsZG1iKJ4IKI1FZU/JDagHKyUs2LbJ0V2r7HHnusiC1atKiI0fmxefPmIjZ58uQi1luor/3qV78qYlS4k+prUNtTPtxs7R76PkL5cEk232f23KfcwfR92eKg9H3UBqNGjSpidC7QOZi9V2QLqdYx0O4pkiSpOdV6YJ798dLMP3Kaedvox3D2j1n6gVxdjn6A0x9xVCiBlsv+0ZotTpB52JApqhmR/+MlK/tHHW1fdVuyhRpPOeWUQ9nELjL9odFWrlxZxOjYUH/I/rFExzDzsIMeBtBxGDduXBEjtB11ip1Ui6fQsaLiJ/SHdvbBFu0/9Wna3szDhjoFgJtFnQcKzXyvkSRJfU+dYuT0u4xi2b+h6Dco/Y6cNGlSETvuuOOKGKG/Leg3PRUhpG3O/i1Eg2/0fbS/9LfPunXritjUqVOL2KxZs4oYtQH9nUeFE2lbqA/R36/ZQcRnnvkX+OxrixgVcM3+3aW8hQsXFrE9e64pYg8/vB6WKwfAqZ/SObh9+/Yili1YS9ch+jt1yZIlRWzChAlFbPz48UWMnmXRvtFzsBUrVhQxQgWJ62htbS1iGzZsKGL093V2gL4n9LW/iU3JIkmSJEmSJElS+MBckiRJkiRJkqSI8IG5JEmSJEmSJEkRYdFPSZIGlEbnjmv09/XFHPDZuha0HOVQpdyeVMyT8oxSXsY6xQWzRf6y+WvnzJlTxKhY4RNPPFHEKIcl5c6k/aVtphi1KaFj9NOf/rSIzZs3r4hRu1Du1oceeqiILViwIPVZyj1K7UI5QCn/JbUV5c3N5rnN1t8htA7Ki0z7QevN9pdsAU36LOUy3bVrV+qzdA7SvmVrBFEbZIsKZ69XdfTE9T6b+5s0+/1IkiQ1xoB6YN7X/gjP/jGRLdBSLVrQ0dFRLEPtQeuk78/+sUzL0Q/1zPfRHxv0hzf9oZL9gU9/vNQp/pdBhTKoyAu178yZM4tYI//orvMHN7U5HefsedndwqV0XGjft2zZklrn1q1bi1idwjrVPjxy5MhiGSpoRNeCMWPGFDE6DvTZOoWdq9/XH4p+1vkDW5Ik6XDL/o1C6G++7N9BhArm0eAWDexQcXtaLxUSpN+cWdnnBfQ7nIryZf8e2rRpU+qz1FbUBvT3Fa2X/kamNqC+wX3tNUXsiCPKv53oGNFAb3YwL7vN9AyD/o6jtqL9nTt3bhGjwf1mQn9XnnLKKUXswQcfLGL33/+JIjZt2u+n1kHtR88yKEb9nq4v1Iee+yYMMJ9Rxqgf1Pn7lfoaFeOlosKEzvPs9Sp7zc4+z2rmZ6iNZkoWSZIkSZIkSZLCB+aSJEmSJEmSJEWED8wlSZIkSZIkSYqIAZbDXJIkdV9P5HFv9rx41AaUY5PqHVBeRvo+yt1K+UgpfyPlR6TaHrTebJG/oUOHFjHK9zl58uQiRtu8Zs2aIkZtSnkjaT+onalf0X7Uqd1AOTY3bNhQxGjfqKgrFQylPLeUp5W+b/ny5UVs7NixRYzaOVunhnK3ZnMgU05g6i/Z+jhUm4T6OK2D2iqbN5fOGaonQttHOUqzeXjps9m2p75L/SB77je7Zr/PSJKk3jegHpj3tR9H9COXfjTTHw60XPUPNPrDhH700o9j+hFNMSpoSOiPFdqHiRMnvujn8Pt/VC63/fztRSy7D/SHBcXoGFaXo/2k71q/fn1quWwfoT8GTzrppCJWLXxJf/xnt4P+2Mz+wZ0tQkGx6nHIFpqhfaCiqvR99MCL+mvmHKE/gilG1zhqSypqQuc5Hdesaptki1D1teu0JElSd9X5LVTnN1Pm9/KhfJYGisi0adNS66XBwY6OjiJGv9VpWyhGRTpp+2iAigbpssVG6e/NJ554oojRADj9XqeigdkCrtm/aan9qF2yA9v09w8NEtdB+0F//1b/1o3ItxX9vUfH45hjjjnodjYDOh70XOD++99cxDZu/EYRGzWqXI76Af2dSoVo6RkVvXxAxUb3vbr8LG0L9QM6vtRW2aKz27dvL2LZop80yL5y5coiRtehbMFV+mz2vkD6w9/1pmSRJEmSJEmSJCl8YC5JkiRJkiRJUkT4wFySJEmSJEmSpIjwgbkkSZIkSZIkSRFxGIp+WsSte6jdKOk+FdygRPyZKvZUOIAKG9D3U/HCUaNGpbaDijtQjPoNFYGoogILG84oC1VSUQjS3WKeEbwP1eOaLXqZLQ5KfYmKSezatauIrV69uohNnjy5y783b95cLEMFKLOFb7KFO7P7Sm2eKbRKRXCpoA/1fepLtB1UJITWWy2iky10Qt9PxUzHjBlTxOoUcMqcI3UKSWX1xr1moN3fBtr+EjrP6Fyhc4/QtZLOA7pnUyxbcPnRRx8tYnSPnTt3bur76qDrIN0HqcAYXQup6Bgtl73GZX5TRfA203LUznQfXbt2bRGj40HFuKlN6V6VvTbTb5Xs9ZrWQd+Xva9RjH6bUHGtSZMmFTH6jTR16tQiRsXX58yZU8TomGeL29P9lPof/YajNsgW76O+S+hY1vnbryf+bmz2v017Yluyf2dmC0XSNtNy48ePL2IPPPBAEcsWmadrO10PaPvo+ke/h++7777UOuicmT17dhGjQqWrVq0qYnReUvvRfZ+Kg1LRVNoPuu7SPYpi1F/oGNHxoM/StZOuf7cUkYgrILZkyZIiRoUdaVuyRRKp39NvlRkzZhQxenaS3beeQH2S+l97e7kctRXJFoSl40G/m478TrktE39nYmpbtm3bltq+7PWKrrGzZs1KfZbQ9YrOD/r7g7Yl+4wr+1uvme6rjeQb5pIkSZIkSZIkhQ/MJUmSJEmSJEmKCB+YS5IkSZIkSZIUET4wlyRJkiRJkiQpIg5D0c/+muw9o06xOvosFc2gxP6UnJ+KNFTXMWHChGIZOn5UHINQgSkqJkDLUSERKhZRLcpBxVVo36nYUrYgW7YgAhWooH3InCNUmIH6AxWwoT5CMSo0QsUzqp+lfaLjR0VkCLVHthAoyRTnoKIyVAiMin7RZ9vb24vY1q1bixi1HR3XarEYOs7UR6gPVou2RvCxpz5N66Dj1d0+nS02R7L95nDfk5q9mFijNXp/s9/X6OWy20LonCXZ4pt0XahTFOg7cA+kglHHHXdct9dRR7bo8tixY4sY3aPouGULbWZ/j2SKSUfkC1RmCzLR/tL9gAqMUTtnf2tSjLYve4+gds6el7Qt9DuP9pdi9JuWCmnRbwKKZe+dtBwdS4pRcUE6P3bu3FnEqJ/SNSdbADJbAC1bKHKgFfjsLXX+RqXjm/2dTkUr6RpL53m2AB/1tTE/Ls/ptne1FbEVK1YUMToXqKgmXdvpWkfos9QGVMiX2qW1tbWIUZFJKiRNfyMOHz68iNF1ks4taj+6ZtPfJ7QthH7TUD+ge2i2CDVdd+m4UZ+ka+fy5cuL2LHHHlvEIllQsrfQfSZ73aW2zxZ/za5j98vLc/C54bnfEXT/pXOfCoXTMym61tVB20e/Baj96LkMPeeg/aVjRHrib7Xe0NxnpCRJkiRJkiRJPcQH5pIkSZIkSZIkhQ/MJUmSJEmSJEmKCB+YS5IkSZIkSZIUEYeh6GdWNdl7TyR6zxY8aeS2ZBPdE0qwTwUoqIACFc2golBVEydOTMVonZs2bSpiVHSA7Nq1K/XZaiERWueUKVNS35UtOFCnIBYVyql+lo4zFRWh76JjQ0VZqGgMFZWjYjDV9U6bNq1YZuPGjS/6uYh80Qjqq9lClRSrHms69lSsg75r8eLFRYz6A+1rtgBKtZgI9a1scVcq7kvoeNF5ky3+VN1/2oc652WjZa7Lda7n/UWdApp1Cr41erk6n6V9o+JB2cLdtF66HmVRMSxCBYComBPdZ7FQFbgFYpdAu2SLEGb7ULagMBXfpN9ZhPpBtnBYtigkxahf0W8GOr60LdkC6NnfPdltpuWoQFb29wD1A9rfbFE+6rt0btG+Pfzww0UsW/yLCotm7+P0u27dunWpz1Lfzd7f6HdBnWt2nd8BzVw4rNnQcSPZwrbZgskdHR1FjIpMZgv+0nk+9wNzi9jm3eXfQtTX6D49evTo1PbRb39qZ2o/Qtc/OldpW+jvwWOOOaaI0XHLFiinAoGLFi0qYtR+JHvdJcuWLSti1PbULnTM6b5Kx5fahfoptSkdtytgm6ldSLatslatWlXEuODva+DTuQLR1FYUy/7+o3VQ8de2trIIMH2WChdTbPXq1UWMfl/RtW7+/PlFjNC1M1uQna6T1P8yzwsj+m+BT+Ib5pIkSZIkSZIkhQ/MJUmSJEmSJEmKCB+YS5IkSZIkSZIUET4wlyRJkiRJkiQpInqx6GdvJHuvs85MwvpGF6TJFrwi2YKD3UVFEcaNG5dajgovjB07tohRobFqIYJssSoqjrl169YiRoU7ssXIsoWuqoW4soVVqCgL7Rf1VSpcQgVOMqZOnVrEqIgWFY2gGO0X7X+2YGim2CgVDKN2o+Ia1M+pT1OMzsHp06cXsWqb0HdR0Roq6JEt6kSoT1MhFlJtTzoXGq3RBR8zxWL7WuGSg6lTvLTZ26DRhVmzxZxItkAlXVOpaNa9995bxKhAIBXWpmsqFSOi6w/dP6kIEhWgousU7Ru1M13P6H5B9zcqMkz3RrpPZQuBUptmiy/R9tH9hrYve9/PFnLOFq6j/kzFRgltM/2mo2NJ5wy1aXY/qHB5TJ5cxsD9999fxLK/62h/6bNUuPPkk08uYvR7iu691PbZwo7UftT22d9A/fm+2uzo/KB+kC22R9dsuh7Qeum6NnLkyCJGfS1bJJvug7S/ra2tRWwyXA+2b99exGh/6Tyn84i+j9qAtjl736LzLfv7he4pxx9/fBGj+0wW/WagYp7btm1LfR/tL/U/2rdsv8/e4+mzS5YsKWJ0378C/vZuNPrdSffG7G8Guh5kn0dki3lmiw/T3/30O5bOQbp303HLPpOi6ws986IC4LS/VHiX7r/UfnQtycp+tj/cz33DXJIkSZIkSZKk8IG5JEmSJEmSJEkR4QNzSZIkSZIkSZIiohdzmB9uPZGPtbqOw50/91A+S7qbQ3jFihVFbM6cOUWMcp9SfqNs7mzK3VSVzQc7ZcqU1PdTTirKk5bN30V5r6p5zSjXG+W3ohyv1OZ0nCnXOaG82FWUf4y+f+fOnUWMjg3tP+XoyuYwz+SknTZtWuq7KOcZ5emlnGctLS1FjHIRUq61qsxxORyoLbM50KrLZY9fM+U6a6ZtkSRJkiRJPaPfPjCXJEk5jS6KXWc50t2B7IN9lpajGA0Q0UBStlghFQSj9f7qV78qYq961auKWFtbWxGjAkA0kJotOEj7sXLlyiJGg4mEBuFpcJIKtGWLQ9Up9JUtekcD4tkiotlidosWLSpihAZoqQ/R8c2eC9lCkdTHs8WhsudgduCW+gsdDyroR6iQFskW/aSC97QcFYHL9kmSLc5Ix4P6KRWLG0gFwfoq6gd1isDReUTnAhXbo3N1KhQ6zL5AQi/t0PlB90YqdEiFd+mlJTpnqJ3pZRlajs4Put7TtYSWI9mCjVnZF+PI3XffXcSyfTL7Uhbdp7OFvUm2wGd7e3sRo22mvksvZmV/D2XbhX5z0flL30fnAqFtzr5QRec+/e6kfp8tUJ69x1Mfohcbafvofp59qZV+t9Nxy17b6ZpIn61TyLc/MCWLJEmSJEmSJEnhA3NJkiRJkiRJkiLCB+aSJEmSJEmSJEVEP85hXiffaVYji34SyrNIebKyuSareag6OjqKZTZv3lzEKLcU5XiifGyUZ3LNmjVFjPJS0WerhTqz+WqzOTCp2CYV/aScXnS8KA9UNe8V5aCl/coWqqT9p3yphPa/eqzpOFOOLsrvRXlqaV8pH+amTZuKGOXjoli1f1F+OMqrls3xRucD5SOjHIMZlGOMcixSn6blTjnllG6vl64HmWtfnRy22RzF5kKVJEmSJEl19dsH5pIk6fBrpgGN7LaQbLHHbMHB7MAqDSRmiyDNmTOniD366KNFjAZcaRBy1qxZRezhhx8uYlSM7a677ipiL33pS1PruO+++1LrWLt2bRGjAko0iEmD07QcDdhSP6AYDbJSYSn6LBXNon5F/Y8GN2lbsgPDtI5GvmASwYW5aFtowD87AEsD19UXLw7mCohRscJscTzq9/SSAqHrBlm/fn0Ro2sJXZuoraiIGb24Qy8hOIDdXLIFJbPFfenlGeovVAiUrnV0LX7wwQdT25ctukj9nq51dM5QG9By2aKGdN8i9HLbjBkzihj9jiDLly8vYvQi0IQJE4pY9sWtLLpeUTvTPYDueXTvpusa9QM6Rtl7XnYd1E/pnpctUFmnwHad8yN7LSHZ30PZlzDppcbVq1cXMTovqV9Nnz69iNFLnXTM6aVTOn/pxVF6yXDBggVFjM7V7O9EWo6uu9l+kP2bqT8wJYskSZIkSZIkSeEDc0mSJEmSJEmSIsIH5pIkSZIkSZIkRYQ5zGvl2unuZ7M5sSg/F6G8R5SLq5p/iL6fcnxSDjn6fsrNddRRR6VihHIiPvbYY13+Tfkfs7kyKW8T5XSlPFWUc4xyOdG+Vtu4vb29WIZyRdFxJrS92RyblLdq+/btXf5NecXoc9n8Y3RsCOU8pONV3d6IMv8Y5Q+jXGaUE3DFihVFjPaVvi97HKo5iOlzGzZsKGJTpkwpYtl8dlQcNJu3l45/RvZaSP2mv+ZJkyRJkiRJvWvAPzCXJEndlx34qFM0sNEDJNnvo+Vo8JIKClExIhqEo4J+2YFEsmjRoiJGbT9x4sTU92WLatLAVrbwZLYgZ7YPUfE5GuinwU7aX/osDc5mi3U1uijujh07ihgNltKLEtT/xo4dW8RoEJi2mQpk0WdpvdSH6GUMKnxFfY3O1WzRz8cff7yI0QsJVGjuhBNOKGLUplnZFyaon9I2Uyxb/JDOS/o+9Qy6ltBLDPTiBZ0z9H0Uo/6SvWbTNZa+j64bhNZBfTdbzDP7ElC2ECidv/R9dF2jawkVH84W3qUihNkirHT/yF6b6LOrVq0qYtQutB/ZIp10P6JtyRY/pP3NFsnOnpe0HLUB9b/serPXe/r9kr1u0L5RjNo02wb0O4KOL92j6JjTS4zZIqd79uxJbR+1Ff1mpW2h85faj44bHXNqF/qNTscj+0Jef2BKFkmSJEmSJEmSwgfmkiRJkiRJkiRFhA/MJUmSJEmSJEmKCHOYp1GuomoOomxOSopRXqFsjPIKZfJLUV4lyml64oknFjEqEEj5kqpFOiMi5syZU8QoZxnFJk+e3OXfVPSRclJRoUbKi0m5KCl3FeX5mj59ehGjvKxV2YKRWWvXri1iU6dOLWLUdpRjs9pPsjnsqH0plkX7QNtLRWqrObooP9fmzZuLGBXWpONM+RipnSi2bNmyIlY9f7N56ig/8rx584oYyeYkzm5L9TrXTPmrJUmSJEmSDsYH5pIkqdsaPaCRLTqW1RMDLplB9QgerKOiOzSgTQWFqLAoDUzSwCwNsJGjjjqqiC1evLiIUeGwlStXFrFp06YVsWwb0AsC9Fl6aYCKf2Xbhdp57ty5RYwGcWkgno5ltpji/wsDmWdD/2ttbS1iq1evLmJUBJPagF4AoEHgLVu2FLHx48cXMSqaSqgf0MAtFT6lwfzly5cXMWp7OqdpOXoJhI7bE088UcToGjFhwoQitm7duiJG/W/+/PlF7NFHH02tN1vgjl4YocH17AtDarzsC1zZYpm0HL0ktX79+iJG1/vsPYpidO2k6y69LEIvCtH1iq5rdYqmZgsdUltRm9I2U9FAOuZ0n6bto/sWXWPp2k5tkH2xirY5uxz1SSqkSr/DqA9RW9G20P2I+gbdF+h3Cd0D6BpL/Znul3R8SZ1jlC2eO/S2su3jt8sQXYeo4CXdk2k/6Pyga0m2IDYd32yxYNoW+iyd59QG2W2pU2ye2mAgGdh7L0mSJEmSJEnS/88H5pIkSZIkSZIkhQ/MJUmSJEmSJEmKCHOYN1Q2D1A2hxDlb6KcRJRfivIZVfNLUS4yQrm+KLcZLUe5pR5//PHUeil3WDVGy1BbUk6qnTt3FjFqXzoOlOuM8klm1CnwScVXKUYFLZ9++uYiNmTIfypi1TamfHAnnHDCC27n4UJtN2vWrCK2atWqLv+mfkO576gvUSFQ6iOUc5TypdF5Uz03Kf8c7QMVJM3K5nDtbi7pbAFR0sjtkCRJkiRJeiE+MJckSSk9UXyzTmG4bMHQOoMwNABLhX1oMJsG16hAERVuogFCWq69vb2I0UBntvAkDSYS2hZCRcyoyBoNsmWPLw1E0iAj7S8NptPgXrYgIm0fFX2i4mm0jiuKSMTEiRNTsWOPPRY+Xcq2Cw3QZ/tptkBWdmCV+imtl7aZtiU7UE1FuKhILF03CLU9rXfr1q1FjK5rtG+E2pn6Lu1H9rPNrk5RtGbX6P2gfkB9l/opXZvofKNrIq2Dzrf77ruviFFRzalTpxYxutaNGzeuiFFx5Oz5lkXbQi/F0fWP7kf0u4ReoKN7VLY4N13H6TdItog3bR/tL/UN2l+6R1H/o89SMW26PtNvkGzxXOpX1H7UVhSjc59+WxA6vtlrO/2ue+rCXBFvimX3g/pV9qUwOpZ0zOm40cuJ2WOe7ffUptRf6NynNs3ez+n8oHVQrD8wJYskSZIkSZIkSeEDc0mSJEmSJEmSIsIH5pIkSZIkSZIkRYQ5zGup5j7K5r7L5pGjHGXZooGZ9VLON0JFDseMGVPEaP+3b99exLI5KmkfqrlZKW8T5XeinEr0/dSWlJ/ymGOOKWKHG+WyWrZsWRGjfFS8/2WbU66y6mfp+1esWFHE5syZU8R6AuUofOKJJ7r8m4q2Zot+UsFb6nPUp6m/0nlezZWYze1YR6ZQcATnFaZ8cdXvy14Ls7nx+kteUUmSJEmS1Fx8YC5JklJ6YqAiO2hSZ3CFPkuDQYQGL6koUFtbWxGjwcUpU6YUsU2bNqW2hQbOVq1aVcSoMBIV36QiVzSYSEWzaDlqUxrApcFEGrCj7Zs5c2YRo4Hd1atXF7ERt5cD+HsvLItD0SAroX5Fg7g94ZbkclRENFsYjpaj4nhUuI4GmilGLzJMnz69iNHxpZcA6PuyBTmpqBedC9R3qYgenZcTJkwoYvTyCfWrn/3sZ0WMBrSzxdOoXWhwnTRTAc1Gv9DUF2ULwmaLBlKbUmG9hx9+uIjNnz+/iNG9jO4fLS0tqfVSP6UXoI4++ugiRvdzahdqPzpXKUb3Lfo+utbR7wh6mSd7fLMvn2XvAVSgct26L8B635Lavmwhbrp20vU5W1Cc7mXUT+sUOKZjRPtL7UL3FGoXOpbZa132/pH9/Uf7kS3cmf2+7P2cXkbLnr/U9rR9dG4ReomVrmG0b/SbZvLkyUVs7dq1RYx+j9N6s4Was+cqaeb7rylZJEmSJEmSJEkKH5hLkiRJkiRJkhQRPjCXJEmSJEmSJCkizGF+2DU6hx/lKKP8XJn1Uo4xQrmHqPgo5aAaPnx4ah2UZ5JU94Hag9o3ux3ZgoaUL6qRaL8eeeRTRezppymXZJmR9Nln/7WIDRpU5q068shy/6s5qqiNqDAsHdOFCxcWsZ5QPYbUpyn/GuUooxx+1Pcpj1c2T291OcrPmM2Lmy3cSdtL51I2N2s1V2Q2hxlpplyoVRYplSRJkiSpf/GBuSRJ6rbeGtCosw76LBXTye7b1q1bixgVE6PBVirIRNtCy1FxHtpmitEAHg3+UVGg7AA5tRUNptJnaZtpcG7EiLJwJxU0pbbfe3E5GElFzMiYMWOKGLVVo2UHQamYZxYVoNq4cWMRo+NLhaVom7MF6aivUfE+ilHRtmxBU+qn1Cdpm6k/07ZQH6KihlRAmNDgOn2WiphRu2T7WrZQZG9ppm0hdC7UecmAPltncL9OgVR6GeXxxx8vYnQdnzNnThGjwo50X6CXwpYsWVLE/s//+T9FbO7cuUUsew+lopB0baJ7FBWXpr6RLc5NvxnoJTt6SYv2jYr8rVmzpojRtXPEiCth+3J9kraPrmvZAuDUN+j40nmU/T1EaL0UyxaxzV6fKUb9hVC7UBtQP6X+R+j40jZnCxfTctR3swU0s+2X/R2RLcyafWGTfm9ki3lm/+6h9dLLs7Rv2aLR2b+3euNFNVOySJIkSZIkSZIUPjCXJEmSJEmSJCkifGAuSZIkSZIkSVJEmMO8oRqdIy+byyeb566a84zyJWUL/2Xz7VEONMq1RPm/MnnHKKdSNndfNn8m5cZrJGrzhx8uC3w+99xlECtzaFFbDhlS5o7jnGtlrJq7i9qI2pxyoG7atKmITZo0Cbajsar9MJsnNZsfjvoS5dej3IaU32zy5Mld/p3Nqbt58+Yils35S8crm/Mt03Z18pU1S/7RbN40SZIkSZLUd/nAXJIkdVujBzSy39cbhV8ieNCMino9/PDDRYwK8VChIPLrX/+6iJ155plFjAbmqLggDXRRER8aJKPlskWaaNCcvo8G9qg4VLbgYLZQVXYgjwonzp49u4iRbGHM1tbWIkbFpuhY0gsLtL+ECutt2LChiFGxW1ovtT29LJAtKNne3l7Esi9e0LGkgWU6p2n7qJAWnYN0nmcL1k6fPr2IUQHD7LmaLSZGx4iOZW8Vfu4vsoXSSPY+mH0BhF6UoHXQS1d1XrCia/sjjzxSxGibafvo2knn6owZM1LLUbvQOU1FNadOnVrE6AUiurbTfYZetCF0TaR9o+VoHRTLXhOp/ehlIGo/KqBO7ULroHstXcPoHkDroGNO+0G/uWg5KvJMMTrf6EUnilER0WyB8jqFn7OfpeXoWNJvR+ov1M50PaDrBh1zug9mi3lSUdcs2j5qv+zLo3Se0zV74sSJRYz6C/226ImXh+t8trvfZ0oWSZIkSZIkSZLCB+aSJEmSJEmSJEWED8wlSZIkSZIkSYoIc5j3Odn8cJT7qZrnifJcLVu2rIhRLkbKNUd5gSiHF+WgyuayrOalyhY8zea8oxxhtL2N9OSTTxaxzs7LixjlhKM8XdRHKE8m5QCtFpuMKHNeZXPIUmzFihWpbaN+U0f1GFJb0nlEbUnnVvYczOYdzBb5rKLziIqPUj67cePGFTHKWUbHv7uaucAnaeZtkyRJkiRJjeEDc0mS1FDZInB1lquzLaTRAzg0AEsFirJFHKnAIpk5c2YR27hxYxHLDn7RoCYNiGaLAY4dO7aIURtQQSFqFyquSuuggXMaKKTB5GxhKRoUpeVoP2jglQpQTZkypYhlB4HrOO6444rYY489VsTomFPBUOovdNzoHKSXJWh/s8XJ6HjQZ6mgKW0fFU2l8436CxXwIsuXLy9iVLiO1jFp0qQiRn0oWzRV9WQL4dWRLYRM25K9XtE9j/pQthgbbTNtH/X7adOmFTEqnkvX2M2bNxcxOn/pGNH30fbRPW/Lli1FjPY3W3CVlqNjTtccuudRMUD6LF2bsi9WUYwKNlKxVkL3X7qOU9/IFk2lY0n9hfpV9t5NL2HR8aVtpmNExUvpN1f2xbxsQU7aDzqP6HpA+5ZdjvpVtuhxtqA9LUeyn81eJ+m40Qt0dYqw0rWJ2jSrt35HdHe9pmSRJEmSJEmSJCl8YC5JkiRJkiRJUkT4wFySJEmSJEmSpIgwh3lTozw7lCeLcj5m825VrV+/vohR/irKaUp50erkVcrI5neifaecXvTZRucBreZ2ozxitE46DtmcV/R9lMON8qFV83nR91PuMsrfRn36kUceKWJnnXVWEaujuv/ZPGPUHnVylFE/pNiaNWu6/JuOPR0Hyr1GlixZUsToeNH3ZYuZmudUkiRJkiT1RT4wlyRJDZUdMGn0ciRb9LNOcVAalKRBWRpEJjRIRoPSDzzwQBE76aSTihgVDKXChNu2bSti2QJjNOhGRcKoABrtG62DXgbo6Ojo9mepwCcNAFOspaWliFEbrF69uohRX6N2phgVcKWBzKlTpxYxGtykwXM6HnTc2traihjtL73sQW1K+zZmzJgiRttMy1G70Asf1IfonKYBbloH9T86F7IF3+jFAuobdM2hwlz0okadAl6qJ9v22cKY2UKb2ULXdV4momsJoXXQfZBesKKCjVT8mvaNCkrSvZHQ9YpePKPrS7ZIdva40fdRP6BrJ90rJk+eXMSyBQKzha7pGkbXxLVr16aWo7bK9j/6bZZ9IZF+N9Fvi7179xaxRhf3JXRfoPtWtk9SsfTs72dqv2wh0Oy9jIqr0u+NbNFZapdGH7fsC5CEfl9lCzrT7yE6FzIv3Ubkr1d19MQ6qkzJIkmSJEmSJElS+MBckiRJkiRJkqSIMCVLr6gzlYCmimTyTkeU01ZoegpNuaCpRjT1lfaBto2m8dE0EdrX6rQi+n5C30/bS8vRVBfaNtqvpUuXFrEtW7a86DqzU4FpOZrulGnLCJ4uRlPrMrI54ukYUo7tBQsWdGs7Isq+TvtZJx0D9ZFZs2YVsZEjRxYxOr+q6ROoLWk7ZsyYUcTo2K9cubKIZacy0nTOzPTf7DRfYj50SZIkSZLUU3zDXJIkSZIkSZKk8A1zSZL6rd4ojnIwjZ5RkC2QVWcd2ZkRNKsrW5gwW8iSipOtW7euiM2ZMye13jVr1hQxKmREsWyhQ5odQzNyqMgQzZihNqAZM9R+NMOHZsfQcrR9tA7qG9u3by9iVDCU2plmoJHsNtOsISoESjO/du3aldo+Oua0bzRrjgpQUVtRMTbqk3R8qSAYodl1tM3UBtm2p2J71J+pkCptHxWdpe3buXNnEcsWXlPvqXPvrnPPyxaQa3SBbbo20fUlWwCyzixX2t/58+cXMSqEvGzZsiJG9w/a32wBa7oHzJs3r4hR8ebsfYbQZ+naRNfEbNHPrOxM92x/oc/Sb5Ds76ZGy57T2YLi2XO/znpJ9npAbZrdFro3Ejq3sjP3aT+yxUuzBZgplm2XbNHeRn+20Xrjb1jfMJckSZIkSZIkKXxgLkmSJEmSJElSRJiS5bBr9HR4mhJBU7hoWmu1QGR2O2jqJk39yk7npPVOnz69iNE0oGqBRGqP7LQbKphJx4vWsXjx4iJG08wplpluRgVZCW0vFbQcN25cEaMpQXQMq1ORaFpdtmhrNr0BTWHftGlTEZs0aVIRI9SXqrLT52g56ku0vZQWgdI2VAucUhvRFDFCU0FpuiStg64t1Ddp2l21T9SZnthMaUUkSZIkSVL/5hvmkiRJkiRJkiSFb5hLktRvNdOb+I3elmxBnDqyxYholkm24NHkyZOLGM1QoZkhNLvnl7/8ZRGjApo0E41mlGXbgGYV0YwvKtqWnX1EBRGrM88iInbs2FHEaN+oTakoGhVIpfXOnDkztQ5qA/o+mtFF30dtT32S2jlbwHX9+vVFjLaZ+iStg2aIZWci0UwnaivqB1S8lGZ6Uf+j2ZV0LGmmHx2P7Iw9mlFG1xdqZ1ouW8RM/Ud2pl52uWwh0GyxuCxaL51HW7duTcUmTJhQxI4++ujUtjz++ONFjIoK0zWbChfTdY2uJXSu0oxmKgpO1wgqrDx+/PjUZwnd9+m4UYz2g65X1E+pXaivZYuRU4zuKdT/soUTSfa3bfZcrfP7gGLUJ7N66zd6tu0zM9IPZVuyxZGzn81mbqjzuz37fbR9PfF3WW/wDXNJkiRJkiRJksIH5pIkSZIkSZIkRYQpWQ677DSE7LQailERSZrWVZ1KlZ12Qt9PU0KyRT/nzZtXxGi61YYNG4pYtQgh7QNNLaPtoKlWtBxNOaH10rRXmvpbXQdNBaa2pD5CRSRpH7LHMFPMlT5H25adhk5T2WiKEE2Jp+KgdLyq66DziLaNYlQIk7aX2pzSLFC/qRbBzRY3JU8++WQRy/bp7LRt6hONnN6dnVpmcVBJkiRJklSXb5hLkiRJkiRJkhS+YS5JkvqJ7GyE7MyD7HI0Y4lmcrS2thax5cuXd3sdNJODYkuXLi1i2UJaNPskW7SNiqzRrCra5hUrVhQxmj1HxdOyM8iosBntL80YollU1F9of2l2VHaWD+0vFSKj2Uu0fVT8lfYtW2SS0Gwt2jc6RjSrq6Ojo4jRMafty86comNOM+La2tqKGBWOpTbdsmVLEaPZdLRv2SLAdQqC0fdZHLRvqnPPyxau64kZfbQtNCOZzg+6HlAh6VWrVqW2ha5XNMuZrrF0L6P7JS1HvyOyRaNpm+naTkUh6X5J9zK6dtJ1nI5HdgYs9T/6Popl7xV03LLFS6lNabnemgWb/Q2XPfezxSh7Aq23TqFhUmfmM/0GyRZMzhbKpfOXZNuqvxbzzPINc0mSJEmSJEmSwgfmkiRJkiRJkiRFhClZ+hya/kDTfmh6RnVqEU2jyk6RbW9vf6HNPGDq1KlFjKaTrFu3rojR9K2JEyd2+TdN86X2oOnK1EY0dYaWoyly2Smz1WM4+kfl1N3jPnJcEaMpbrt37y5iNKWPppXRtFpqz+o6stN3aHobtVt2ml6dwq2ZZWj6UnZaNZ03dI5kp1xVp1ZS/124cGERowKf1ObZ6VbZfaXitpnpanWmivfW90mSJEmSpP7NN8wlSZIkSZIkSQrfMJckSUnN9MZ+bxUdIzRjh2atUJE/KuI4bty4IkYzsSZMmFDEaJYRbR/N0MnKzubYu3dvar30fWPGjClimzZtKmI0i4ZmstFMK9oWKviWLXS4fv36IkazcrIFlLKz2bIxQjOOaH9HjhxZxKioHPUDmiFHRV3pnMkWqaP+QmhmHs30om2m5Xbu3JlaL52/dIxovdQ3sn2SZuxlr5P0WWp7NZeeKH5dpzhodh3Zc5/WQbM8qQBztiAnFdWk/aXrPX0fXevonKZ2oW2h60a2UHh2ObJ58+Zub0t2Zna2cGL22kTHjYo3U4zu8XQ8stvS6HO10er85u+Jfcv+FqVjTrLFQbPrzW4L3Wvpd2x2++gcJD3Rr5rp78bu8g1zSZIkSZIkSZLCB+aSJEmSJEmSJEWEKVkOu+yU8ezUjuzUBJqK0dLS0uXfVEiQ1knTsmiaCE2VPvroo4sYTYnbuHFjEVuwYMGLbl92alSm+GZEfkoMoWl9pDoNb9Z/nVUss29YefyyxTEnT55cxGiKPh0HmuJWnTKXLSBa7W8REfPmzStiDz/8cBGjad7ZaUiZKaG0DB1nSp9A0yhpuiWlVKDimPR91Vhra2uxDE0DX7JkSRGjY0P7T+1LRU9pCiW1XfX8yk5VI42+PmbXIUmSJEmSBh7fMJckSZIkSZIkKXzDXJIkJTXTm/jNtC0044FmBQ0dOrSI0Uwh+izNRpk4cWIRo1kgNKOsjmxxt2zhMPo+KppFM2GosCgtR/2FZsdkZ+UQWm+2yCkdX9o3KlqZnYFFaEYixWgmFc3ey8w4isjPdKpT9C5b0DTbN2iGXbafZouDNrpAVrZN6xQWVe/prcJ6da452f5HsrMVaaYqFdim6xBdX2jWbJ3zg7aFriXZWc40gzXbBnSfoZmtNAuZ7kfZawkdN9pfKjhdp8Dn9OnTixjNqKXCsRTL7ltfVKfgL/22pd852fXSNYf6JJ2r9H0dHR1FjPofHUvaN/o9ROj3ffZ3bDbWTJrpb7Xu8g1zSZIkSZIkSZLCB+aSJEmSJEmSJEWEKVmaWrbQHcVoCnB12tSYMWOKZWiKE8WyU2mXLVtWxE444YQiRlPEaPpqtQgjTWvJTs3LTgWj6T6EpufQFOEpU6Z0+Te1269+9asiRlOBaR8mTJiQ2rbx48cXMZpOtHjx4i7/zha43L59exGj9s1OS8z2uYxscVdqD9o2mgpJ66AiuNn+VUV9hND0N5o+SMc1ey5lrlXZ6xnJFvisMwWy0UVEJUmSJElS3+Qb5pIkSZIkSZIkhW+YS5Ik1ZItGpgtMknFQemzS5cuLWI0MyY7q6ZOIcbsjB/a30YXJqSZNbQOKliWLfCZLXzV6AKL2eUaXTySZszVQX2cjhGtl9q50cVQaVYfzd6k823jxo1FjLY5W2CW+jPNTqNZl9QG2dmqFKPj1ui+obzemgmXvYZlC0Rn15HdX+qT69atK2I0y3XUqFFFjM5Bml1KxSOpuODUqVOL2KxZs4oYzeqkoqTZWZe0HG0zFdqk9sv2A/oNQn2Drmv0fVTwvK2trYjRcaN7xSOPPJJaLvtbpbfUubZTW1Es269o5jL1KyrSSVkQsoU7qcg9nTPZ38DU/2i2ORXurDOTus7vYjWWb5hLkiRJkiRJkhQ+MJckSZIkSZIkKSJMyXLY9cQ0iex0y+oUFSoOmZ3mmp1isn79+iJG01yz09Sry2ULNdIUr+xUXdrX7LRDOv7VaWnZqay0D1QYtb29vYjR9MLW1tYiRvt68sknd/k3TUOqFmM9WIym5NEUOtovkk0LUG1jOlYUo2NDfS47jXzNmjVFbO7cuUWsitqciqpm0ynQFDs6b2gfaFu6Kzutr9EFPrMyhUslSZIkSVL/4hvmkiRJkiRJkiSFb5hLkqQBJju7oc73ZQttUqE+Ki7Y6MKdJFtQstGys4VoOZoxQ6it6LhRoSo6RjRbjmbqUNG27PGo0y7Z/pydcUXqFJWjz9Lssjr9uU5xQTrmtC20XLbts0XWaGYXbQv1XULFybLbTH28TiFfNb9GFxDOFkymc4H6Wvb7SPbcoqKB1C5UhDBbgHT69OkH3c4XQ/cj+s1A20LXDZrFSjOYaR3ZezJ9lmZTU6Hh+fPnpz5Lv6WoWHq2uHS2iG323l2nMHV2Wxo9Izf7W5Tu55S1gGak030mey2hNqXfINSv6HpAfYj6S53fZv2lwGej/7aqs95G9nvfMJckSZIkSZIkKXxgLkmSJEmSJElSRJiSpVdkpytkpxJkv686dSc75ZsKRtKUmuy0WSp82NLSUsRoili1cOm4ceOKZaZNm1bEVq1aVcSyU0iz+0XfR9PSqtN3aZoQHb/qvkdwG2UKvkZELF68uIgtWLCgiFX3P1sElqbG0WdpWtOcOXOKGE0hXr16dRGjPlzdPtpemv5F66Tldu3aVcTo2KxYsaKITZ48uYiNHDmyy7/vvvvuYhk6f6mvUr+kqZs0NY2myXW3UGedaVqNng7cEwVDJUmSJElS3+Qb5pIkSZIkSZIkhW+YS5KkAabRBT4JzbKh9VKxpNbW1tRy2YKSNNOk0TMtst+XnbWVnR1Cn60zk4/aOTNbLIJnTNF66XjQ7CiaqURFqbLrzaozs462OTsLi2ZwZWdTNbp4Wnbf6LN0PKhf0TGi40szPWm91dlpEdzHaftof2kd2YJldc7BRs8+a/ZCaVk9Uaw6q05R4TrqHHM6t+geSrMwx4wZk4qNHz++iGVnINP20SxkupbQjFBaB83Mnjp1ahGjIqf0m4baL1tIlbaP7oOzZs1KbQvdZ37xi18UsWwfou3L/lbJ3o+y19jsNtcptpy9D2Z/T9Ixohidb9THadY39TXavmyfpO0jtBxtM623TtHjbOHYOr+H6uite+3hLpDqG+aSJEmSJEmSJIUPzCVJkiRJkiRJighTsvSK7BSBOlNyMtNVaWoLFf6jGE2ZoikwND2FptdSUUqamlud/kbTyGiaNLUb7RdNnclOb6IYrWPTpk1d/k3tRvu1cePGIkZT9bLTand9tZy61vGBjiJWPTbUllOmTCliTz75ZBGj6cwbNmwoYtOnTy9iNCWZjg1Nk6r2OTpWVCy2vb29iO3cubOIUV+lY0/TpanvV6en0TLZ6fF0rtK2ZQv5ksz5kL2eZVMukOz3EYuDSpIkSZKkCN8wlyRJkiRJkiQpInzDXJIkNZFGFzbrCXWKENK+0Yya7IwtmilVZ7ZEneNRpxBPnaKBpM7MOJItdEgzq2iGH20LFTGjQlDZIlLZGXNUhIuKu9E6aN9om7PFPEmdIpM0w4qKYNL30bmV7S+E+iStg9pq7NixqXVki41S22eLwJE6179Ga/RssUZvc53tyxYNbPS1OHt9yW4foRmj2ZnAhPaD+vPWrVuLGO3HjBkzUstRQc4VK1YUsZaWliJGbXDssccWsXnz5hWxjo5y1jC1Aa2XrhtUCJTanq73NBOXZt1SwXPy+OOPF7Fly5alPpvtB1m9Nbu1TrHR7P2Sfr/Q7wNaL/0WoCwAdH+rU7yU7qt0XtJ+0G9v6s80K5v6EPV7aoM6xWTrtFW279aZ/d1bGnnf9w1zSZIkSZIkSZLCN8z7hexbA9WRFhoJo9FAGqmj3NnZ0TZ6k4feGJg0aVIRmzx5cpd/jxs3rliG9p1G82idFKPP0ggk7Re9yVQdkaPvpxF2Wo72lUZCaRRwzyvKt8X2HVXuQ/UYUq5v2s9qvvkI3l56S4reYqA3Jeg4UDtV39Cgfk7bRm+PZEeGacQ8O/pcfUuClsm+GUlvB9By1H8b+RZps+cSN1+5JEmSJEmK8A1zSZIkSZIkSZIiwgfmkiRJkiRJkiRFhClZJEmS0ihtEKVIovRN2eJpVCSM0itROihKuUSp1bJFfBpdLI7WQYXN6hR2pBilXKMUZoTWS6nIskUrCbVBtmBj9njQNlMbUEq1bOqxLDpnaPuoDRpd6JVkixo2ukAbpZrLpjek1Ih1+iR9tk6R00Zr9mLQWXUK9Na5FpOeuN6TRheP3LFjRxHL9l26V2zYsKGIbdy4MfV92XOQtpnSVra3txcxSuVJ1/E66QezBZ3pGNG2TJw4sYhRX3viiSeKGLVB9v5B19NsocNscWSSLdiY/Ww2hSVtX3Y5ki3YTeuoptY9lPVS4fFsAXW6X2Z/C0yZMqWIUdpgKtJO14jsdagnrvd11lun3zd6Pw53oVLfMJckSZIkSZIkKXzDvM9p5IhMdpSFRu9GjBiRiu3cubOI0Wh09g27Y4899qDb+UJoFJveFsi+5Zbdf3rDsPrGU/btKdo2Gh2lNqcirdk316qjstk38mjbsu1L7ZaVfbOwas2aNUWMjim95dnS0lLEaBSdYnR+LV68uMu/qT9kR6ipfbMFQ0mzFMfsiVF2SZIkSZI08PiGuSRJkiRJkiRJ4QNzSZIkSZIkSZIiwpQskiSpidRJt9NbReCy20JpqahAFqVI6ujoSK2DZIv4ZFMYUWqu7LZki0ORbNHFbCEyaudsYa5sEccs2pY67Uzpw+j7qJDW1q1bixgVzcqmaMsWjKJjVKd4LqFjRJ+tU+QqW1Qu+1laLxVeo23OpAaMyPc1+myd4q/KF5+rI3u9r1NkLXtNzK43mz4wu9465zTJ7ls29SZt3/bt24vYypUrixgVC6Z7QBZ9Xzat5oQJE4rY7NmzixgV6Vy6dGkRo2KKjS5kmVXn90H2HGz077U6hUWzxdIp9Sulw6XPUrFbWi6bEjV7j6J9o/5Cv33oNxL9lqJC3Nnjkf1dQvubTcdaJ21pNn0wbUujC9pn27S7++sb5pIkSZIkSZIkhW+YDyiZURoaeaEROFqORvloNJCKHNL30Yjetm3buvy7tbW1WGbJkiVFbOTIkUWsvb29iGVHQ2kkld6MqhZvjIgYO3Zsl39TkU4auaS2pBE1Wo6OPbUJvv343a77ten1m4plqI/Qd2XfuKSRYOojtBzFqscm29+yb3lSgVrqm7NmzSpi9KZItc9l34KkIqXU5nXedunu27913qzIyo4c99ZbyJIkSZIkqfn5hrkkSZIkSZIkSeEDc0mSJEmSJEmSIsKULJIkqYk0Oj1QT8imHKIiPpQiiVJL1SnYky18lS3Gli2cmC2QlS16V6fQHKFiSZSeK1vIqE7RT2qDOsecClBl+wH1yWyRtWyaNfo+Ku6W7WvZtsr20zpFEusUcczuR53iWtm+mz3feqsAc531NlPR6N5SZ38b/dnsfYY+my2AW6eYYvbaXqf4dfb3webNm4vY2rVri9i8efOKWJ1Ckdm0j5Tyk75v3bp1RYzSVNYpzljnntwT6txn6lzrsttCqE0pvSwdD/rNRSltqSAspcjNFgbO3vNoP+j6Quto9O9Y+r5sofU61yvaPjpudD2g3yWUirnRRT+zuvt9vmEuSZIkSZIkSVL4hnmfkx0lJNXR2DpvatGoIY3A0egTjSrR99Ho8V133dXl31REsVpUM4JHp2lEc8qUKanlaASNtpcKkFaXo1E7QqOKtM6nnnqq299Ho4D7Xtn12Ox9rlwmO5KZfYOJ2pe2d/fu3UWMis9W25hGkOn7qS3HjBlTxAh99vHHHy9iy5YtK2LVdqI+km03aqNGj4LXuS41cjvIQHtjTJIkSZIk1eMb5pIkSZIkSZIkhQ/MJUmSJEmSJEmKCFOySJKkHtDowkPNJFukiVI/UUqyOgUMKZZtU0rrRGnUSPb4UoqpbLEf2jdKuUayx6NOP210Ea4625ItRkltkC36Sd+XTVmX7ePZGK23TruQbFvVSWeWXUe2/w0fPryI0TmTLTqbLcTYaD1xzqj3UL+ne0W272avG3UKYpPseZktfk3rpTagApqUonT8+PGp7cte12i5bdu2FTH6nUPpMbPX4jqpJesUlu+JYsF1CoGSOvtLsttHhccJnfvZFK51igVnC5ln1fntXef6QucCbUu24CoV7aVjRKmY6ZjTcs3+O6LKN8wlSZIkSZIkSQrfMO8XGjkSWaegHxWMJDRylR3prL5BtHbt2mIZ2t7W1tYiRsUhR48eXcRolG7jxo1FrK2trYgdddRRRWz79u0v+v00Yk9vCdFoIcWyI7UTJkwoYpMmTery79WrVxfL7Nixo4jR217ZUdDsCCqtgwrNVvdh69atqe+n45AtIEv7tXLlyiKWOedoO2iddJyzx76ORhYCbfRbZJIkSZIkSYfCN8wlSZIkSZIkSQofmEuSJEmSJEmSFBGmZJEkqd/qiQJFWQOtyFo2VROlecoWNRw1alQRyxZaolRPFKNiP9ligNmihiRb0JTWQUWGaN9oP+i4UcGy7H7UKQBZp2BUtsAdHbesbPo00ugiVxQbOnRo6vuyGl0QjFID7t69u4hRX8sWEsxuC6H2o+3rCQPt/tETGn3NqVNYL1u0N4s+W+eanT0HSbZobxal6KR7FMWy5z59ltqF0oBSsVFK20rpMbPFULP7QXri/ptdL8n2U0LtQtfxbGHHp556qohlfyfSOU39JbteKhKbLU6b/T2ZLVpO/ZTaOZvaldZB5xEdt2zxTfptS9tM+9be3p5ark6B8p7Q3fX6hrkkSZIkSZIkSeEb5v1CdlS0GqvzdgAtR6NKNKJMI+M0wpV5Q4xGJbds2VLEpk2bloqRTZs2FTEqNkojn1RstFoclI4DvcGTfSMjO9JNn6X2rO4rjfpmR/9pNJK2g95soJHR6dOnF7HMSDX1LYrRd2XfsFi2bFkRI5lzlc4Z2jY637p7fYjIn/vNopm3TZIkSZIk9Q2+YS5JkiRJkiRJUvjAXJIkSZIkSZKkiDAliyRJ/ZZpanoGtTOlSKJYVraoVx3Zgo11NHqbG71eSs9VR/YcbHSBu0af+9m+Qen5ekuj+xoV68qi40FF75pJnf1V82v0NYdkC9HWuTdm9cT9LavOfSHbVr/61a9SsUZbtWrVYV9HnWLVpCfu01l1+il9lq7jvVW8uZlk2yqLCm3W0RPXROX4hrkkSZIkSZIkSeEb5gNKpuhntkBg5vsjeASYikaOHDmyiFEh0OqIKH0/jZrSaDcV5Ny8eXMRo2KeVICSlqN9rRZrpDango5UgDE7mkmjqLTedevWveg6du3aVSxDbUnHhmLZ/kUjrfQmG62D9ivzXVTMk0afly9fnvq+bGHNaoyWye57tn3rnPvZAq8ZdbZXkiRJkiSpLt8wlyRJkiRJkiQpfGAuSZIkSZIkSVJEmJJFkqR+wQKfkiT1Ld67JUlqTr5hLkmSJEmSJElSRAzqTA5rW2Ct52XbnA5hpoBfne/Pyq6DimhSrLotVOCTUAHCSZMmFbEpU6YUMSpyOWrUqBfdtoiIjo6OIlYtJPn0008Xy1ARRSrcSdu7adOm1GdpHUOGDCli1f0/6qijimWoPbZs2ZLaNlonFbSkfaDiqNRvqm1M6yRUtLW9vT31WdoHavPMPuzbt69YhmLZoqKEPnu4i372xH2lmd/cauZtkyRJkiRpoPINc0mSJEmSJEmSwgfmkiRJkiRJkiRFhA/MJUmSJEmSJEmKCB+YS5IkSZIkSZIUERFHvvgi6i3ZYp7ZwnmZ5RpdCJS+j4oBUuFLKvI4YsSILv+mAo/PPPNMEaPCijt37ixie/fuLWKtra1FbOLEiS+6bRERd911VxGbO3dul3/T9mYLS7a1tRUxakv6PioYSp9dt25dl3+PHz++WIb2YePGjUWM2ogMHTq0iFEBVSreSftQPf7UL6k/UOFS6kvUp6kQJu0XfV+1oGe2gCidM5kCwHV19/uy17g6eqK4sSRJkiRJ6j98w1ySJEmSJEmSpPCBuSRJkiRJkiRJEeEDc0mSJEmSJEmSIsIH5pIkSZIkSZIkRYRFP5va4S7ASZ+jQoJ1CvPVKaRHhQ6rBR2zRRSpmCcVqtyzZ09qO6jYKBUHJdOmTevy76eeeqpYhgpmVgtBRvBxoEKgmzZtKmKTJ08uYqtWrSpi1cKa9F3UHrt37y5i1EbU5kceWV6aqLDm8OHDi9iYMWOKWLWgJ7Xvjh07ili2wCehzxLqh9TnMt+f3TY6L+ucq40smFnnuxpdMFSSJEmSJA08vmEuSZIkSZIkSVL4wFySJEmSJEmSpIjwgbkkSZIkSZIkSRFhDvOmls0dTrFMjuJsvt/sctn8ySSbt7ia2/m5554rlhk5cmQRo5zYlCea8kLTck8++WQRo9zZlJ+7mj97woQJxTLVPOcREVu2bCli1fziB4tRvvJqXu+DraOaO33dunXFMrTvFKO85tW89AdD+fXHjh1bxCgX+cqVK7v8m44pfT/1r2xucsrrTnnoKYd55lyqk687e81oZnX2v6/tqyRJkiRJ6jm+YS5JkiRJkiRJUvjAXJIkSZIkSZKkiPCBuSRJkiRJkiRJEeEDc0mSJEmSJEmSIiJiUGey+lmdAmt6cdkifHWW6y76fiqQWKfoZ1Z1v2jbqMDn8OHDixhtLxVgpMKPI0aMSH2WYtVikPRdVDCT9nXmzJlFjIpobt68uYi1t7cXscxxpXYbN25cERs9enQRo4Kk1Ea0r0899VQqlimOSkU/s6gvEVoH7WvmXKX2oOKjmWK/h7KOZtYT18zD3SZ9rc0lSZIkSRoIfMNckiRJkiRJkqTwgbkkSZIkSZIkSRHhA3NJkiRJkiRJkiLCB+aSJEmSJEmSJEVEExX97I2Ca83ucLd5nfbNbluzHEMqNknFQWl7qXgj7T8Vw6RYtbDm0KFDi2WomCWhwpq0bfR9+/btS62jWlySiqBSAcpsMVPaXiqOuWPHjiJG7UvfV91XWob2IRujtqRtq1OUM7MdtF91Cpw2s2a5Z9TZjma5PkqSJEmSpP/LN8wlSZIkSZIkSQofmEuSJEmSJEmSFBE+MJckSZIkSZIkKSJ8YC5JkiRJkiRJUkQ0UdHPga5O4bhmKX5XLWYZkS/KeLi3l9ZJRT+pkCLtF30fFcOk4pVV2TaibcsWlqTlSOY8p++nNqLtzX5fdh9oucx6aXupPxAq8Fmnzbsr2+aZPqjeYdFPSZIkSZKaj2+YS5IkSZIkSZIUPjCXJEmSJEmSJCkifGAuSZIkSZIkSVJE+MBckiRJkiRJkqSIiMhVudNhV6f4W7MUjmuW7SDZoozZApxUIHLIkCHd2hZaJxUQpVi2GG8jC61mv4uKTWa3lwpaZouNZgqmUnHMbNHWbJFS0sjj0MznmyRJkiRJUl/lG+aSJEmSJEmSJIUPzCVJkiRJkiRJiggfmEuSJEmSJEmSFBE+MJckSZIkSZIkKSIs+tlvZYoLZgswkmzhQyrKSEUuewNtB+0DtRMViHz22WeLGO1/teAkFRDNFqDMtm92H0im39Q5ztm+NHTo0NR6qaBn9djs27fvRZeJaGyRzrqflSRJkiRJ0uHnG+aSJEmSJEmSJIUPzCVJkiRJkiRJiggfmEuSJEmSJEmSFBE+MJckSZIkSZIkKSIiBnUmq9DVKRCpvi1b+DDbR/pD4cM6bVKNUeHKIUOGpL6LCoZScVCSLQ6aOV51imNSkU6SLdJKxUyrn23kvveWbPHVbHFX9bxm7l+SJEmSJA1UvmEuSZIkSZIkSVL4wFySJEmSJEmSpIjwgbkkSZIkSZIkSRHhA3NJkiRJkiRJkiLCop8DSvUYZotU1inwSfpDobtGFoikQo1UlJKWyxZ+zMYyx5X2k76LlssW7nz22WdTn82qrqOvFfgkA/0c7A88DpIkSZIkNR/fMJckSZIkSZIkKXxgLkmSJEmSJElSRPjAXJIkSZIkSZKkiPCBuSRJkiRJkiRJEWHRz35R/O9wozaiWLYoY6O/b6CoU5C1TuFWKujZXbQd2dhAPlezx8p261s8DpIkSZIkNR/fMJckSZIkSZIkKXxgLkmSJEmSJElSRPjAXJIkSZIkSZKkiPCBuSRJkiRJkiRJERFxZG9vQG+z6FpXdYoL1vk+dZVttzoFHbPLZYqvZrejzvnmudqV7SFJkiRJktR4vmEuSZIkSZIkSVL4wFySJEmSJEmSpIjwgbkkSZIkSZIkSRFhDnNVmIu6OdTJEV8nr3lmW3qij6irOm3Z6LoEkiRJkiRJ/ZlvmEuSJEmSJEmSFD4wlyRJkiRJkiQpInxgLkmSJEmSJElSRPjAXJIkSZIkSZKkiLDopyqyBR0Pd2HJg6mud6AXKuyN/R/obd7XNPpclSRJkiRJ6s98w1ySJEmSJEmSpPCBuSRJkiRJkiRJEeEDc0mSJEmSJEmSIsIH5pIkSZIkSZIkRYRFP5VARQOz6hQXtFhhc/K4NC86NkcccUQR279/f09sjiRJkiRJUp/jG+aSJEmSJEmSJIUPzCVJkiRJkiRJiggfmEuSJEmSJEmSFBE+MJckSZIkSZIkKSIs+jmgUYHAbEFHimU/W6eIaPWzdbYj8/2H8tmBwvbI6Y2+NHhwOQbq8ZIkSZIkScrzDXNJkiRJkiRJksIH5pIkSZIkSZIkRYQPzCVJkiRJkiRJiggfmEuSJEmSJEmSFBGHUPTTwnGSJEmSJEmSpP7MN8wlSZIkSZIkSQofmEuSJEmSJEmSFBE+MJckSZIkSZIkKSJ8YC5JkiRJkiRJUkT4wFySJEmSJEmSpIjwgbkkSZIkSZIkSRHhA3NJkiRJkiRJkiLCB+aSJEmSJEmSJEWED8wlSZIkSZIkSYqIiP8P8Abm4y7EyGEAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Function to convert .mif files to temporary .nii.gz files\n", "def convert_mif_to_nii(mif_file):\n", " with tempfile.NamedTemporaryFile(suffix=\".nii.gz\", delete=False) as temp_file:\n", " result = subprocess.run([\"mrconvert\", mif_file, temp_file.name, \"-force\"], capture_output=True, text=True)\n", " if result.returncode != 0:\n", " print(f\"Error in mrconvert: {result.stderr}\")\n", " exit()\n", " return temp_file.name\n", "\n", "mif_file_path1 = 'sub-02_den_preproc_unbiased.mif'\n", "mif_file_path2 = 'voxels.mif'\n", "\n", "# Convert both .mif files to temporary .nii.gz files\n", "nii_file1 = convert_mif_to_nii(mif_file_path1)\n", "nii_file2 = convert_mif_to_nii(mif_file_path2)\n", "\n", "# Load the converted images\n", "data1 = nib.load(nii_file1).get_fdata()\n", "data2 = nib.load(nii_file2).get_fdata()\n", "\n", "# Choose the middle slices \n", "slice_idx1 = data1.shape[2] // 2 # Axial slice \n", "slice_idx2 = data1.shape[1] // 2 # Coronal slice \n", "slice_idx3 = data1.shape[0] // 2 # Sagittal slice\n", "\n", "# Prepare RGB overlay for tissue types \n", "rgb_overlay = np.zeros(data2.shape[:3] + (3,), dtype=np.float32)\n", "rgb_overlay[..., 0] = data2[..., 0] # Red for CSF\n", "rgb_overlay[..., 1] = data2[..., 1] # Green for GM\n", "rgb_overlay[..., 2] = data2[..., 2] # Blue for WM\n", "\n", "# Plot the slices using Matplotlib\n", "fig, axes = plt.subplots(1, 3, figsize=(15, 5))\n", "\n", "# Axial slice\n", "axes[0].imshow(np.rot90(data1[:, :, slice_idx1, 0]), cmap=\"Greys_r\", vmin=0, vmax=700)\n", "axes[0].imshow(np.rot90(rgb_overlay[:, :, slice_idx1]), alpha=0.7)\n", "axes[0].axis(\"off\")\n", "axes[0].set_title('Axial Slice')\n", "\n", "# Coronal slice\n", "axes[1].imshow(np.rot90(data1[:, slice_idx2, :, 0]), cmap=\"Greys_r\", vmin=0, vmax=700)\n", "axes[1].imshow(np.rot90(rgb_overlay[:, slice_idx2, :]), alpha=0.7)\n", "axes[1].axis(\"off\")\n", "axes[1].set_title('Coronal Slice')\n", "\n", "# Sagittal slice\n", "axes[2].imshow(np.rot90(data1[slice_idx3, :, :, 0]), cmap=\"Greys_r\", vmin=0, vmax=700)\n", "axes[2].imshow(np.rot90(rgb_overlay[slice_idx3, :, :]), alpha=0.7)\n", "axes[2].axis(\"off\")\n", "axes[2].set_title('Sagittal Slice')\n", "\n", "# Adjust layout and display\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fiber Orientation Density (FOD)\n", "We will now use the basis functions generated above to create Fiber Orientation Densities, or FODs. These are estimates of the amount of diffusion in each of three orthogonal directions. As described in the introductory chapter, these are analogous to the tensors that are used in traditional diffusion studies. However, MRtrix allows for the estimation of multiple crossing fibers within a single voxel, and can resolve the diffusion signal into multiple directions.\n", "\n", "To do this, we will use the command ```dwi2fod``` to apply the basis functions to the diffusion data. The “-mask” option specifies which voxels we will use; this is simply to restrict our analysis to brain voxels and reduce the computing time. The “.mif” files specified after each basis function will output an FOD image for that tissue type:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "dwi2fod: [100%] preloading data for \"sub-02_den_preproc_unbiased.mif\"\u001b[0K[0K\u001b[?7h\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\n", "dwi2fod: [100%] performing MSMT CSD (4 shells, 3 tissues)\u001b[0K[0K\u001b[?7h\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\n" ] } ], "source": [ "! dwi2fod msmt_csd sub-02_den_preproc_unbiased.mif -mask mask.mif wm.txt wmfod.mif gm.txt gmfod.mif csf.txt csffod.mif" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In order to view these FODs, we will combine them into a single image. The command ```mrconvert``` will extract the first image from the wmfod.mif file, which is the image with a b-value of 0. The output of this command is then used as the input into an ```mrcat``` command which combines the FOD images from all three tissue types into a single image that we will call “vf.mif”:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "mrconvert: [100%] copying from \"wmfod.mif\" to \"/tmp/mrtrix-tmp-yYcy0U.mif\"\u001b[0K[0K\u001b[?7h\u001b[?7l\n", "mrcat: [100%] concatenating \"csffod.mif\"\u001b[0K[0K\u001b[?7h\u001b[?7l\u001b[?7l\n", "mrcat: [100%] concatenating \"gmfod.mif\"\u001b[0K[0K\u001b[?7h\u001b[?7l\n", "mrcat: [100%] concatenating \"/tmp/mrtrix-tmp-yYcy0U.mif\"\u001b[0K[0K\u001b[?7h\u001b[?7l\u001b[?7l\n" ] } ], "source": [ "! mrconvert -coord 3 0 wmfod.mif - | mrcat csffod.mif gmfod.mif - vf.mif" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "####\n", "The white matter FODs can then be overlaid on this image, so that we can observe whether the white matter FODs do indeed fall within the white matter, and also whether they are along the orientations that we would expect. The command for mrview is: \n", "```javascript\n", "mrview vf.mif -odf.load_sh wmfod.mif\n", "```\n", "\n", "####" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, we will use the python libraries **DIPY** and **Matplotlib** to visualize white matter FODs of an axial slice. However, we will also need to reconstruct the fiber orientation distribution function (fODF) with DIPY in two steps:\n", "\n", "- Estimation of the fiber response function\n", "- Use the response function to reconstruct the fODF\n", "\n", "#### Reconstruction with Constrained Spherical Deconvolution with DIPY\n", "\n", "A straightforward method for estimating the fiber response function is to focus on regions of the brain known to contain single, coherent fiber populations. For instance, placing a region of interest (ROI) at the center of the brain can capture single fibers from the corpus callosum. The ```auto_response_ssst``` function computes the fractional anisotropy (FA) within a cuboidal ROI, defined by radii specified in roi_radii. It then derives the response function based on voxels with an FA value exceeding 0.7." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "data_path = convert_mif_to_nii('sub-02_den_preproc_unbiased.mif')\n", "\n", "# Load the data\n", "data, affine = load_nifti(data_path) \n", "\n", "bvals_ap = np.loadtxt('sub-02_AP.bval')\n", "bvecs_ap = np.loadtxt('sub-02_AP.bvec').T # Make sure to transpose bvecs for correct shape\n", "\n", "# Create the gradient table\n", "gtab_ap = gradient_table(bvals_ap, bvecs_ap)\n", "\n", "response, ratio = auto_response_ssst(gtab_ap, data, roi_radii=10, fa_thr=0.7)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To ensure the response function is accurate, we can visualize its orientation distribution function (ODF):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "scene = window.Scene()\n", "evals = response[0]\n", "evecs = np.array([[0, 1, 0], [0, 0, 1], [1, 0, 0]]).T\n", "\n", "response_odf = single_tensor_odf(default_sphere.vertices, evals, evecs)\n", "# transform our data from 1D to 4D\n", "response_odf = response_odf[None, None, None, :]\n", "response_actor = actor.odf_slicer(response_odf, sphere=default_sphere,\n", " colormap='plasma')\n", "scene.add(response_actor)\n", "\n", "csd_response = window.snapshot(\n", " scene, fname=os.path.join(os.getcwd(), 'csd_response.png'), size=(200, 200),\n", " offscreen=True)\n", "\n", "fig, axes = plt.subplots(figsize=(3,3))\n", "axes.imshow(csd_response, cmap=\"plasma\")\n", "axes.axis(\"off\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After estimating the response function, we can begin the deconvolution process by fitting the CSD model to the data." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Set up the scene for visualization\n", "scene = window.Scene()\n", "\n", "# Fit the CSD model to the data\n", "csd_model = ConstrainedSphericalDeconvModel(gtab_ap, response)\n", "\n", "# Fit on a smaller portion of the data (optional for faster computation)\n", "data_small = data[30:60, 40:70, 29:30] \n", "csd_fit = csd_model.fit(data_small) \n", "csd_odf = csd_fit.odf(default_sphere)\n", "\n", "fodf_spheres = actor.odf_slicer(csd_odf, sphere=default_sphere, \n", " scale=0.9, norm=False, \n", " colormap='plasma')\n", "scene.add(fodf_spheres)\n", "\n", "csd_odfs = window.snapshot(\n", " scene, fname=os.path.join(os.getcwd(), 'csd_odfs.png'), size=(600, 600),\n", " offscreen=True)\n", "\n", "fig, axes = plt.subplots(figsize=(7,7))\n", "axes.imshow(csd_odfs, cmap=\"plasma\")\n", "axes.axis(\"off\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using DIPY's ```peaks_from_model```, we can identify the peak directions (maxima) of the ODFs. We will then visualize the ODFs and their corresponding peaks in the same space." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "csd_peaks = peaks_from_model(model=csd_model,\n", " data=data_small,\n", " sphere=default_sphere,\n", " relative_peak_threshold=.5,\n", " min_separation_angle=25,\n", " parallel=True,\n", " num_processes=10)\n", "\n", "scene.clear()\n", "fodf_peaks = actor.peak_slicer(csd_peaks.peak_dirs, peaks_values=csd_peaks.peak_values)\n", "scene.add(fodf_peaks)\n", "\n", "fodf_spheres.GetProperty().SetOpacity(0.6)\n", "scene.add(fodf_spheres)\n", "\n", "csd_both = window.snapshot(\n", " scene, fname=os.path.join(os.getcwd(), 'csd_both.png'), size=(600, 600),\n", " offscreen=True)\n", "\n", "fig, axes = plt.subplots(figsize=(7,7))\n", "axes.imshow(csd_both, cmap=\"plasma\")\n", "axes.axis(\"off\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Normalization\n", "\n", "In order to make the comparisons valid across subjects, we will need to normalize the FODs. This ensures that any differences we see are not due to intensity differences in the image, similar to how we correct for the size of the brain when comparing volumetric differences across subjects.\n", "\n", "To normalize the data, we will use the ```mtnormalise``` command. This requires an input and output for each tissue type, as well as a mask to restrict the analysis to brain voxels:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "mtnormalise: [100%] performing log-domain intensity normalisation\u001b[0K[0K\u001b[?7h\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\n" ] } ], "source": [ "! mtnormalise wmfod.mif wmfod_norm.mif gmfod.mif gmfod_norm.mif csffod.mif csffod_norm.mif -mask mask.mif" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we’ve correctly estimated the FODs for each tissue type, we are ready to begin laying down the foundation for our tractography analysis. The next step will be to determine the boundary between the grey matter and the white matter, which we will use as a starting point for our streamlines." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Creating the Tissue Boundaries\n", "We are almost ready to begin our streamline analysis, in which we will place seeds at random locations along the boundary between the grey matter and the white matter. A streamline will grow from each seed and trace a path from that seed region until it terminates in another region. Some of the streamlines will terminate in places that don’t make sense - for example, a streamline may terminate at the border of the ventricles. We will cull these “error” streamlines, and be left with a majority of streamlines that appear to connect distant grey matter regions.\n", "\n", "To do this, we will first need to create a boundary between the grey matter and the white matter. The MRtrix command ```5ttgen``` will use FSL’s FAST, along with other commands, to segment the anatomical image into five tissue types:\n", "\n", " Grey Matter;\n", "\n", " Subcortical Grey Matter (such as the amygdala and basal ganglia);\n", "\n", " White Matter;\n", "\n", " Cerebrospinal Fluid; and\n", "\n", " Pathological Tissue.\n", "\n", "Once we have segmented the brain into those tissue classes, we can then use the boundary as a mask to restrict where we will place our seeds." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Converting the Anatomical Image\n", "\n", "The anatomical image first needs to be converted to MRtrix format with ```mrconvert```." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "mrconvert: [100%] uncompressing image \"sub-CON02_ses-preop_T1w.nii.gz\"\u001b[0K[0K\u001b[?7h\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\n", "mrconvert: [100%] copying from \"sub-CON02_ses-preop_T1w.nii.gz\" to \"T1.mif\"\u001b[0K[0K\u001b[?7h\n" ] } ], "source": [ "! mrconvert sub-CON02_ses-preop_T1w.nii.gz T1.mif" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "5ttgen: \u001b[03;32m\u001b[0m\n", "5ttgen: \u001b[03;32mNote that this script makes use of commands / algorithms that have relevant articles for citation; INCLUDING FROM EXTERNAL SOFTWARE PACKAGES. Please consult the help page (-help option) for more information.\u001b[0m\n", "5ttgen: \u001b[03;32m\u001b[0m\n", "5ttgen: \u001b[03;32mGenerated scratch directory: /home/jovyan/Git_repositories/MRtrix_preprocessed/5ttgen-tmp-EHE0O1/\u001b[0m\n", "\u001b[03;36mCommand:\u001b[0m mrconvert /home/jovyan/Git_repositories/MRtrix_preprocessed/T1.mif /home/jovyan/Git_repositories/MRtrix_preprocessed/5ttgen-tmp-EHE0O1/input.mif\n", "5ttgen: \u001b[03;32mChanging to scratch directory (/home/jovyan/Git_repositories/MRtrix_preprocessed/5ttgen-tmp-EHE0O1/)\u001b[0m\n", "\u001b[03;36mCommand:\u001b[0m mrconvert input.mif T1.nii -strides -1,+2,+3\n", "\u001b[03;36mCommand:\u001b[0m maskfilter /opt/fsl-6.0.5.1/data/standard/MNI152_T1_1mm_brain_mask_dil.nii.gz dilate mni_mask.nii -npass 4\n", "\u001b[03;36mCommand:\u001b[0m standard_space_roi T1.nii T1_preBET.nii.gz -maskMASK mni_mask.nii -roiFOV\n", "\u001b[03;36mCommand:\u001b[0m bet T1_preBET.nii.gz T1_BET.nii.gz -f 0.15 -R\n", "\u001b[03;36mCommand:\u001b[0m fast T1_BET.nii.gz\n", "\u001b[03;36mCommand:\u001b[0m run_first_all -m none -s L_Accu,R_Accu,L_Caud,R_Caud,L_Pall,R_Pall,L_Puta,R_Puta,L_Thal,R_Thal -i T1.nii -o first\n", "5ttgen: \u001b[03;36m[100%]\u001b[0m \u001b[03;32mGenerating partial volume images for SGM structures\u001b[0m\u001b[0K\u001b[0K\u001b[?7h\u001b[?7h\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\u001b[?7l\n", "\u001b[03;36mCommand:\u001b[0m mrmath [mesh2voxel_*.mif (10 items)] sum - | mrcalc - 1.0 -min all_sgms.mif\n", "\u001b[03;36mCommand:\u001b[0m mrthreshold T1_BET_pve_2.nii.gz - -abs 0.001 | maskfilter - connect - -connectivity | mrcalc 1 - 1 -gt -sub remove_unconnected_wm_mask.mif -datatype bit\n", "\u001b[03;36mCommand:\u001b[0m mrcalc T1_BET_pve_0.nii.gz remove_unconnected_wm_mask.mif -mult csf.mif\n", "\u001b[03;36mCommand:\u001b[0m mrcalc 1.0 csf.mif -sub all_sgms.mif -min sgm.mif\n", "\u001b[03;36mCommand:\u001b[0m mrcalc 1.0 csf.mif sgm.mif -add -sub T1_BET_pve_1.nii.gz T1_BET_pve_2.nii.gz -add -div multiplier.mif\n", "\u001b[03;36mCommand:\u001b[0m mrcalc multiplier.mif -finite multiplier.mif 0.0 -if multiplier_noNAN.mif\n", "\u001b[03;36mCommand:\u001b[0m mrcalc T1_BET_pve_1.nii.gz multiplier_noNAN.mif -mult remove_unconnected_wm_mask.mif -mult cgm.mif\n", "\u001b[03;36mCommand:\u001b[0m mrcalc T1_BET_pve_2.nii.gz multiplier_noNAN.mif -mult remove_unconnected_wm_mask.mif -mult wm.mif\n", "\u001b[03;36mCommand:\u001b[0m mrcalc 0 wm.mif -min path.mif\n", "\u001b[03;36mCommand:\u001b[0m mrcat cgm.mif sgm.mif wm.mif csf.mif path.mif - -axis 3 | mrconvert - combined_precrop.mif -strides +2,+3,+4,+1\n", "\u001b[03;36mCommand:\u001b[0m mrmath combined_precrop.mif sum - -axis 3 | mrthreshold - - -abs 0.5 | mrgrid combined_precrop.mif crop result.mif -mask -\n", "\u001b[03;36mCommand:\u001b[0m mrconvert result.mif /home/jovyan/Git_repositories/MRtrix_preprocessed/5tt_nocoreg.mif\n", "\u001b[03;36mCommand:\u001b[0m 5ttcheck result.mif\n", "5ttgen: \u001b[03;32mChanging back to original directory (/home/jovyan/Git_repositories/MRtrix_preprocessed)\u001b[0m\n", "5ttgen: \u001b[03;32mDeleting scratch directory (/home/jovyan/Git_repositories/MRtrix_preprocessed/5ttgen-tmp-EHE0O1/)\u001b[0m\n" ] } ], "source": [ "! 5ttgen fsl T1.mif 5tt_nocoreg.mif" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "####\n", "This command will take about 10-15 minutes. If the segmentation has finished successfully, you can visualize the images with\n", "```javascript \n", "mrview 5tt_nocoreg.mif\n", "``` \n", "\n", "\n", "The output from 5ttgen fsl anat.mif 5tt_nocoreg.mif will be a single dataset with 5 volumes, one per tissue type. Check this image with ```mrview```, using the right and left arrow keys to toggle between tissue types. The tissue types are: GM, WM, CSF, subcortical GM, and pathological tissue. If no pathological tissue is detected, then that volume is blank.\n", "\n", "#####\n", "Here, we will use **AnyNiivue** for interactive visualization of the different tissue types:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "9bae53c4cd5544429772b597066909b4", "version_major": 2, "version_minor": 0 }, "text/plain": [ "VBox(children=(RadioButtons(description='Select Tissue Type:', options=(('GM', 0), ('Subcortical GM', 1), ('WM…" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Extract frames from 5tt_nocoreg.mif into temporary files\n", "output_files = []\n", "mif_file_path = '5tt_nocoreg.mif'\n", "\n", "for idx in range(5):\n", " with tempfile.NamedTemporaryFile(suffix=\".mif\", delete=False) as temp_file:\n", " temp_file.close()\n", " output_path = temp_file.name\n", " output_files.append(output_path)\n", " subprocess.run([\"mrconvert\", mif_file_path, \"-coord\", \"3\", str(idx), output_path, \"-force\", \"-quiet\"])\n", "\n", "# Initialize the AnyNiivue viewer\n", "nv = AnyNiivue()\n", "\n", "# Define the volumes with opacity set to 0 initially for all but the default layer (opacity = 1.0)\n", "volumes = [\n", " {\"path\": output_files[0], \"colormap\": \"gray\", \"opacity\": 1.0}, #default layer\n", " {\"path\": output_files[1], \"colormap\": \"gray\", \"opacity\": 0.0},\n", " {\"path\": output_files[2], \"colormap\": \"gray\", \"opacity\": 0.0},\n", " {\"path\": output_files[3], \"colormap\": \"gray\", \"opacity\": 0.0},\n", " {\"path\": output_files[4], \"colormap\": \"gray\", \"opacity\": 0.0},\n", "]\n", "\n", "# Load all the volumes into the viewer\n", "nv.load_volumes(volumes)\n", "\n", "# Function to update the opacity of the layers based on selection\n", "def update_layer(change):\n", " selected_index = change.new\n", "\n", " # Set opacity to 0 for all volumes\n", " for volume in nv.volumes:\n", " volume.opacity = 0.0\n", "\n", " # Set opacity to 1 for the selected volume\n", " nv.volumes[selected_index].opacity = 1.0\n", "\n", "# Create RadioButtons widget for selecting the active layer\n", "layer_selector = RadioButtons(\n", " options=[\n", " (\"GM\", 0),\n", " (\"Subcortical GM\", 1),\n", " (\"WM\", 2),\n", " (\"CSF\", 3),\n", " (\"Pathological Tissue\", 4),\n", " ],\n", " description='Select Tissue Type:',\n", " style={'description_width': 'initial'}\n", ")\n", "\n", "# Observe changes in the RadioButtons widget\n", "layer_selector.observe(update_layer, names='value')\n", "\n", "# Display the RadioButtons and the AnyNiivue viewer\n", "display(VBox([layer_selector, nv]))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "
\n", " \n", " !\n", " Note:\n", "
\n", "

\n", "If the segmentation step fails, this may be due to insufficient contrast between the tissue types; for example, some anatomical images are either very dark across both the grey and white matter, or very light across both tissue types. We can help the segmentation process by increasing the intensity contrast (also known as intensity normalization) between the tissues with a command like AFNI’s 3dUnifize, e.g.

\n", " 3dUnifize -input anat.nii -prefix anat_unifize.nii

\n", " The difference between the image before and after may be subtle, but it can prevent a segmentation error from being thrown.\n", "

\n", "
" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.6" } }, "nbformat": 4, "nbformat_minor": 4 }