dbfold_161_2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# protein 161 dbfold test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import dbfold\n",
    "from dbfold.dbfold import Protein as Protein\n",
    "import dbfold.analyze_structures as analyze_structures\n",
    "import dbfold.load_data as load_data\n",
    "import dbfold.compute_PMF as compute_PMF\n",
    "import dbfold.folding_rates as folding_rates\n",
    "import dbfold.kinetic_model as kinetic_model\n",
    "import dbfold.nonnative_states as nonnative_states\n",
    "import matplotlib.pyplot as plt\n",
    "import importlib"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First thing to do is visualize native substructures and also count total number of contacts. This is useful so that we can properly set up umbrella biasing in the energy minimizaiton simulations. Make sure to set the d_cutoff and min_seq_separation values to whatever you plan to use in the simulation!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "hi\n",
      "This protein has 102 contacts\n"
     ]
    }
   ],
   "source": [
    "unequilibrated_protein161 = Protein('protein 161','cts1_2uy3_len_161/cts1_2uy3_len_161.pdb')\n",
    "print('This protein has {} contacts'.format(unequilibrated_protein161.count_contacts(d_cutoff=6, min_seq_separation=8)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Running simulations and extracting data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step is to run minimization simulations, to equilibrate the protein in the potential. This can be done with the MCPU Monte-Carlos package, which is included as part of dbfold--detailed instructions are available in the MCPU Readme. Upon obtaining a minimized PDB structure, we count its contacts (for use in setting up umbrella biasing) and identify substructures. It may be worth experimenting to find parameters that produce sensible substructures for a given protein. We will both visualize the substructures on a contact map, and also produce a .pml file, which contains a series of commands to visualize these substructures directly on the equlibrated structure in PYMOL. This .pml file should be be opened along with the equilibrated .pdb in PYMOL.\n",
    "We also save the substructure as a file 1igd_data/Substructures.dat\n",
    "\n",
    "We note that the d_cutoff value used to generate substructures need not be the same as the one used to define contacts in the umbrella simulations, and likewise for min_seq_separation. In general, there may be two definitions of contacts: one that is used duruing simulation to perform umbrella biasing, and one that is used to produce substructures during analysis.\n",
    "\n",
    "Finally we set the value of f to 1.7 as in the paper\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "hi\n",
      "This protein has 228 contacts\n"
     ]
    }
   ],
   "source": [
    "protein161 = Protein('protein 161','../Olivia_data/cts1_Emin.pdb')\n",
    "print('This protein has {} contacts'.format(protein161.count_contacts(d_cutoff=6, min_seq_separation=8)))\n",
    "protein161.f = 1.7\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "hi\n",
      "This protein has 228 contacts\n"
     ]
    }
   ],
   "source": [
    "%matplotlib qt\n",
    "full_protein = Protein('protein F','cts1_Emin.pdb')\n",
    "print('This protein has {} contacts'.format(full_protein.count_contacts(d_cutoff=6, min_seq_separation=8)))\n",
    "full_protein.generate_subs(d_cutoff=6, min_seq_separation=8, contact_sep_thresh=6, min_clustersize=8, savepath = 'cts1_full_substructures.dat')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we run equilibrium and unfolding simulations. As above, details for how to run these simulations in MCPU are available in the README.\n",
    "We then extract data from these simulations for further analysis using the following steps:\n",
    "\n",
    "   1. We run the script read_log_files.py to extract energies, number of native contacts, RMSD, and other variables of interest from the log file for each MC step at which these values were printed. It is only necessary that this be done for equilibrium simulations. The output is saved in a .dat file (for example, 1igd_data/Equilibrium_log_data.dat).\n",
    "\n",
    "   2. We run get_snapshot_configs.py to obtain substructure scores from PDB files, which are necessary to assign each PDB snapshot to a topological configuration. This script can either generate the substructures de novo, which requires that we input the parameters to create the substructures used above, or it can be given a file Substructures.dat, which can be produced when generating substructures above if savepath is set to something other than None (the default value). The output is saved in a .dat file (for example, 1igd_data/Equilibrium_scores.dat and 1igd_data/Unfolding_scores.dat). This step should be run for both equilibrium and unfolding simulations.\n",
    "\n",
    "\n",
    "IMPORTANT: The dbfold analysis package requires that the MCPU parameters MC_PRINT_STEPS and \n",
    "MC_PDB_PRINT_STEPS, which are set in the cfg files, have the same value. This ensures that the MC steps printed in the log files correspond to those for which PDB files are produced.\n",
    "\n",
    "This algorithm can also be used with alternative simulation packages, including MD software, but the above scripts must be adapted to work with the PDB naming conventions and log file format for whichever program is used.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Equilibrium properties"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's load the test data for protein 161. First, indicate the directory where the equilibrium data is stored, which in our case is 1igd_data. By default, this assumes the log data filename is Equilibrium_log_data.dat and the scores filename is Equilibrium_scores.dat, but these can be changed via the optional arguments log_data_filename and score_filename, respectively"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The equilibrium log datapath has been set as ../Olivia_data/log_cont.dat and the scorepath as ../Olivia_data/scores_cont.dat\n"
     ]
    }
   ],
   "source": [
    "protein161.set_eq_path('../Olivia_data', log_data_filename = \"log_cont.dat\", score_filename = \"scores_cont.dat\"  )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We note that this file contains the data for the last 300 milion MC steps of a simulation that was run for a total of 1.2 billon MC steps.\n",
    "Now, we assess convergence for the equilibrium simulations. To do this, we plot the simulation energy as a function of MC step, averaged over a sliding window of 50 million MC steps. We will do this at a number of temperatures, T = 0.8, 0.9, and 0.975"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#fig, ax = plt.subplots(1,3)\n",
    "#plt.subplots_adjust(wspace = 3)\n",
    "compute_PMF.plotTrajectoryTimecourses(protein161.log_datapath, 'natives', \n",
    "                          ['0.800_40.', '0.800_30.', '0.800_20.', '0.800_10.', '0.800_0.'], \n",
    "                          colors = ['b','g','r','c','m'], \n",
    "                          labels = ['S = 40', 'S = 30','S = 20','S = 10','S = 0'], sliding_windowlen = 50, \n",
    "                                     legend_fontsize = 15)\n",
    "\n",
    "compute_PMF.plotTrajectoryTimecourses(protein161.log_datapath, 'natives', \n",
    "                          ['0.900_40.', '0.900_30.', '0.900_20.', '0.900_10.', '0.900_0.'], \n",
    "                          colors = ['b','g','r','c','m'], \n",
    "                          labels = ['S = 40', 'S = 30','S = 20','S = 10','S = 0'], sliding_windowlen = 50, \n",
    "                                     legend_fontsize = 15)\n",
    "\n",
    "compute_PMF.plotTrajectoryTimecourses(protein161.log_datapath, 'natives', \n",
    "                          ['0.975_40.', '0.975_30.', '0.975_20.', '0.975_10.', '0.975_0.'], \n",
    "                          colors = ['b','g','r','c','m'], \n",
    "                          labels = ['S = 40', 'S = 30','S = 20','S = 10','S = 0'], sliding_windowlen = 50, \n",
    "                                     legend_fontsize = 15)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'cts1_2uy3_len_161/log_cont.dat'"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "protein161.log_datapath"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the above plots, it appears that beyond 150 million MC steps, these energies generally cease changing substantially. Thus we will compute all equilibium thermodynamic properties starting with this timestep.\n",
    "We do this below. As in the paper, we define a substructure using a value of f=1.7, where f is the maximum value of the average distance between residues assigned to a substructure within a snapshot\n",
    "dived by that same average distance in the native file such that a substructure is considered formed\n",
    "\n",
    "The k_bias value should always correspond to the one used in umbrella simulations. This assumes biasing was performed along native contacts\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Need to compute PMFs using MBAR!\n",
      "Computing PMFs vs natives..\n",
      "# loading data...\n",
      "# calculating potential...\n",
      "# Computing normalizations...\n",
      "Computing state free energies...\n",
      "Computing free energy at T=0.4\n",
      "Computing free energy at T=0.425\n",
      "Computing free energy at T=0.45\n",
      "Computing free energy at T=0.475\n",
      "Computing free energy at T=0.5\n",
      "Computing free energy at T=0.525\n",
      "Computing free energy at T=0.55\n",
      "Computing free energy at T=0.575\n",
      "Computing free energy at T=0.6\n",
      "Computing free energy at T=0.625\n",
      "Computing free energy at T=0.65\n",
      "Computing free energy at T=0.675\n",
      "Computing free energy at T=0.7\n",
      "Computing free energy at T=0.725\n",
      "Computing free energy at T=0.75\n",
      "Computing free energy at T=0.775\n",
      "Computing free energy at T=0.8\n",
      "Computing free energy at T=0.825\n",
      "Computing free energy at T=0.85\n",
      "Computing free energy at T=0.875\n",
      "Computing free energy at T=0.9\n",
      "Computing free energy at T=0.925\n",
      "Computing free energy at T=0.95\n",
      "Computing free energy at T=0.975\n",
      "Computing free energy at T=1.0\n",
      "Computing PMFs vs topological configuration..\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/amirbitran/Dropbox/CurrentStuff/Harvard/Shakhnovich Lab/dbfold/dbfold/load_data.py:114: RuntimeWarning: Mean of empty slice\n",
      "  mean_substructure_distances.append(np.nanmean(x[np.where(x)]))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# loading data...\n",
      "# calculating potential...\n",
      "# Computing normalizations...\n",
      "Computing state free energies...\n",
      "Computing free energy at T=0.4\n",
      "Computing free energy at T=0.425\n",
      "Computing free energy at T=0.45\n",
      "Computing free energy at T=0.475\n",
      "Computing free energy at T=0.5\n",
      "Computing free energy at T=0.525\n",
      "Computing free energy at T=0.55\n",
      "Computing free energy at T=0.575\n",
      "Computing free energy at T=0.6\n",
      "Computing free energy at T=0.625\n",
      "Computing free energy at T=0.65\n",
      "Computing free energy at T=0.675\n",
      "Computing free energy at T=0.7\n",
      "Computing free energy at T=0.725\n",
      "Computing free energy at T=0.75\n",
      "Computing free energy at T=0.775\n",
      "Computing free energy at T=0.8\n",
      "Computing free energy at T=0.825\n",
      "Computing free energy at T=0.85\n",
      "Computing free energy at T=0.875\n",
      "Computing free energy at T=0.9\n",
      "Computing free energy at T=0.925\n",
      "Computing free energy at T=0.95\n",
      "Computing free energy at T=0.975\n",
      "Computing free energy at T=1.0\n",
      "Finished computing PMFs!\n"
     ]
    }
   ],
   "source": [
    "protein161.obtain_PMFs(eq_step=200000000, k_bias = 0.02)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now plot a melting curve that shows fraction of native contacts (using the contacts definition in the umbrella simulations) as a function of temperature. There are two ways of doing this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib qt\n",
    "protein161.plot_melting_curve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "compute_PMF.plot_natives_vs_temp(['cts1_2uy3_len_161'], ['protein 161'], norm = 'Fraction', d_cutoff =6, \n",
    "                     Native_structures = ['cts1_2uy3_len_161/cts1_2uy3_len_161_0.100_20_Emin.pdb'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "importlib.reload(dbfold.compute_PMF)\n",
    "\n",
    "compute_PMF.plot_substructure_PMF('cts1_2uy3_len_161/PMFs.dat', [0.4])\n",
    "compute_PMF.plot_substructure_PMF('cts1_2uy3_len_161/PMFs.dat', [0.6])\n",
    "compute_PMF.plot_substructure_PMF('cts1_2uy3_len_161/PMFs.dat', [0.7])\n",
    "compute_PMF.plot_substructure_PMF('cts1_2uy3_len_161/PMFs.dat', [0.9])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Nonnative contacts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To produce nonnative contact maps, we must run the script get_snapshot_contact_maps.py, which extracts contact maps (including both native and nonnative contacts--the former can be filtered out in subsequent analyses) for snapshots of interest (for instnace, those assigned to a specific subset of topological configuration at a specific temperature range). These contact maps are saved in a .dat file (e.g. 1igd_data/Nonnative_maps(b)c.dat), and they are analyzed below.\n",
    "\n",
    "For the case of protein 161 consider the nonnative contacts among all snapshots assigned to topological configurations c or bc at physiologically reasonable temperatures. We can visualize the average nonnative contact map among such snapshots by accessing the data stored in 1igd_data/Nonnative_maps_(b)c.dat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "nonnative_states.visualize_nonnatives('cts1_2uy3_len_161/cluster.dat', protein161.native_structure, d_cutoff = 6, min_seq_separation = 8 ,filter_distance = 2,cmap = 'Reds')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also cluster these snapshots into structurally distinct, recurrent nonnative states.\n",
    "You can experiment with the thresh, which is the maximum allowed Hamming distance between the nonnative contact maps for two snapshots such that they are clustered together"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean number of nonnative contacts is 33.05731225296443\n"
     ]
    }
   ],
   "source": [
    "clusters, percentages, mean_maps = nonnative_states.cluster_nonnatives('cts1_2uy3_len_161/cluster.dat', protein161.native_structure, d_cutoff = 6, \n",
    "                                                                       filter_distance = 2,min_seq_separation = 3 ,cmap = 'Reds', thresh = 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['cts1_2uy3_len_161_0.675_20.150500000',\n",
       " 'cts1_2uy3_len_161_0.650_0.85000000',\n",
       " 'cts1_2uy3_len_161_0.650_10.78500000',\n",
       " 'cts1_2uy3_len_161_0.675_0.93000000',\n",
       " 'cts1_2uy3_len_161_0.675_10.113500000',\n",
       " 'cts1_2uy3_len_161_0.650_10.87500000',\n",
       " 'cts1_2uy3_len_161_0.625_0.27000000',\n",
       " 'cts1_2uy3_len_161_0.650_10.68500000',\n",
       " 'cts1_2uy3_len_161_0.625_0.47500000',\n",
       " 'cts1_2uy3_len_161_0.600_0.30500000',\n",
       " 'cts1_2uy3_len_161_0.675_10.114000000',\n",
       " 'cts1_2uy3_len_161_0.675_20.105000000',\n",
       " 'cts1_2uy3_len_161_0.625_0.44500000',\n",
       " 'cts1_2uy3_len_161_0.625_0.29500000',\n",
       " 'cts1_2uy3_len_161_0.650_0.85500000',\n",
       " 'cts1_2uy3_len_161_0.650_0.84500000',\n",
       " 'cts1_2uy3_len_161_0.675_20.104500000',\n",
       " 'cts1_2uy3_len_161_0.675_10.137500000',\n",
       " 'cts1_2uy3_len_161_0.650_0.79500000',\n",
       " 'cts1_2uy3_len_161_0.600_0.61000000',\n",
       " 'cts1_2uy3_len_161_0.650_0.86500000',\n",
       " 'cts1_2uy3_len_161_0.675_10.140000000',\n",
       " 'cts1_2uy3_len_161_0.650_10.87000000',\n",
       " 'cts1_2uy3_len_161_0.625_0.46500000',\n",
       " 'cts1_2uy3_len_161_0.625_0.71000000',\n",
       " 'cts1_2uy3_len_161_0.650_0.91000000',\n",
       " 'cts1_2uy3_len_161_0.675_0.92500000',\n",
       " 'cts1_2uy3_len_161_0.625_0.72000000',\n",
       " 'cts1_2uy3_len_161_0.650_10.89000000',\n",
       " 'cts1_2uy3_len_161_0.625_0.28000000',\n",
       " 'cts1_2uy3_len_161_0.675_10.138500000',\n",
       " 'cts1_2uy3_len_161_0.675_10.114500000',\n",
       " 'cts1_2uy3_len_161_0.675_10.112000000',\n",
       " 'cts1_2uy3_len_161_0.650_0.83000000',\n",
       " 'cts1_2uy3_len_161_0.650_0.94500000',\n",
       " 'cts1_2uy3_len_161_0.675_20.103000000',\n",
       " 'cts1_2uy3_len_161_0.650_10.78000000',\n",
       " 'cts1_2uy3_len_161_0.675_10.141000000',\n",
       " 'cts1_2uy3_len_161_0.675_10.129500000',\n",
       " 'cts1_2uy3_len_161_0.625_0.28500000',\n",
       " 'cts1_2uy3_len_161_0.650_10.68000000',\n",
       " 'cts1_2uy3_len_161_0.600_0.31500000',\n",
       " 'cts1_2uy3_len_161_0.625_0.27500000',\n",
       " 'cts1_2uy3_len_161_0.600_0.60500000',\n",
       " 'cts1_2uy3_len_161_0.625_0.71500000',\n",
       " 'cts1_2uy3_len_161_0.650_0.84000000',\n",
       " 'cts1_2uy3_len_161_0.650_0.86000000',\n",
       " 'cts1_2uy3_len_161_0.625_0.80500000',\n",
       " 'cts1_2uy3_len_161_0.650_10.69500000',\n",
       " 'cts1_2uy3_len_161_0.675_0.75500000',\n",
       " 'cts1_2uy3_len_161_0.625_0.65000000',\n",
       " 'cts1_2uy3_len_161_0.675_0.15000000',\n",
       " 'cts1_2uy3_len_161_0.625_0.43500000',\n",
       " 'cts1_2uy3_len_161_0.650_10.90500000',\n",
       " 'cts1_2uy3_len_161_0.650_0.91500000',\n",
       " 'cts1_2uy3_len_161_0.625_0.45500000',\n",
       " 'cts1_2uy3_len_161_0.675_10.136000000',\n",
       " 'cts1_2uy3_len_161_0.650_0.93500000',\n",
       " 'cts1_2uy3_len_161_0.675_10.140500000',\n",
       " 'cts1_2uy3_len_161_0.675_10.139000000',\n",
       " 'cts1_2uy3_len_161_0.675_0.72500000',\n",
       " 'cts1_2uy3_len_161_0.650_0.83500000',\n",
       " 'cts1_2uy3_len_161_0.675_20.104000000',\n",
       " 'cts1_2uy3_len_161_0.600_0.49000000',\n",
       " 'cts1_2uy3_len_161_0.625_0.81000000',\n",
       " 'cts1_2uy3_len_161_0.675_10.130000000',\n",
       " 'cts1_2uy3_len_161_0.600_0.59500000',\n",
       " 'cts1_2uy3_len_161_0.650_10.70000000',\n",
       " 'cts1_2uy3_len_161_0.625_0.64000000',\n",
       " 'cts1_2uy3_len_161_0.650_0.41500000',\n",
       " 'cts1_2uy3_len_161_0.675_20.102500000',\n",
       " 'cts1_2uy3_len_161_0.625_0.66500000',\n",
       " 'cts1_2uy3_len_161_0.650_0.42500000',\n",
       " 'cts1_2uy3_len_161_0.675_20.103500000',\n",
       " 'cts1_2uy3_len_161_0.650_0.94000000',\n",
       " 'cts1_2uy3_len_161_0.600_0.48000000',\n",
       " 'cts1_2uy3_len_161_0.625_0.47000000',\n",
       " 'cts1_2uy3_len_161_0.650_10.88500000',\n",
       " 'cts1_2uy3_len_161_0.625_0.82000000',\n",
       " 'cts1_2uy3_len_161_0.625_0.45000000',\n",
       " 'cts1_2uy3_len_161_0.650_0.82500000',\n",
       " 'cts1_2uy3_len_161_0.650_10.77000000',\n",
       " 'cts1_2uy3_len_161_0.650_10.110500000',\n",
       " 'cts1_2uy3_len_161_0.650_10.79000000',\n",
       " 'cts1_2uy3_len_161_0.625_0.62000000',\n",
       " 'cts1_2uy3_len_161_0.650_10.88000000',\n",
       " 'cts1_2uy3_len_161_0.625_0.64500000',\n",
       " 'cts1_2uy3_len_161_0.650_0.67000000',\n",
       " 'cts1_2uy3_len_161_0.650_10.90000000',\n",
       " 'cts1_2uy3_len_161_0.675_10.113000000',\n",
       " 'cts1_2uy3_len_161_0.650_0.43000000',\n",
       " 'cts1_2uy3_len_161_0.650_10.76500000',\n",
       " 'cts1_2uy3_len_161_0.625_0.81500000',\n",
       " 'cts1_2uy3_len_161_0.675_10.137000000',\n",
       " 'cts1_2uy3_len_161_0.650_0.42000000',\n",
       " 'cts1_2uy3_len_161_0.650_0.92000000',\n",
       " 'cts1_2uy3_len_161_0.600_0.48500000',\n",
       " 'cts1_2uy3_len_161_0.675_0.74500000',\n",
       " 'cts1_2uy3_len_161_0.625_0.61500000',\n",
       " 'cts1_2uy3_len_161_0.650_0.80000000',\n",
       " 'cts1_2uy3_len_161_0.625_0.46000000',\n",
       " 'cts1_2uy3_len_161_0.675_10.141500000',\n",
       " 'cts1_2uy3_len_161_0.675_10.115000000',\n",
       " 'cts1_2uy3_len_161_0.650_0.67500000',\n",
       " 'cts1_2uy3_len_161_0.650_10.89500000',\n",
       " 'cts1_2uy3_len_161_0.625_0.30000000',\n",
       " 'cts1_2uy3_len_161_0.650_0.76000000']"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "clusters[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['cts1_2uy3_len_161_0.675_10.235000000',\n",
       " 'cts1_2uy3_len_161_0.600_0.215500000',\n",
       " 'cts1_2uy3_len_161_0.625_0.208000000',\n",
       " 'cts1_2uy3_len_161_0.625_0.185000000',\n",
       " 'cts1_2uy3_len_161_0.650_0.190500000',\n",
       " 'cts1_2uy3_len_161_0.650_0.218500000',\n",
       " 'cts1_2uy3_len_161_0.625_0.172500000',\n",
       " 'cts1_2uy3_len_161_0.675_0.178500000',\n",
       " 'cts1_2uy3_len_161_0.650_0.173000000',\n",
       " 'cts1_2uy3_len_161_0.625_0.217500000',\n",
       " 'cts1_2uy3_len_161_0.625_0.184500000',\n",
       " 'cts1_2uy3_len_161_0.650_0.195000000',\n",
       " 'cts1_2uy3_len_161_0.650_0.211500000',\n",
       " 'cts1_2uy3_len_161_0.650_0.210000000',\n",
       " 'cts1_2uy3_len_161_0.625_0.212500000',\n",
       " 'cts1_2uy3_len_161_0.650_0.191500000',\n",
       " 'cts1_2uy3_len_161_0.625_0.171000000',\n",
       " 'cts1_2uy3_len_161_0.675_0.181000000',\n",
       " 'cts1_2uy3_len_161_0.650_0.205500000',\n",
       " 'cts1_2uy3_len_161_0.650_0.182500000',\n",
       " 'cts1_2uy3_len_161_0.675_0.220500000',\n",
       " 'cts1_2uy3_len_161_0.625_0.213500000',\n",
       " 'cts1_2uy3_len_161_0.650_0.211000000',\n",
       " 'cts1_2uy3_len_161_0.675_0.176500000',\n",
       " 'cts1_2uy3_len_161_0.650_0.209500000',\n",
       " 'cts1_2uy3_len_161_0.675_0.179500000',\n",
       " 'cts1_2uy3_len_161_0.650_0.183500000',\n",
       " 'cts1_2uy3_len_161_0.675_0.219500000',\n",
       " 'cts1_2uy3_len_161_0.625_0.172000000',\n",
       " 'cts1_2uy3_len_161_0.650_0.191000000',\n",
       " 'cts1_2uy3_len_161_0.650_0.207000000',\n",
       " 'cts1_2uy3_len_161_0.675_0.223500000',\n",
       " 'cts1_2uy3_len_161_0.600_0.214000000',\n",
       " 'cts1_2uy3_len_161_0.650_0.209000000',\n",
       " 'cts1_2uy3_len_161_0.650_0.174000000',\n",
       " 'cts1_2uy3_len_161_0.650_0.206500000',\n",
       " 'cts1_2uy3_len_161_0.650_0.192000000',\n",
       " 'cts1_2uy3_len_161_0.650_0.174500000',\n",
       " 'cts1_2uy3_len_161_0.650_0.193000000',\n",
       " 'cts1_2uy3_len_161_0.625_0.207500000',\n",
       " 'cts1_2uy3_len_161_0.625_0.186500000',\n",
       " 'cts1_2uy3_len_161_0.675_0.176000000',\n",
       " 'cts1_2uy3_len_161_0.650_0.167000000',\n",
       " 'cts1_2uy3_len_161_0.675_0.175500000',\n",
       " 'cts1_2uy3_len_161_0.675_0.221000000',\n",
       " 'cts1_2uy3_len_161_0.675_0.181500000',\n",
       " 'cts1_2uy3_len_161_0.675_10.226500000',\n",
       " 'cts1_2uy3_len_161_0.675_10.231500000',\n",
       " 'cts1_2uy3_len_161_0.625_0.171500000',\n",
       " 'cts1_2uy3_len_161_0.650_0.167500000',\n",
       " 'cts1_2uy3_len_161_0.650_0.184000000',\n",
       " 'cts1_2uy3_len_161_0.625_0.186000000',\n",
       " 'cts1_2uy3_len_161_0.650_0.206000000',\n",
       " 'cts1_2uy3_len_161_0.650_0.194500000',\n",
       " 'cts1_2uy3_len_161_0.675_0.175000000',\n",
       " 'cts1_2uy3_len_161_0.625_0.212000000',\n",
       " 'cts1_2uy3_len_161_0.650_0.195500000',\n",
       " 'cts1_2uy3_len_161_0.675_0.179000000',\n",
       " 'cts1_2uy3_len_161_0.675_0.178000000',\n",
       " 'cts1_2uy3_len_161_0.650_0.170500000',\n",
       " 'cts1_2uy3_len_161_0.650_0.219000000',\n",
       " 'cts1_2uy3_len_161_0.600_0.214500000',\n",
       " 'cts1_2uy3_len_161_0.625_0.213000000']"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "clusters[3]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['cts1_2uy3_len_161_0.600_0.232000000',\n",
       " 'cts1_2uy3_len_161_0.600_0.248500000',\n",
       " 'cts1_2uy3_len_161_0.600_0.265000000',\n",
       " 'cts1_2uy3_len_161_0.600_0.229000000',\n",
       " 'cts1_2uy3_len_161_0.600_0.229500000',\n",
       " 'cts1_2uy3_len_161_0.600_0.257500000',\n",
       " 'cts1_2uy3_len_161_0.600_0.254500000',\n",
       " 'cts1_2uy3_len_161_0.600_0.226500000',\n",
       " 'cts1_2uy3_len_161_0.600_0.227000000',\n",
       " 'cts1_2uy3_len_161_0.600_0.234000000',\n",
       " 'cts1_2uy3_len_161_0.600_0.263500000',\n",
       " 'cts1_2uy3_len_161_0.600_0.271500000',\n",
       " 'cts1_2uy3_len_161_0.600_0.235000000',\n",
       " 'cts1_2uy3_len_161_0.600_0.234500000',\n",
       " 'cts1_2uy3_len_161_0.600_0.256500000',\n",
       " 'cts1_2uy3_len_161_0.600_0.272000000',\n",
       " 'cts1_2uy3_len_161_0.600_0.247500000',\n",
       " 'cts1_2uy3_len_161_0.600_0.228500000',\n",
       " 'cts1_2uy3_len_161_0.600_0.248000000',\n",
       " 'cts1_2uy3_len_161_0.600_0.255500000',\n",
       " 'cts1_2uy3_len_161_0.600_0.274500000',\n",
       " 'cts1_2uy3_len_161_0.600_0.246500000',\n",
       " 'cts1_2uy3_len_161_0.600_0.244500000',\n",
       " 'cts1_2uy3_len_161_0.600_0.264500000',\n",
       " 'cts1_2uy3_len_161_0.600_0.256000000',\n",
       " 'cts1_2uy3_len_161_0.600_0.233000000',\n",
       " 'cts1_2uy3_len_161_0.600_0.227500000',\n",
       " 'cts1_2uy3_len_161_0.600_0.243500000',\n",
       " 'cts1_2uy3_len_161_0.600_0.249500000',\n",
       " 'cts1_2uy3_len_161_0.600_0.245000000',\n",
       " 'cts1_2uy3_len_161_0.600_0.249000000',\n",
       " 'cts1_2uy3_len_161_0.600_0.275000000',\n",
       " 'cts1_2uy3_len_161_0.600_0.244000000',\n",
       " 'cts1_2uy3_len_161_0.600_0.264000000',\n",
       " 'cts1_2uy3_len_161_0.600_0.228000000',\n",
       " 'cts1_2uy3_len_161_0.600_0.226000000',\n",
       " 'cts1_2uy3_len_161_0.600_0.247000000',\n",
       " 'cts1_2uy3_len_161_0.600_0.233500000',\n",
       " 'cts1_2uy3_len_161_0.600_0.231000000',\n",
       " 'cts1_2uy3_len_161_0.600_0.246000000',\n",
       " 'cts1_2uy3_len_161_0.600_0.230500000',\n",
       " 'cts1_2uy3_len_161_0.600_0.231500000',\n",
       " 'cts1_2uy3_len_161_0.600_0.232500000',\n",
       " 'cts1_2uy3_len_161_0.600_0.254000000',\n",
       " 'cts1_2uy3_len_161_0.600_0.245500000']"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "clusters[6]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clusters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Folding rate calculation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us now proceed with calculating folding rates. Note that this expects files from unfolding simulations to be titled ProteinName_temperature_TrajectoryNumber.MCStep, for instance 1igd_0.500_0.50000000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The unfolding path has been set as cts1_2uy3_len_161/Unfolding_scores.dat\n",
      "No unfolding data path has been set! Need to set unfolding path with self.set_unfolidng_path, or provide the directory using the optional argument unfolding_dir!\n",
      "The following parameters will be used to run HMM: \n",
      " f = 1.7 \n",
      " s = 0.01 \n",
      " m = 0.1 \n",
      " starting_state = abcdef \n",
      "Loading scores...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Olivia\\dbfold-master\\dbfold\\load_data.py:74: RuntimeWarning: Mean of empty slice\n",
      "  mean_substructure_distances.append(np.nanmean(x[np.where(x)]))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training HMM...\n",
      "Computing kinetic distances...\n",
      "From this resulting distance map, you may choose to cluster by choosing an adjancy threshold T_A that results in a significant separation of exchange timescales\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATkAAAF1CAYAAACTe53mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOydd3iUVdbAf2eSQAgttNBCbwECCQIiCAIivSg2FKwrLpa1L6yfFVfdRRAURcGOCiqCgIhU6R1BQpVOgNB7SUjCJOf7475JJsmkQMYA2ft7nveZzH3vfe87JWfOPfcUUVUsFouloOK60jdgsVgsfyVWyFkslgKNFXIWi6VAY4WcxWIp0FghZ7FYCjRWyFkslgKNFXIWi6VAY4WcxWIp0FghZ7FYCjRWyFkslgKN/5W+AYvF8r+FiCRfzjhVvSylzAo5i8WS3yQDkl+T2eWqxWLJV1TVX1X9PA8gEJiNEX7vZDzv9LksrJCzWCxXFBHxAyYCXYBE4EURecZX17dCzmKxXDFERIBxQC9gFRAGbANGiMg9PpnD5pOzWCxXChH5EngIWAfcrKpnRKQqsBwoC/RQ1d/yNIcVchaL5UogIh8DjwGbgXaqesLjXCNgERDgnFt7ufPY5arFYsl3RGQ4RsBtB27xFHAAqroRuA3jATJDRGpf7lxWyFkslnxFRN4EngP2AB1U9Yi3fqq6GOgHlMbsvF7efHa5arFY8hMRSQIOAG1UdW8u+j8OfGSdgS0Wy7XCYcwmQ44CDkBVR4tI5cudzGpyFoslXxGR+qr6Z77NZ4WcxWLJT0SkWm77qupeEfEHKqc8v+T5rJCzWCz5iWOTy1Xsqqq6RKQu8CdGXl2yXc7a5CwWS34zlksL0D9zGWNSsZqcxWIp0Fg/OYvFUqCxy1WLxZKvXOrGQ57ns8tVi8WSn1zqxkNe57OanMViyW/G4l3IBQORQDVMcH60LyazmpzFYrlqcPLLvQw8BVxvl6sWi6VAIiLLgX2qmufEmXZ31WKxXI2sADr44kJWyFkslquR0kAxX1zICjmLxXJVISKdgD7AJl9cz+6uWiyWfEVEFmRxyh+oAlQFkoDBPpnPbjxYLJb8RETceHchSQZOAauBYaq6yCfzWSFnsVgKMtYmZ7FYCjRWyFkslgKN3XiwWCxXBBGpB/wDuAEoCZwAVgIfq+oOn81jbXIWiyW/EZF+wJd4V7TcwGOq+pUv5rLLVYvFkq+ISEPgC0zd1R5ACefUu87zaOBTEWnhi/mskLNYLPnNCxjZ01NVZ6pqrNMep6ozgbZAPDDQF5NZIWexWPKb9sCsrOxuqnoYmAbc6IvJrJCzWCz5TQVgew599mPiV/OMFXIWiyW/iQUCcuhTBzjii8mskLNYLPnNQUx8qldEpB3QC5jhi8msC4nFYslXRGQ0cBdQUVUvOm1JmBxyRzEC7gjQTFUP5XU+q8lZLJb8ZgJGmHm6iCjQEuNC8jPQwhcCDqwmZ7FYrgJEpA1wHtju4VLim2tbIWexWK42ROQmoL2qvpHna1khZ7FYrgQiUhQoBfh5Of0M8CxQE7OUvexC01bIWSyWfEVEagLjMDa57IpMq+f5yy00bbOQWCyW/OY94HpgMbAPk+o8I5HOMTavk1lNzmKx5CsicgyYrar3ZdPnNWDw5WpvnlgXEovFkt+UIuewLp9hhZzFYslvBO9L1L8Ea5OzWCz5zcPAuhz6TMXklcsz1iZnsVgKNFaTs1gs+YqIVMttX1XdKyL+QOWU55c8n9XkLBZLfuIE42fnH5eKqrpEpC7wJ0ZeXfI+gtXkLBZLfjOWXAo5hzOXMSYVq8lZLJYCjXUhsVgsBRor5CwWS4HG2uQsFssVQUTqA48DNwDBwGlgFfCxqv7ps3msTc5iseQ3ItIf+BjvaZaSgCdV9TNfzGWXqxaLJV8RkRbAaOAs8DxQHyjhPD7ntH8sIs19Mp/V5CwWS34iIhOBnkATb8tSEQnDhH1NU9U+eZ7PCjmLxZKfiMhhYFF2AkxEJgA3qmpoXuezy1WLxZLfBAM7cuizAyjri8mskLNYLPnNcaBMDn1CgJO+mMwKOYvFkt9sAJrl0KcZsNEXk1khZ7FY8pvPgbMiUtnbSRGphPGZ+9wXk9mNB4vFUqCxmpzFYinQWCFnsVgKNFbIWSyWAo0VchaLpUBjhZzFYinQWCFnsVgKNFbIWSyWAo1NmmmxWK4IItIKeAi4DiiJKVizDvhSVVf4bB7rDGyxWPIbERkJPAV4E0ACfKiqz/hiLrtctVgs+YqIPAH8A9gJ/A2oBRRzHh/GZCB5SkT+4ZP5rCZnsVjyExHZjEm3FK6qp7ycLwlsBk6ranhe57OanMViyW9qApO8CTgAVT0DTMJodnnGCjmLxZLfnAAu5tDnotMvz1ghZ7FY8pvJwK0iEujtpIgUAm4FpvhiMivkLBZLfvMyJuvvAhFpJyIuABFxiUhbYCEmn9zLvpjMbjxYLJZ8RUT2AIWAik6TG7M0LUOa7+4hIDHjUFWtfqnzWWdgi8WS3yQDCUB0hvaYHMZdlkZmNTmLxVKgsTY5i8VSoLFCzmKxFGisTc5iseQrIlItt31Vda+I+AOVU55f8nzWJmexWPITEUnCBOHniKq6RKQu8CdGXl3y6tNqchaLJb8ZSy6FnMOZyxiTitXkLBZLgcZuPFgslgKNFXIWiyVfEZG/i8gpEemcof05EVnjHH/32Xx2uWqxWPITEfkJaAWEqmqS03Yf8DUQh9krKAzcpqrT8jqf1eQsFkt+EwEsTRFwDk8Ah4EqQCgmdtWmP7dYLNckFfCIWxWREkBzYIKqnlbVE8AvQH1fTGaFnMViyW/iSe8O0gojixZ5tB0HSvtiMivkLBZLfrMXaOPxvLfzuMyjLQQfZQa2zsAWiyW/+Q4YJiK/AEcwtVcXqepxjz6NMVW78owVchaLJb8ZA9wNdHOenwQGppwUkdKY5ewEX0xmXUgsFku+IyJ+wC1AELAkgxbn27mskLNYLAUZu/FgsViuOkTkJhF53SfXspqcxWK5EohIUaAU4Ofl9DPAs5hC1AqXl0sOrJCzWCz5jIjUBMYBLcg+fZJ6nr+cXHJgd1ctFkv+8x5wPbAY2AckeekT6Rxj8zqZ1eQsFku+IiLHgNmqel82fV4DBl+u9uaJ3XiwWCz5TSlge35NZoWcxWLJbwTvS9S/BGuTs1gs+c3DwB859JmKR6aSvGBtchaLJd8RkUaYHHLVgV3AcFXd85fMZYWcxWLJT0TkBmAhEODRfApooqr7RaQGsBl4QlXH5nU+a5OzWCz5zVsYH7heQAlMFpKSwIsAjka3Gejhi8msJmexWPIVETkBzPR0IRGR+UBFVa3vPP8SaKuqtfI6n9XkLBZLfuMPHMjQtgGo6vH8EFDRF5NZIWexWPKb9UB4hrYzQBGP52WAi76YzAo5i8WS3wwDOmeou5rqNyci/piEmht9MZkVchaLJb85C/wITBeRMSJyD8aVBBHpDczClCUc7YvJ7MaDxWLJV0QkCRP1kCJ8PP8GcANDVNUn+eRsxIPFYslvBuM9xVICZkNinqoe8tVkVpPLAyLSBRiJSfr3uaoOycO1vsT4BR1V1XCnbTDwKHDM6faSqs4QkWBgItAaY3KIAV4FXsYU5N3t9J8BvAB8DlyH+VH7RlX/61y/KSaVTRGn7zNqvxCWAobPbXIislBErsg/ii/nFpFOIrJcRE6JiIrI1Azn/YCPgK5AA+BeEWmQhynHAl28tL+nqpHOMcNp+wAIA24AimMKgqwHbgfOA28DdZzjLaCwqjYCmgIDRKS6c53RwN89+nqb32K5pvG5JiciCzFOfNll/PxL8NXcjhD4EzgNfI8xlG5V1R88+rQMCPBf3jjc7IQfOnwEgIoVyme+oMtbdufMJCQksHPnTho2bAjAwYMHcblcVKhQIbVPUlISmzZtonjx4tSsWTPd+IsXL7Jx40bq1atHUS5y8uw5jp8+h0uEWqEVSEpOZmv0AcKqV0ZV2bbvIKFJ5vM/l5xMPEq5LO41ODIiV6/hamXt2rXHVbVcXq5RRfw1nvT/L8dJnq2q9sfhKsba5LxzCxAIvKCq32XRp3KJEiVYs3QhAN9+9wOr1qxl1IhhmXsGlczVpNHR0fTo0YM1a9YAMHjwYMaOHUtiYiLNmjVj+PDh7N27l1tvvRV/f3+2b99OoUKFePbZZ3nllVdYs2YNHTt25Ntvv6XJxRiWRG3hnW9+okTRIOav3UBcfALvP9efR2/rxJo/d/LS6G95ctdhALZcTGRqQhwvFQv2em+3O/d0rSIil1UfwJMElLsomq5tNOfK5vW6/2uISLXc9r3cug6eWBcS71RyHg9m0yeTtig+1l0ff/xxdu3aRVRUFBUrVuSFF17A7XYTExPDhQsX2L17N3fffTdjxoxh3rx5eNPKz8TG4ufnYv+0L9g5aQzv/fAzuw8c9trXkj0C+IukOyyXxW5gTy6PPJMrISciD4nITyKyW0QuiMhZEVkmItmlLy4sIm+JyB4RSRCRXSLyuogUyqJ/mIh8KSLRTv+jIrJERB730vceEVnr3MtREflWRCp5u67HmM4iMkNEjnvczzDHiJ/Sp51j03vDaVrg2ONURNpluGTMxcQ0h+yYAwepVMEnUSiplC9fHj8/P1wuF48++iirV68mNDSU0qVL06lTJ8qWLcu9995LUFAQf/zxB6GhoSQkJKSOP3D0BGfOxdG5RRMC/P0JKR1Mq0ZhrN26i9CQMhw4eiK17wlNprTY37zsEIRCrvSH5bIYm8WRkkNOMVlKxvpistwuV0cDWzCFJw5hQi66Ad+KSD1VfdXLmB+B5sAkTHjGrZit42Yi0stzF09EumN2CwtjHAG/B4KBCGAQHk6BIvIcMAJjL/vGeewMLMeEhmTCyRf/BnASmA4cBRoD/wS6iUhLVT2LeYPfANoBbYGvSUvcF016fo9PSGBPdDSVK1Xih0k/8d1Xn3ub/rI5dOgQFSsawTllyhTCw8OpUKECtWvXZvXq1cTFxTF37lwSEhJo0KABFStWxM/Pj40bNxJZL5hvZy2geYM6LFi7kX5d2hIXn8Cqzdt5uk9PKpYtTfGgImw/FUcdP38WJcTTNbBIDnf0v43R5K70XVz7qOojWZ0TEcF4CTyFSa7pkwlzPIBaXtoKAfMwAqyyR/tCjCTeDpTyaA8EVjjn7vdoL4sRTomYTYOM84R6/F0d40tzEqju0e4CfnKurRnGt3falwPBGc495Jx7L0P7YKe9XXbvS+2aNbVO7Vpas0Z1fev1V1RjT3s/csE999yjFSpUUH9/f61cubJ+/vnnet9992l4eLg2atRIe/bsqQcPHlRV1XXr1mmNGjW0cOHCWrx4cX3qqad08uTJWrlyZQ0ICFA/Pz8tUriQPnFHVz01d7ze0b6lNqhRRetXD9UhTz6g7uVT1L18iq78YphWcflpeZdLuxQuopOCy+lPpUK8Htc6wBrNxXc9u6OC+OlLAcHpDl9c1x5eZc5y4AdfXCtPu6sicrsjXB5U1W+ctoUYLegBVf02Q/92wAJgoaq2d9peAN4FPlDVZ3KY72WMS8S/NYM3tFPLcQfgUo/dVRGZAtwGhKvqZi/XXIcR0iEebYOB14H2qrowq/tpdl0TTdl4yJZcbjz4kqQVU3PuBPzcfUCu+t1+8khebueKIyJrVbVZXq5RyeWvAwKKp2sbnHg6z9e1ZEZEhmNkSJ52xCGXy1URqQr8C+iASYeScV1T2cuwRV7almBCNpp4tN3gPM7Mxa1cl9W1VXW3iOwHMu7ctMRom3eJyF1erlkIKCciZVT1hJfzFguQtvFgyRdKA8V8caEchZyjIa3GlBFbAszBLC+TMMvHBzG2tIxk+ulX1SQnYV6IR3OK4T9jfilvpKhEWakVh8ks5MpgXmdOcXDFgL9EyM2aNYtnnnmGpKQkTpw4QWhoKH5+fvj7+7NmzRoGDhzIp59+SmxsLEFBQezbt4/g4GDWr1/PY489xvnz56levTrjx4+nRIkSnDhxgjvvvJPff/+dhx56iFGjRqXO1aVLFw4dOoT73ClaR9Tnwxf+jp9f7vz0LNljbXL5g4h0AvoAm3xxvdxspz2PERSPqGo7VX1aVV9V1cHA7GzGZfKKdaIEQkjvfnHaefSmDWYkZWPBi8ctABW8tJ0BTjlL2EUAqipejjz743gjKSmJJ598kpkzZ7JlyxZiY2MZM2YMUVFRqf5wHTt2ZOrUqaxevZrChQvz3//+F4D+/fszZMgQNm7cSO/evRk2zPjgBQYG8uabb/Luu+9mmu/HH39k/fr1rB83kmOnzzJp/vK/4mX9TyJAIZF0h+XSEZEFWRxLRCQas/kYgLGN55ncLFdrO48/eTnXNptxbYFvM7S1IbN/2UrgTkx41Kwc7uUPTOhSW2C+5wlH46ziZcxKoLuINMzh2peOyy9He9vqFSuoXaM6NSuUAfcFigYFMfPnqdzYJC1nYKfWLQDYE72XIoGFiYnejcaeZtu2bbRp2hiNPc0tN15Pl6FDefPNNylatCitW7dm586dmeYrUaIEAO4TR0mMi0PizsPprBXUa93Wlp+I2OWqj/AmBwCSMQVtZgDDVNWbyeuSyY0mF+08tvNsdBLe9c9m3KsiUsqjfyDwXy/9vsaETT0uIjdlPCkioR5Px2Psa095xF8iIi5MIj5vr+c95/EzjP0t4/WLOtWD/hIOHDhAldA0JdXf359PvvqKpje25dMvx2bqf/r0Gbp0ugWA8AZhTPvVmConTvmZ/Qdys6KHzp07U7Hv8xQvEsgdra1N3FcIUMiV/rBcOqrqr6p+Xo4AVQ1R1R6+EnCQOyH3Mca9Y6KIjBeRoSIyA7NRMCmbcX8Cm0XkA2enZBNmk+GkZydVPQ70xdj4FojIdBH5j4iMEpHFGDtgSt9oTEWfUsA6J+HeOxgNrykmTzwZrj/PGXMD0ALAeQ0fi8ivGPve4Fy8D5dFxt3rV18cyD133s7MKZP46JPPWLx0Weq5UZ98CgL9+twNwBcfj+LjTz+nWet2nDt3nkKFAsgNs2fPJmbccBIuupm//k/fvZj/cQSxEQ/XIDkKOVXdgPE1W45xAH4cU0bsdmCM0y0yJSICo4oClMNsWPQE/uHMNRjHmOgZEQFMxqQT2ojZef0ncBfGVy0lLVCYk47oacwObTFMBo0nnGu2wqi6mSIiMPnke5O2sfCsc/3KwKfAK54REcArTr/HPCMiLofQ0FD2x6RpYLGxcVSqUJGQkHL07tWD1WtMIfGvx33H/EVLCK1cGXH+ecLq1WX2tMmsWbqQe++6g1o1auR63sBCAfS8IYJfVkbl5fYtHgjgJ+kPy+UjIvUdJWi1iGx3Hj8Ukfq+nCdXLiSquhy4OYvTIiIX8B4R0RJ4Sz0iIkSkvfOnt4iICExEwqVERFTUtNJm7ZyIiO/JHBHRxrm3cqqabtmaTUREH6CRR0TEJdO8eXN27NrFnuhogkuWZPyEH5nwzVfExsYyZ94CXntxELPm/MY7743km8/G8NDf06LYjh49RkhIOZKTk3l76LsMeCR7B/Dz589z7tw5KlasiDspiZm/b6R1eN3LuW2LF0Swmw0+QkT6Y1aJGbf+U9KBPamqn/liLl9lIQlX1V2eDU6M6kzgRREZo6oZDUr1gYaqmqJ9vYxxFO4B3IezaSEiZYHvnHu9OeNa3dNm59jphmA0uuuc5S0i8n8YIXl7xht3hO4bmGiMbqp62uPcQ8BXzvnncvtmeOLv78+o4cPofOsdxMfHc/HiRfo+3B93UhJ9776TLp1uoXajJhw8dJiW7TvidrspHlKZ94cN4fz5WD7+zISK9e7Vg4fv75d63erVq3P27FkSExOZOnUqc+bMoUyZMvTq1YuEhASSzpygfUQYA7pltzdkuRSMJmeFXF4RkRaYUM0zwJsYOXEAs7LqArwGfCwiUar6e57ny0vEQ44XL6ARESk0a9ZM1+QmBVGc15DaTOTms5CiuVs9J838Ilf9/LpmGUZYoPBFxEMt/wB9p3jpdG13nT5qIx4uERGZiDFjNVHVTEZjEQkD1gHTVLVPXufziSZnIyIs/wsI4LKanC9oA/zsTcABqOpWEZkG3OiLyfK8Ce5oSH8Aj2EiDj7HaFNvYNxD4BIiIjCbAyU8mjs4j76KiMhIGYzj4etZHCkpzXMMMalevTqNGjUiMjKSZs3Mj/vAgQMJa9KcRs1bElyxKp17mRXzxMlTadjsBlzFSrHmj3UAzJr7G2FNmlO1XjjVwsIJb9aSxte3Ij4+HoAJkybTuHFjGjZsyKBBg1LnTUpKokmTJvTo0QOAqKgoWj33H5r+4w1aPP0mq7ft9rxNPpj6GxGPv0bjx17j/fffN/czcSINGzbE5XLhqZ2uXr2ayMhIIiMjiYiIYMqUKTm9DQUau/HgE4Ixq6rs2IFJ3pFnfOHp4+uIiDIYv7kU3M7jXx4RkfEgfYREriIiFixYkCmaYdPvK3j4/n5UqliBnbuNwAlvUJ/J333LTa1bAUZQ/eP5gUz78XuCS5YgsHAgP377FQtmTicgIIATJ04y6JXXmDdvHps3b+bIkSPMmzcPgJEjR1K/ftqG1KBBg3i1b0/Wjnqd1++/lRe/TPP02RR9gC9mL2bFey/zx0evM336dHbs2EF4eDiTJ0/mppvSuyqGh4ezZs0aoqKimDVrFgMGDMDtdvO/iAB+SLrDclkcx/yfZ0cIGdzNLhdfLFd9HRHhj1mPp3AW447yl0dEeLPJ5ZVOnTqxf/tmfp01hz533s6nX45FVQmr5+x6qrHFrfp9DbVr1mTX7j1ENGpEw/rGEfj/Bj4PwJ69e6lbuzZBf7+VWKDNvmP88OhsqtStzLS1OxlYL5RROw8S27sVyeu3cL7fPbhuup1zh5XKYQdw3XQHANuOT6VFuw4U69QXgLbrjjBlypR0mqEnQUFB5jb3buTC/hhEk9C9G1F/718dqdbIh+/e1YbYjQffsAHIyY7ZDONSlmd8oclFO4/tPBvzGBHxlUe/lKXnXx4RIVlkF85tRISI0KlTJ5o2bcqnn36a2v7coJd45603mDV3HuXKec8cc+DgIaqEVmb7zp2IwHc/TmLEqI8YOmIkALVr1mTr9h3sjY3Hnaz8cugkB+ISGbQxmrfCq6V7Ye80qs6gV16jar2GDHz5Nf7zxmup58Ib1GfJsuWcOHGSuLg4ZsyYwf79+7N9XatWrSK8Y28ad76D0W+9in8WAq6gI2KXqz7ic+CsiHhdnTn/h6edfnnGF9/WjzEZPCeKyE8Y21k4Ziv4R4yvmTeKAcccx9ckTNxaIPAr6TW8lBzjScBCEYnHbAbEYlxFFKgBJiJCRF4EhgPrReQQUBFTtk+dMekqkajqPBH5HrgXOCAiCcAu4HcgRS0ZTC7K9S1btoxKlSpx9OhROnbsSFhYGGfPnqVcubLMmvsbfi4/Klf0niI9ZWfV7U5i6YqVDHzmaTZs3syUX6bTtEkkHdq3ZfTI4Tz490dxAS1KF2flybNEBBejSXAxFh9L28H9fM8RRgx5hztu68WPP02h/xNPM3e6yS9XP6weg557hk69elOsaFEiIiJyFFotWrRg09wp/LlzNw+98Apd27UmMNCbmbVgk7JcteQNVZ2MCQDI6vxBsvbLvWTyrMnlMiLCG4eBzZg6of4Y4z/ABvXuS7ECiMcIQ8EIq2qYX4TUb56qjgCGYgRbPYzg3AKMwwjKdDjOw/cC5zAFNi5ifPjuB4Kcbq9kHOeNSpWMIhgSEkLv3r1ZvXo1y5YtY8Kkybw5ZBg7du1iweIl3P/I3zONDa1cif0xBwitXIm2rW/kzNmzVKtShW6dO/JH1HoAenbrysK2jZjfthF1ihdBVZhx6BQNZv/BQ2t2sOj4WR5Zs4Pv9h3j9lt7AnDX7bexeu0f6eZ65MH7WbtsEYvmzKB06dLUqVMnNy+P+rVrUrRIETZtz5wY4H8Fl6Q/LFc/PgkxVtXlqnqzqpZS1eKq2lpVp6rqQsdoP9ijbzugtqpep6oRqhqsqn4YgTIfGJiFGhuK8Vcr5vQvhrGnNcY4DwOpzsOPYTYs2qlqYVUNV9UHVbUkHna5DM7DjVW1lqoWxwjdqRgBiarm6AwXGxvLuXPnUv+eM2cO4eHhtG3blgrly7N/22Z++PpLbm7bhm+/+DTT+OZNr2PHrl2E1a3D+o2b+O7HiXTr3JFFS5bRoH49wERAAJxKdPPZ7sN81qw227s0ZUvn6xjbrA5ty5bgi2Z1qBBYiEVLTEzs/IWLqVMrfX3WlOvs27+fyZMnc++992b5uvbs2ZO60bA35iDbdkdTPTTbmkEFFluty7eISD0njOt3J6xrhYi8JyK5+9XNJVfEuJIxOsJpSxSRjzBqagdMSJYnb6ZERzj9451IhgXA30hb4j6I0SQ/8JbJQFVjPJ72wyx9h6RERzh9kkVkIMZJOFc/BEeOHKF3794AuN1u+vbtS5cuXahduzYJ8Rfo1Ks358+f56IjMKZMm87T//wXx44fp8cdfYhs3IhRw4dx9wMPc/rMGQTo97e/061zR7p36QzAMwNfZN08E4v6Yr1Q6hTzXnhmVJOa/POlV3C73QQGBvLJh++nO39nvwc4cfIUAQH+fPTRR5QqVYopU6bw1FNPcezYMbp3705kZCSzZ89m6dKlDBkyhABNwuUSPnrzZcqWLuV13v8FXHa56hNEpB/wJellUC3geuAJEXlMVb/yOvhS5/orIx6ynDRn5+GXVDUlMH8hZre0mqruy3AdP8wSNlZVg522iTj56VQ1291Yx4Z4O9BBVed7OR/tzOv1m53biAeNPZ1jH2e+HPvE9uuaq2sFjZuRuzlzGUGhe3O30XW17q76IuKhfkAhHVsmvXfSDUdibMTDJSImt+NazKblc5iY97OYyKaFwPtATaC1qq7K63z5rsnJtZ9O3fI/jI148AkvYFZIPVV1B6T+wMep6kwnlHIHMBCjsOSJK5H273Kdh+NFZKxnQxbOw75Mpx7qzDM2i/OZog2yihwA2LBpE61u7uQ1kiEsshmBpcsTUq02kTe0pkSFKrw/6mNe/fdbNL6+FZE3tKZTz95M2H+UJqR0EHAAACAASURBVL+to+GctYTP+YPyv6zi+fVpEQ1x7iTuWPEn9ZtcT3izlrz42mC2bd9Bk5ZtUo+SFavy/kejGfjyq4SFhdG4cWN69+7N6dPpNc6Mr+3V4aOI6HIHTbreRef7B3DwyFEAZi1aRtGw5pRseANNmzZl/vz5fPjhh9SrVy9TdMaGDRto2bIlDRs2pFGjRqnvwbWAiODKcFgui/bArBQBlxFVPQz4LKzrStRTnIVx5yju5dwY59xgj7aFTpsCYzP0b+e0L/Boe8FpG5mLe3nZ6fuGl3M1MdplpnlTjqZNm+rw4cP13nvv1e7du6uq6pYtW3Tr1q3atm1b/f333zWFixcvaqNGjTQqKkpVVY8fP65ut1uPHz+uVapU0aNHj+r521pq3yrldGqr+hpSOEC3dLpOD3Zvrudva6nnb2up74RX0+LFi+muTVF6MmaP1qpRXV97cZA+OeDR1BqvsccO6vwZ03QAxbQ/RbUCLu1KoA6gmA6gmD5KUS2CaF+CtBuBmnjqqCafO6kDn31aBz77tCafO5l6vPufN/Xeu+7Q7p07qarqmTNnUl/PyJEjdcCAAaqq+vzzz+ttt92m3bt31w2L52iZUqW0w0036oWY7Zp8bK8e3rJWk4/t1cRDu7RRgzBdt2CmJh+N1mNb1+nFQ7s0+Wi018OX4IP6qA0CCmlUlWrpjtxcl7Sq8N6Ow1mMaYVJA34SiMM40D4L+GUzz4OYVdJ5zA/4QqBHNv2LYMIvt2HMPkcxbl/18/pe5fB+XADezdCWBLzu8XwIkOCL+a6EJhftPLbzbMyF8zB4pC/Pxnn4r06nnkpiYiK//vor/fun3Xb9+vWpV69epr5z5syhcePGREREAFCmTBn8/PzYvXs3devWTXUSbl+uJGN2HaJm0UCqBhWmRECaRWH7+XhKlihJzRrVKVWqFI88+ACb/9yabp6goCDatzUv2w+hLH7EkmZ3PUASJRCK46IK/qk+cjc0b8aBgwdT+8UcOMCM2XN55MH7U9tS6keA2UEWEWJiYoiKiuKpp54CoGFYXc6dP8/zj/encGFjdQgpZ0IQ5yxYTOMGYUSEm3DgMqVLXVuVxARcLkl3XAJnMAIl45GpGpGI3IqxU90ETAE+wnz33wN+8HprIu8CYzF+oZ9hXKYaAb+IyD+89C8MzMWkNToLjAR+wySXXeOkQ/qriCXNZSwr6pC1GemSuBK7q5frPAxwm4h8QFqSzVpkcB5W1eMi0heTjHOBiMzE/AqWwLibVMG78/A6EZmA+TJ2xtj2/sT4zHll//79fPXVV6muI9mxfft2RITOnTtz7Ngx7rnnHgYNGkTt2rXZunUr0dHRlHIiGTaeieO5umluGoO37OP7/cbto22vW1PbQytX4reFi6hQwftqOwFlL24aeezr7MJNbS8f+1ffjufuO3qnPn/uXy/xzpuDOXf+fLp+L7/8Mt988w0lS5ZkwYIFPP744wwdOjT1PfjplxkUKlSIlWujeGPYSAILF2bYGy/TvEkE23ftQUToctf9HDtxgj639WTQU4/l+N5dTVyiYPPktHq4UmWFiJTACKkkjAvUGqf9VYyL1Z0ico+q/uAxphVmBbMLaK5pORqHYQz874rIdPXwIMCYjW7E/J/0UdVkZ8wEjPvUlyLSKKXdxxzEbDh6xUm51gvIXb6wHMh3TU4v33kYjDHydozaXgMj7EoBHTPM8Ssm9m0pcAvwf8CTmO3pYyLSzKPvCIwLSiIwABjkXHss6W196RCRHgEBATRt2jQXr9q4lSxdupTx48ezdOlSpkyZwrx58yhVqhSjR4+mT58+dFqyidAihTiacJHeldLilwc3qMq2zk1pVboEf27dlv4+spgvGWUe8YQTQAnnY05yhF7NDELu7WHD8ff3p18fk2lq+szZlCtXjqZNIjNd9+2332b//v3069ePp59+mpCQkNT34Ny5c7z45hAqlg/h1OkzrJg1laGDX6JP/ydQVdxJbpau+p1xY0ay5JdJTJ0xm3mLl2Wa42pFBPz9JN3xF3AnJlb7B/Xwz1TVeNKc0h/PMCbll+JtTe9mFY3RAgtjFAvApPL2GDPIU5Cp6s+YDcEGZB97nheWAW1EJKM211FEJmM0yqMYTTfPXCk/uWzTqWfo285ZSu7BbCpUxQjIpRjVvA8wU0T6quoESP0QB2LU/eMY1f0YZiOhPSb7cMovZGHMlyYEWI/Z/AgGnsEsGbLixtOnT1O9enXi4+M5e/Ys9913H+PGjfPaOTQ0lLZt21K2rFm6devWjT/++IMOHTrQs2dPevbsSWzvVjy7fjcVAwtRPjBTYTF6Vy7No5vTXAxjDhwkuKT3koiLSaAkLhp7FCjbTxJl8SPI47ft6/Hf8+vM2fw2fWqqC8uylav4ZcZMZs6ZS3x8AmfPncv02vr27cuIESMoXLgwM2bMIC4ujuPHj3NL29a4xMXtPbogIlx/XSQul4vjJ04SWqkibVveQNkypUGVrre0548Nm+hwk2/sy381guC6fMFWWETuw3x/YzGri8Vq0ot5kvJ/4c39aTHGPtdKRAqrakIuxswEXnX6pCSSreXcx3ZV3ZPFmDbOmAU5vbDLYAJGgLbA/B+DsU+2xGiwPwPPqOohX0x2rRVVuwn4XFXbqOr/qepDmA8jGRjjqPoAj2KMsL9joiv6O/3vx7iEfOJxzRcwtSYmY1Km/0tVB2AScLYhC1T1/xo3bkx0dDQ//PADN998c5YCDkyZwA0bNhAXF4fb7WbRokU0aGBsU0ePml3KU4luJu4/zt9qpHnE7Dx/IfXv/XGJqCp7oqNJTEzkh0k/0SSicaa5XnnjLRJRWmWowLgTN7U8ftf24WboeyP5ecJ3qdlGAP77xmvs37aZPZvX8/3Yz7n5pjaMGzeOHTvSNsOmTZtG27ZtU21yxYsXp0mTJsyeOI5bu3Vi/hJT1Hr7rt0kJl6kbJnSdG7flg1b/iQu7gJut5vFy1fRoJ5Pndv/WgT8/FzpDqCsiKzxODLH7BkqYMwqb2P8wOYDO0Qko7aUYtDdnvECqurG/Nj7YzbGEJGiGE+C81kIhZQPzbPYR5ZzZDPGZ6iJhGqgqks9mttjVl+lVPUOTe+0nyeutXQSZ4B/ezao6hoRGY8Rar0xGw9POacHqOqZDP2TMAVtUngYIyQzqu17HPtfujTqOZFV5ECpUqV4/vnnad68OSJCt27d6N69OwDPPPMM69evx71/D25VHq2Rlvbutc372HH+Ai4RqhYpzOj3R9D51jtISkri5MlTDP9wFImJF5n6y6/MmTaZEiWK8/bQdwlG+AkjIBsSQG38icFNGw8XxGUkUOL8eTrdahJ5tmjejDEjR2T52l588UW2bduGy+WiWrVqjBljrAujRo3iwIEDHDt2jCbtupKcrNStVZ1GbTpSKCCAsaOGIyKUCi7Jc4/35/pOPRERunZoT/eOPovD/ssRwC+zTe645uwM/BVmCbgZEyNdE1PB7u+YVUhLVV3v9E1RzbPKmZ/SnuIPeqn9L3fMX4qqLsm51+VxRSIeLhWP5eoCVc30XyFpBWdGYtxCzgNHVNVbkkzPccUxdrf9qprJECppNSe+drTGdOS6xkMuie3dKsc+RcfnJgs8PFY0NOdOwOhz+3LuBEix3IVy6fHcXY9cfO+knO/8sH0R8dA4KFB/rVM9XVvVDdsu+7rOjugLwFRV7e20bcfsLNZR1UyZEERkOWZZ11JVVzppiQ4AB1Q104fu2L0SMe4YgU5bX4xnwXhNq3TnOaYTxmwzW1VzzL5ztXOtaXI5pTUvSQ4RD2Jqt/bAGDZTPkA/ETmAsdvhnIskzbG4m4jsxNgA9zjHsMDAQGrUqMHFixc5ffo0Dz30EKNGjUqdq0uXLuzZs4e9e/dSrFgxQkND2b17N4MHD2bo0KHExcVRs2ZNihQpQkJCAu5d21Juhh3nL9CvajlGRJjg+jh3Evf/vp0VFasSGxtLqeBgju0z9rmBL73KLzNnUSgggFo1a/DVmI/Yh5vlJOBGAaEQggC9KUIyMM3R8la2uomYAwfp1b0re6L3cvjIEVwuF48+/CDPPJG28/nuu+8ycOBAwsPDqVatGtOnT+fVV1/l559/xuVyERISQunSpZk/7zfKlSlNs4jG/LFxEydPncbtdnPsxElWzZnGwNff5tCRo2hyMvsPHaJWtWr4+bl4feCz9O5uPo7vJ//Mf0d+jPgXolKlSowbNy7VlnklERH8/H262TAGI+Q8XZ1StCjvxta00gBnMjxm1d+b1napc1zTXGs2uZzSmncEUtbyWUU8jCVNuKV8iMWA91Q1UlUjMTUqupCWsWQFJlPJPZhar1OA76tVq8bGjRsJDQ3loYceyjTRjz/+yLZt27hw4QJt27Zl0KBBBAUFcfLkSdq0acNNN91EVFQUp06dokKFCqy4OYIX64WSpMrb4Zm1mKdrV+KXid+zbN5s4uLimDl7rnnRN7dn0+8r2LB6OXVr1+btYcNZRgJdCKQQgh/QkUB6UgQXUAjhToK4kyDWLV9MtapV6NGlM+/+5022rF3Fivlz+PjTL9iy1fjg7Y+JYe7cuQQHB1O3bpqZZuDAgWzYsIGoqCh69Ohh7umHrzl77hwJiYlsWDyHaeO+xN/fn+uvi0gdN27MSNYtmMnZ3ZvZuGQOMyd8w2MDX8btduN2u3n2lX8zf/L3bNiwgcaNG6f74bjS+Lkk3ZFHjjqPnjkOU7bPM9nDRMQfs/PvxqQFQ1VjMT/oxUTEW7LCFKOnp/0tyzmyGXPNcq1pcteJSHFVzeiY1s55TGnfBISLSBNV9UyljqouTnH8VdVzjoZWC4+c8x59Ura4T2EMo9NVVUVkPVC2mLgpsvE3nup6I5/9PJuGNauStGJq6lxFMVtFF48fISFmJxt/Hk9oiUCWT/+J/+vTnfenzCVp5hcknTtF4qFogpYuhx8nET5zFiU63Iz/H+sIGjEMMHmougKPF6vCOZJJ4gJDbr+Pn51bTMlAuAc3W7jIDR3a88ATj+H/4yQa1g8DSE2l7sm2eb9w9NBBbm0Yigjo7vUUA8KqlCdm7TLqF0rguRde54WnnuHWRYu4787efPHNOIg7Qwl/Usstxp46QeXy5ShTqyG4/Il1K0nBlageWZSixUvgH1gUCa4AAYWR4Aok9e9HypbK8QsJcPYM5+67CwWSz5zhSP/7Kf3ZZ5zZt4OaDeqRvHWlt+8DAK6wXCVuzjMi4Odbt5GWzqNnpaH5mOw4XTBF0j25CfNVWOyxs5oy5n5nTMbMHV09+qSwC9gH1BWRGl52WL2NuWa51jS5khgP7VQcn7d+GK3suNP8gfP4iYiUzNDfhfFDSuErjE35aRHZICJfiknLHkpade+OmPcq5QvUz+MclUPKcOpceqfZFLo++wYV+z5P8SKBHDp5GndSMkP+dicul7At5jDVHxgIwJ/7D3lNV+6NBJSLQOVMxcdhKxcJxpVtKnVPfpg5j7s7t0+XASX6wGGitu6kRaP6TFu4jMohZfn0q68JLhWMy5X+K/Py4DepUrch4ydM5N+vvASYyIiiRYtSsWJFqlatyj//+c902YcffvhhWq/8k2e27KXFii20Wvkn74VVxd8lBLiEEWFVaLXyT0Jv6smfO6N55I6e2b4f+cmlRjyISEMRKe2lvRqQoqJ6bstPwnyP7/H053QifN5yno7OcLkU/9KXJX1JgeoY/9AEPISfGkN8ypihzv9EyphbMV4FW/BeNvSa41oTcouB/iKyWET+6wTOL8G8jgGkZf79HJOPrjlmm/4zEfmP0z8aI6RSGI4pnBOE2WUNwyxPf4XUeKjFQGtMxMV2jJKWwRPc+xd+5vuvEzNuOBcSEvlxyRoialahqWO8rhdagehvhlG2ZDFa1a/Fvm2bGTHkbfo/8XSWb0AyynISKIykOvmm8AeJuIDyTntKKvXH+/+NO27txZRfpjNvQfrv7YTZC7ina9pezvm4C9z1wmuMGPgk/n5+/PezcbSKaEhIuXIULpTZd+/twa+yf/tm+vW5i1GfmGSgFy5cwM/Pj4MHD7Jnzx6GDx+eGog/fvx4Nm7cyMxmddkbn8iz1cqz4Pp6jIg+THxSMheTlS9ijrO4RX1iFv9Co3q1GPJpxtSCVwYRcPlJuiMX3AUcFJGZIvKxiLwjIpOArZgiUDPwCO1S1bMYFyg/TLr/z0VkKBCF0fwmYfzM8BizHBiBWZFsEJN48iOML2hp4J8Zoh1w+i/HOB+vEpEhIvKdc/044G9/UbRDvnPVCzln5zTll+tGzBKyIfA0cDemQle3FEdgh0LAfsxGQllMRMNAjAPiEkzMHgCO2v8hxj4Sgfki1cH88qXUl0hQ1SiMQPwMI0BTnTgPHD1BqeLpSkekI7BQAFVDSuPvcjFv3Z/Ueuhf9HvnUxZs2MoDwz5j054D/LnfeLV4S1fuyWISKI6LwAxCdRsX2YubmwmkGK4cU6kDrN+wEbc7iaYNjNvUxYtu7nz+Nfp2u4Xbb7mJXTEH2XPgMAPeHM6nX45l77793N73fuYtXMR9f0vvDta3z538NPUXAM6cOUOXLl0ICAggJCSEG2+8MTXsq3JlYyot7u/HXRVKs/ZsHPWKFqGon4stsRfYcC4OgJpBhRER7urSgeXrfFK0yQcI4u+X7sgFCzA23BpAX0w4VVuME+yDmAD6RM8BqjrV6bMYuAPjEnXRGXuPenGJUNUXgIcwm3B/Bx7AuKz0VNVMRk3ne38LxiUrGJPXrSMmpKu5+iCPm68QkaCce2XNVS/kMAIuBOP/Ngz4DpPFIAgYrqo3qmrGFE0/YgTbD5ig5hQHyk2YVOmppQfF1Hj4EiMYP8XYIc5ifuFSvkxPOEbdtph41gHA4fMX4lFVvp21gMi66VOMn4+7wKHjpmykOymJqcvX0eOGCPZ+O4xdY9/hnUfuon3jML4Z+ChBgYUo7QhJb+nKU0hx8m2aIbZ5H26iSKQLRQhACMGVYyp1gO8n/pSqxakq/QcPpX7Najz3wN0ANKpTk8MLp3Bq2a+4z56gWtUqTP7uWzq0a8u4Lz9lx8606Itpv84kzHHsDQgIYP78+agqsbGxrFy5kqCgINxuN8ePG4vCzth4Zhw7Tf1igey7kMCO2ASqBRamUmAA22IvcDzR/L78tnw19WtV9/p+5DsC4udKd+SEqi5S1XtVNUxNqv8AVS2nqh1V9RtvAssZt0xVu6kpKVBEVRup6nteIiQ8x3ytqs1VtaiaMgRtVXV6Nv0vqOrrqlpHTZmAcqp6l6puydX7kX/8U0SO5tzNO1e9n5yI1NIM6dJFpBAm9OQmoLqqHnDaF2IE0Q6ghaYFKgdiflFvwNSFqI7R8I47j8kYQZiIWc4mY4L0A50pD2J+EAo7YyYDkwILBfxeqWxpjp0+i7+fH4luN8HFijLz/dcpU7I4t/7zbRIuXuRiQgJ7jhxn37fvUqZEMQDaDXyHqN37qFG+LMWKBHLuwgVcxUoRGBjI/pgY4hMSSEy8SHDJksz++SdKlChO1XrhBGC21hQjlVtSmHUkkgSp2l0ILt796Tue/df/paZSL1e2HN06d2To22m+1DUbRjD9/X8TVqMqS//YSNuHn6ZRnZqptqa3nupPtzbGqO+q3YTq9RsxcugQPhv7DdN/msAdfe9n2/aduFxCtapV8HP5sWrtHxw/fpyAgACCg40N7/z581y4cIGSJUsSHx9P9erVOb5tK/HJyVQOLIRL4F81KtIjxHj/fBFzjDH7jlI4NJSqlSrw1X9epUyprLwdcrfx4As/uetKFtVFLRukaysxe43NDPwXIyKvYdKvXZZSdtULuawQkdsxBa0fVNVvnLaFGCH3gKp+m6F/O4ygW6iq7Z22lzHG3H+r6usZ+tfECEuXZpX+vH5tXfVlpkw5mTl9IlevKaUAdHY8XsxbfezMjInNXVRM8s51OXfCCLlcEZS1MPLkbMecExsU+/CjXF0rP4Xc4tbh6dqKz1hthdwlIiIPY7KjbM5OM/Xonychd9W7kEjO9SC8+cN52xVaglGCPP9br8uqv6ruFpH92PTnlhREcrVEteTI55idungR2YAReGsxGyW5EnyXwlUt5OTy60HsFZF0oVjqUQ/C2d0KJ83ZN0JEPsTkjrte01LcnAKqicgpjEH3G00rsNM0sHAh6t31OF1bNuW95x7JVSGarOj/3lf8+sCLhJQry6Tx33DPg3/j4KHDnD13losX3TSoH0bbNjcygVhcCIEYS7QbvEcy3NCamIMHua/P3bw/bEiW885atprn3hlFUnISj/Tuzr8e6ZvtfcbHx3NTp24mQiMpiTtv68Ubr7zEyZOn6PPAw0Tvj6F69eqEhIQwd+5cQkJC2LRpEwB9+vRh2zbjh3pq25+U9Pdj6Q1Zpuu7+hCQACvkfMAnmGD8RhgPiOakuSd4E3x5klNX+yd2ufUgMiFp9SDcmPzyYR7XOIXJU5cxtVJKPOvPQFNgQIojMTC6WoVybP3xY3bEHGTWyqx3RHPDA7fcyMypkwCoV7cO61Ys4dsvPuH0wX2UK1uWG2+4nqNHj3EXQdxBEU6STClc3E2Q10iGqJVLqValSmqRaW8kJSXx1H9G8uvHQ9g0ZSw/zJrHll3R2d5n4cKFmT9jGutXLSNqxRJmzZ3HytW/M2T4e3Ro15YdO3bQoUMHVJVZs9Jn/pkwYQJRUVFERUXRKySYniH5Fv/tGxxN7lI2HiyZUdUnVPV6TAH4ppjd4E8wAk0wQu8xjMYXRS6Lu2fFVa3JYfyIwNjeMnKpCf3aYF5vPGkZR9diMpfUUNWxnpqYo0WmOFYKZpmcCJx1dlpLuGOOMK3HY9RPjOeDP94lsWhaenBPbj+Zcxbn9l0hOjoaXH6pZQI797qdOXPmUKt2bTp378WkSZP45Nx+Zsyey7b3RhJauRKjvRSqhrRIhhvLF8nS7rZy/WZq1wujdjuTDeWevvfxy4Y9hHfunbmzY2sToJhzfxfj4riYlIwUKc7PM2axcOFCNPY0D9x1G+2//IJSgf6QnJSpJKOqMjUpgPm/TKNE7VpZvym5tO/lFyLk1m3EkgvUpI6Kco4vIDV0LRwj/Jo5j43wvmLLFVf7T1G089jOs1HyVg/iIPCVmLJnjci+xkOK1OuDCYN5V1VPYuyAqZb9MuLi5F/kN/nDDz9w77338uWXX9K1q4m22b5zJ7v27GHDpi00bd2Ooe99kHmcl0iGjBw4epwqoWkmzdDKlThwKOc8hUlJSURGRhISEkLHjh1p0aIFR44coWJFEzpZsUIFjh47luX4JcuWUz6kHHWyE3BXJQJ+rvSHxaeoqltVo1T1C1V93NH4/pOXa17tn9LHGO1pooiMF5GhIjID4z4yKYext4nIt05saizGfeQsxpY3WlWbYLKaLMRobFsx2t5KjMd3e2Cvc60JGGfOFxwNL19q0SUmJjJt2jQOHDhg0pP3M4Eas3+bz6lTp5n361SWzJnB1F+mM29h9pEM3vC2s54bs6Kfnx9RUVHExMSwevXqVJtbbvl+4k/ce1fOO8lXHZfhJ2fxCXn6f7uqPyXNWz2IYxgP85oY4RaF0e5cGEEHRlAmYQKWCzvn3BitzURIpd6KHsXkpm+G0eJSc3ed0GRKi+/fypkzZ1KhQgUWLlzI+PHjERG+Hv89u3bvoXevHpQrW5agoCC6Zoxk2LgpXSRDVoSWL8f+mLSMVDEHDlKpgrdEFt4JDg6mXbt2zJo1i/Lly3PI0QIPHT5MSLlyXse43W6mTJtOnztvz/U8Vw1yWREPlivMVS3kwMTlqerNjud3cVVtrapT1aRQFvWogKSq7XAqcWE2DUaoqssZ2wSjqSnwqZMqvQMmDrUWxsF3EdBaVetiXE1SVR0nzfQNwFY1aabPxTt1HRclxNO80GWbDLJkxIgRnDx5kmnTphEUFMSsWbMY+t5Ifv1pAlu2bktNpb546XIahIWljvOMZMiO5g3D2LFrV7p06r26d812zLFjx1KLUF+4cIHffvuNsLAwevXqxddffw2YuhFZXee3BQsJq1uH0Mq5qf199WE1uWuPgvwpeU2VDkzHaIObMYkxK2MiHFph4lZ/FZHZTvqZ952ht2LqRXzlaJcAjx9NTuLJsyco7+fHdf6Zg9cvhXvvvZeWLVuybds2QkNDGT16NEuXLsXlctGxY0ciIyPp06cP586f5+4HHubEiZNUqRdOk1Y30SSiMd27dEq91sQpU7mna4cc5/T392PU8GF0vvUO6l93PXff0ZuGDbJ36Th06BDt27encePGNG/enI4dO9KjRw9efPFF5s6dS92Ipvw2fyG7du+h1c2d2LZjJ1XqNuSLr41v9oRJk7nnWlyqgl2uXiJOUow7Pe3dV+Q+rtWIh6y4xFTpr5KH9Oe1/QN0aIlMWXQykZvd1dyi50/l3AmTFy43+DqSIeNOalbkyqfQh7urvoh4aFo+WFf1SV+vPODDX2zEQxY4GXtqYWxqJzHZflL839aq6u5shntepyrQRVW9uxLkwNXuQpIXcpMqvWQu+1osOOW6rvRNXDOoal3HLHQdxhWkEybduwAuETkNRGWhjAzAlG1cr6r7MMkzLouCrG/nlCr9DGnpz7P6Jc62EE52TI+P49kzJ3jmzAnef9+seidOnEjDhg1xuVx4FsBZvXo1kZGRREZGEhERwZQpU3I1x/6YGG7u1osGTVsQ3rwlIz9Ovxcz/OsJ+EW05/ipvKfqnzVrFvXq1aN27doMGZJ1BMWlsm37DiJvaJ16lKhQhfdHfcyr/36Lxo0bExkZSadOnTh48KDP5rxsBCPkPA9LtqjqWVVdiPGISPF7vRPj5/oGWdRiwaRSW4zxS90mIj+KyEsi0l1ELsmgW5A1uZxSpa9Tk/78AlDEW7YTMvjn5ZZ9SW5+S7jAOyVK4w98Mn063bt3Jzw8nMmTJzNgwIB0/cPDQzvglAAAIABJREFUw1mzZg3+/v4cOnSIiIgIevbsmS6brjf8/f159z9vcl1kBOfOnaNZm5vpeHM7GoSFsf/wUeauWEPVilnJ+tyTlJTEk08+ydy5cwkNDaV58+b06tUrtW5sXqhXtw5RK5emzlO5dn169+pBqeBg3hxiUr9/8MEH/Pvf/04tgXjlsJrc5SAi3TEp0qYDH2G0sqaaTRlCVW3o+Lc2wtjOI4HuwP9hNgtdTtbvt1X1H9nNX5CFXEqq9IEpDRlSpaeoS4cxO7LviMjd6mRDFZEamF+TLAmOjOB2LyUJJ06cSKfZs7n3888B2Prmm0yZMoVBgwZ5vU5KYWfdu5EL+2MQTUL3bkS9CDmp1ij170q1S1GpttkoKFGsFPUbNuTgqfM0LFaKF0aPZ9iI97i1T19cNRrhKlsm07XM5DnbvVavWEHtGtWpWaEMuC9wz+238fPECTTwUjMiJVrjcpjnRHdUCzOvMcX+eP7kcXAnZmmPzG25xDwjgKsgL358j4i8AAzBCKPBTlsj4CcRaa0ZEoZ6oqrxmA2/3zNcM6XQzkVysdoqyEIuJVV6C4x/W0VM5IILU3T6rNNvP0bI3QH8ISKzMQKyj3ONXpc6cXh4OC+//DInTpygSJEizJgxg2bNsrdNr1q1ir/d34+9Bw7yzYj/5KjFZSQ6Opp169bRokULpk2bRuWKFYlo3CjngbngwIEDmSIjVq1Z65Nre5IS3ZHCy2+8xbff/0DJEiWY/+s0n893yYjAJX4u/8uIyBcYX9UHVNWzKM9jGN/Xj8kQuSQiq4BZwB5VHevtuqq6w3mMwyx9s6Ug/yxVwxSZvhF4EVPNaD+ZU6WnbC9/gHEcHogJGE7C7NJeMvXr1+df//oXHTt2pEuXLkREROQotFq0aMGmuVNYPe17hoz+gvj4hGz7e3L+/P+3d+ZxVs/7H3++pymhRUnrVCOlpmna0yKtWsWlFClbXW64LtkuQriEFCLL5Vrjxk+FoqKLtFBKhUSlmjbRJqppm+nz++P9+c75zplzZs6cOTNTM5/X4/F5nHM+38/3s3yX93l/3us+BgwYwFNPPUV8fDwPP/wwD957dzRTD4loPSPyAs+7Y+DAgZl1D4++h00/reSyQQOZ+OJLsR0wCgiClCqVpTjkiH7AuUEEzgu9fhHwFxG5Juic21Fj/1diNYliR+SMJuw4COwCPgXGonkZtqG5GzqGOTURJYoT0ZDpe4CbUDnC1Xmdx/Dhw1m2bBnz5s2jcuXKNGjQIPeTgKT69Tj5xBNZuSZb8vSQOHLkCAMGDGDIkCH079+fdevWsWHDBpq160hiUgpbtv5Cy7M78+uv0ZuxJCQk5MszIhLMmjWLli1bUq1adhniZYMuZtoHM2I6XlRwioe8op0xZmGoA1ZjOpiALapXPw/NtZIYq0kUV967SQ4h0+8UkRe8kOk+JAHJvpDpo1AbuX5oXohJ5AHbt2+natWqbNq0iWnTpvHVV1+FbbthwwZq165NKWDjll9YvT6VxISauY5hjGH48OEkJSVxyy0qH0tJSWH79u2Z+VATk1JYOn8uVcLJ5CJAmzZtMj0jatWsydtTpvLfV/8TdX+hMHny5Cxb1bVr11K/RhUAps+cRaMzI/uTKFg4xUNeYLLncw0+/j8ReShE/VHUtTImKJZELoSWFGPMYZumrRvqzhWc5+5fHoGz7Q+KyF0ooRtGHoncgAED2LVrF6VLl+bZZ5+lUqVKvPfee9x4443s2LGD8847j+bNm/Pxxx+zYMECHn30UUqbDOLihGf/NYoqlXMXpi9cuJBJkyaRkpJC8+bNARgzZgx9+/bNy1RzRXx8fKZnREZGBsOuGJqrZ0RekJaWxpw5c/j3v/+dWXfnnXey+sdVxMXFUbd2bZ6fMD5m40UNp3iIObwgtAWJYufxABGFTL/bF+F3LmqzU9ey0P5+SqFb3/3GmGxqw9atW5ulIbSr0cJszD31nl+7miPSIrSNi9SrINb9RYBIvTsi0a7GwuOhde2qZvFNg7LUxd/+rPN4iBHs7mm81arGDMXub8mGQlqGanB+RaOLPoQaHr5um4Xyps8mtLKx5nehvq4FgoyMDFq0aEG/fv0AuHf8RJr1HkCLPgPpdfnf+OU3zcQ2Z/5XtO53CSkpKbRq1YrPPvss4jEOHjzIWZ260azt2SS3bsfohzQ81+7dv9OjRw8aNGhAjx49+P33yIiKh81bttC1Tz+SWp5Fcut2jLhxZBaD4XHjxiEimWkIPWzatIly5coxblzWJECbN2+ma9euJCUlkZycnGnc/O5779OkTXtKVTiVpcsCAUAPHz7MsBE3kJKSQrNmzZg7d26e5p93iJPJFSxGAPk3vgxCcdyueiHTrw5WQYvIYDQvRChUI0gOYDm5aviikcQaEyZMICkpiT//VIuW26+9in/dqraNT7/6Fg9O+DcvjLmXKpVOYfrLz1DrrO6sXLmSXr16sXVrOGPxrPBClpcrV44jR47Q8dze9OnZg2kfzKB79+7ceeedPProozz66KM89thjEc89vlQ848c8RMsWzdmzZw9V69Znxocf0rVrV5o1a0blypWpUyebSzAjR47MDACapb/4eMaPH0/Lli3Zu3cvrVq0oEe3LjRJSmLqW28w4qasdnkvvaYSh++//57t27fTp08flixZQlxBbSmdCUlBYwbwoYhMRW3jVgCrbAThqFHsODmiD5ke6tg59jNye448YMuWLXz00Uf89a8BU6EK5ctlft+fdiDTVKNFkyRqVqsKQHJyMgcPHuTQocimJSKUK6f9HjlyhCNHjiAifPDRTK68Umn+lVdeyfvvv5+n+deoUZ2WLVQW+OPqNZl5VsuUKUOpUqVo1apVNkf8999/n3r16pGcnByivxq0bKkJ1MqXL09SwzPZ+ss2kho1pGEIxcOqn1bTrYvetqpVq3LKKacQS/FBSMTFZS0OscSNaD6HGsDDqEP/PhFZJiKvishNItLZejpEjOL4t5RqP7ug/wxARCHT7xWRD4MSUntC0X2xnybcfPPNjB07lr171fPMk7eNGjWKN954g4oVK/L5558jvgCUZucmpkz/iBZNkiiz9zeyOa1ZSJWsHFRGRgatWrXi559/5oYbbqBt5+78tn0H1SuciNm/h+oVTmT79t/CRhHJzZNh+Y9rSTtwgLYpjZj+7mRqVa/K0cMHwRyFtD8hrTT79+/nsUfGMGfGezz+xJOYUgazJ7RpS+qmTSz/fiXtupyLlLPSglLxyEkVMmVwzVufxfSP53DpReezectWvlm6lE0/raRNw8Rs/bVq3iz3RK+5QZx2tSBhxUOv2IKIVCbg0tUcGA40BEqjbl0VgDuMMTkmuimOf0XRhkz/EfhBRJ4WkfHASjRI5gE0fHpM8eGHH1K1alVatcr+7j388MNs3ryZIUOGMHHixCzHfvhpDXf+61FeGJc3pVR+Q5bnhH379jFu3DjatWunxshjx9P/gguyGQyPfugRRv79+kyuMnx/+7n4ir/y1FNPUaFCeHHosGHD1Je2S09G3nUvHdq2Ib4go/WKQHzprMWhwGCM2W2M+cwY84Qx5gpjTFMCGb5ADfZzVfoUO07OGPOdiHRFlQ190TV+i1pR70HdtULhczRw5t9RY4E/UU4u50SkUWLhwoVMnz6dmTNncvDgQf7880+GDh3Km2++mdnmsssu47zzzuOBBx4AdHvb/8preX3iE5xxenQ5r7OFLP/1V2pUr55jyPKc4Bkj9+vXj9WrV7Nu/QY2pG7kzvvuB2Dvvn20PLszX3/xKYuXfsOU9z/gjnvuY8+eP4iLi6PsCSfw92uHZ+nv4iuHc9nA/vTvn3OI9Pj4eJ588knMH6qcObvneTQ4o16e15AnxDlOrihhfV2X2+/7gd65nVMcObloQqY/jXo5HEE9Hp5Ac0T0wJf1K5Z45JFH2LJlC6mpqbz99tt069aNN998k7Vr12a2mT59Oo1sWPM9e/Zw3nnnMeaeOzi7bZs8jZVjyPK31OMmp5Dl4eA3Rn7iiSdYu3Yt5cqdzJa1q6hTO4EF/5tNQq2aLFv4BdWrV2P+nFmk/vg9qT9+z03XXcNdt/wjC4EzxvDXG0fS6MwG3HLDiFzHT0tLY/9+ZbLnfP4F8aXiadwo57wW+YKnePAXh2MeJf4uiUgHNNrIOuAsm3LQ7/FQg0DWrgLHnXfeyerVq9UItm7dzPBCEydO5Oeff+ah8c/w0PhnAPj43UlUPa1Krn1u27aNK6+8koyMDI4ePcqgQYPo168f7du3Z9CA/rzyxpvUSUjg/ya9lqe5+o2RvQAEHc/tTdmyZaMyGF646GsmvTOFlMZJtDinO5SKZ8yYMRw6dCikEfX27dvp1asXcRhq1ajOG/+emPsg+YGTyR2XKJbGwHmBiLyEKiSGGWNeDTrWBSV0G40xicHnxtoYOBKYnZF5uwQrHsL2F2m48khDKEVoNGwO527vKadEFgvP267mhDZderJ0+Yp8hRVoXS/BLH44a/St+Mv+6YyBj3GUeE4ODc0MmqkrGPPRFIUODs5O7jhFsZTJ5RFh8zz4PB4KDMEeDxkZGSQkJFCuXDmSk5Pp0KFDZhjw9u3b07r7eTTpeC4p5/Tg8utuommnnjTu0I1Hnno2s8+uf7mEhg0bZoZUHzx4MFWrVqVJkyaZbS655BKaN29Oi/bncHrjprRof06WeW3esoVufc6nccu2NGndngkTJmQ5HtabYfNmylWtxbinngm75tmf/I9Gbc6mQct2PPpkoN3mLVvpdn5/Grc9hybtO2WOef/991OrVq3M9cycOTPznO+++44OPfrSpF0nmnbozMGDMfUIygpxHg/HI9zfUiDPwz4RyZKVy3o8nEr4OPT5RrDHw4033kh6ejqdOnVi5syZrFu3jjPOOIP09HTq1KlDh9YtePeV53nx9bf4bP6XfDfvE9LSDpDc8VwG97+AxDq1AXjrrbcy5WTz5s3j9ttv54orrsgc9513NKSe2b+HW++6h4pBphrx8fGMe+ShQGj1Tt3p0aMHjRs3ZvPmzcyZMye0N8M/76ZPz3PDrjcjI4MbbrmNT6a9Q0LNGpzVrTcX9OlJ40YNdcyH7qdls6bs3buP1t370KNHD+135Ehuu+22LH2lp6czdOhQ3nhuAs1Sktm1ezelSxewWYfTrh53cJyc+rmGwzkU4B9BsMfDli1bmDp1Krfeemuma9IZZ5wBwCeffELVqlWpcqqmQCxfvhxpBw6Qnp7OgYMHKVO6NBXKlw85TqdOnahcOXTqRGMM7057j8FBuVBrVK9Oy+bN7FjlSUpKynQjGzlyJGPHjg3tzZCYSHJSI8Lh66XfUL9ePeol1qVMmTJc0v9CPpj5sR2zGi2bNc1cn3/MUPjkk09o2rQpzVLUe+LUypUpVZDclbOTOy7hiBy85vueaS4S5PFQIPA8HjyCdvPNN1OxYkVWrlzJggUL6Ny5M0uWLGHUqFEMGTKE1NRUVq9ZS6tufUnduJmTTzqRmk3aULdFe269/hoqn1IRjAEMV3Y5h6blT+KeM2ryx7kt2Tu0H0dT1/Fnj1ZZyoIvv6Ja1aqc2aA+IhKybNy0KWto9Vq1aNasWZa17N+/n8cee4z7brsZk5GOyTiCOXwwW9myaSMJNasjZcoiZcpSu04dftm+I/O3VzZu+y1zTFDtctOmTRk2bFhmIIE1a9YgIvTufwmtOp3L2Kee0fWHKjFxP3bb1YKAiNQQkWYi0k5EGojISbHsv8QTORu59DX788Igj4d4NKJwzBHs8fDbb79RtWpVypQpw969e+nYsSOPP/44gwYN4qGHHmLUqFGICK2aN2X+jCn8d+r77Ni5m63fLWb9kvk88fx/WJ+qmtc3n5/AV+0bM6v1mXy5Zx9vb9sddh6T352ajYvzY9++fQy47IqsodUffDBbu9GjRzNy5EjKlTs5x3VHEkrdP2aFChW47rrrWLduHStWrKBGjRrceuutgG5XFyxYwJsvPsv8WR/w/kez+PSLsAmg8g8RKBWftTjkGSJSWkQuFZF3ReQ3VBy0HM37sBr4Q0SWiMijvqQ1UcPdJcUDwFWoO9gI1F1EUBu5UqinREwR7PGwa9cufvjhB4wx/PLLLxw8eDDT4X3nzp0kJCTQvXt3Zn/2BY/fP4qTTjqJ06pUpnTp0lQ9rQodzmrF0m+/o15iHWrVqM5eoHx8KQZWr8w3f6bRoVJ2V6r0o4ZpH8zgm4VzQ87xyJEjDLjsCoZcMpD+/fvz/fffa2h1y8Vt2bKFli1b8vXXX7N48WKmTJnCHbdlsOePP0N6MwAk1KzJlq2BHKrBodSDxwSyhES/5pprMpU0CQkJdO7cmSqnatTjPj26s+zb7+jeOasSJWbwwp87RAUROQEYiaYV8KL77AQWoAq+g0Al9L1rhlo+3C4in6IxIKOy1yrxnFwQvkX9VJeiHhDvojfiNBEJ5w4WFYI9Hnr16kVaWhpPPvkk/fr1o1u3btx3333s37+fKlWq0KtXLxYvXkz90xNJT09n9+972Ln7d4wx7N+fxuJvltOovioodu5Szu3IUcPsnX+QVK5syDnM3f0njRo2IKFW9ly9xhiGX/d3khqeyS3/0NBPXmj11NRUUlNTSUhIYNmyZVSvXp358+eTmprKhu+WhvRm8NCmZXPWrlvPhtRUDh8+zNtTpmZ6WoQaE9SY2cN7772XqSXu1asX3333HWlpaaSnpzNv4Vc0bnhmdDckIogqHvzFISKISEeUS3sYDaJxA9DAGFPNGNPJGHORMWawMaa3MaYZavVwLpqjtTWwWESCo3lHBMfJZUUnYJwxxp+rdSLwFfCCiMzypTIsEAwbNoypU6eyaNEiLr30Us4880xSUlKIi4ujSpUqrFq9lhbd+nJej+788uuvpHTqiTGGqy4dSNPkJPbvT6P3JVdwaP3PZBjoUrk8C37fy2Prt7HrSDpJ87/nrno1uKJWFab+9juDr74h5DwWfrWISZPfISW5Mc3bdYS4UjEJrR4fH88zY8eEDKUebszJkyezYsUKRITExMTMMOmVKlXilltu4azuvRGEPj26c16vHvmaX47wFA8O0aArykT8xRjzbW6NjTEHUEP8z0XkZnSHdVc0A5d4jwcAEUlE0w/+AdQ2JmsAIxF5DQ22eZUxxosuXDQeDzsi8zDbe1nOzu0eKnwQYYThCMOahwudFAwpE5q7LKgx23TtydLl3+bP46FRfbP45bFZ6uI7DnAeDxFARE60hKvQ+3Db1axYFkzgLObazxaFOBeHYw0iSKnSWYpDZIiUOFmGI+TvaImkI3JZEY4l+NV+xi5LSxCGDRuWzSvh3XffJTk5mbi4uCwRb79etoIWXfvQomsfmnfpzXsfzc48NnnaBzTt3IsOi1bRf/nP7DocmVfa6jVrad6uY2apUL02T018jnsffIimZ3WgefPm9OzZk19++SX3zmKE2bNnZ8kZcUzAmZAUGESkNvCdiCTY3zXt78gcscPAEbmsCOcRXt1+RpiyKu+46qqrmD17dpa6Jk2aMG3aNDp16pS1vlFDlsyZwfLPZzHrnTcYcfso0tPTSU9P5+Z7HuSzaZP5sl1jksuV5cXNuTuvAzQ8swErFi1gxaIFfLPwC0468UQuuqAft9/8D777+ktWrFhBv379QpqPFAQyMjK44YYbmDVrFqtWrWLy5MmsWrWqUMYOC3GKh4KEMWYzmnjqflv1APCf4Cx6eYVTPGRFSxEpH2LL2sV+LqeA0KlTJ1JTU7PUJSVlD1Ukp9XFb4l2aN9RJK4UclpdjDEYhLQTK1P56YkcfGAsSY0bUm7QheEHDiH3+vSTTzijfn3qNgqkPzT7fmff7p2QfjjHVIGRRg7JDV9/9RX1T0+kXvVTIf0Al/a/kA/efYfGt9+SrW0kY36z4rtv8j0p56BfGBgNrBKRi4A+xCB7l7tjWVERuA/wa1dbA0NQLu69IppXNixevJhhw4axceNGJk2aRLx9+Z5//nlSUlI4+YQyNKhbm4n33pZLT9nx9ttvZ8lmP2rUKN54/TUqVqjAZx9Nj9kacsLWrVupnRAwbUmoVZPFS/NPp/INx70VCERktO/neuD/UKuGkZ77oDHmgWj6dttVxcX28yBwm4iki8gWEZmHhluKA/5W0OYjeUHbtm354YcfWLJkCY888ggHDx7kyJEjPP/88yxfvpwt82aQ0vAMHn0xb6ZFhw8fZvr06QwcODCz7uGHH2bTTyu5bNBAJr74UqyXEhKReEYUOkSgVOmsxSFWqOsrJ6Pv3MlB9VHBETnFQ/ZzK/ASsAWoiTro7wT6GmPeKaK55YikpCROPvlkVq5cyYoVKwB16hcRBvbuzpfLv89Tf7NmzaJly5ZZvAw8XDboYqZ9MCPEWbFHQkICm7cEnPODPSOKBoKUKpWlOMQGxphhxphhqD1cRdQrogJwne9YVHBETpFscz/UN8Zca6MAlwU+Q5UOsUttFQNs2LCB9HTVmm7cuJHVq1eTmJhIrVq1WLVqFTt27ADgf19+TdIZiXnqe/LkyVm2qllyTsycRaMQ+U8LAm3atGHtunUhPSOKDIJTPBQ8/gl8a4yZiBoP35HfDp1MDjDGrAtRd1hEngW6Ad2BqFxKIsXgwYOZO3dupp/qAw88QOXKlUPmNliwYAGPPvoopUuXJi4ujueee44qVTTXw+jRo+nUqROljx6hTs3qvDrm3ojnkJaWxpw5czI9CsCXcwJD3dq1eX7C+JivPRTi4+OZOP7xkJ4RRQfn8VCQEJFawK1AW1t1F7BERF41xmyJul/n8QDWDuefKDGrA5wY1ORuY0y2sEtF4fEQKY7+tCiidnGN2kXULieNqh9e4ueYIMJ8EZF4RohIvj0TWjdNNktmvJ2lLi6xqfN4iCFEpIMx5kv/b+Arkw9CVeI5ORGpB3yNRj+YD3yCalIzgETUneuEopqfw7EEcVvUAoafwIX6HQ2cTA5uQUOcXwNcAdxjjLnX5mb92LYpcXuUzZs307VrV5KSkkhOTmbCc5oa8d333qdJm/aUqnAqH83+mG59L6Bxq7Ykt25Huy7nkpKSQrVq1TjttNPC5mRo3749ycnJpKSkxCQnQ6F6Rkhc1uJwzKNE3yVRAxwv+Ni/0fyqf4jIVhs483J7rMRxvPHx8YwfP54ff/yRRYsW8dyLL7Pqp59okpTE1LfeoNPZHShVqhTjxvyLVd8sZtgVQ/lpzVreeecdrrjiCk444QSWLVvGihUrMiOXeDkZXnjhBX744Qfmzp2b75wMheoZIUCcZC0OeYaIjBaRv0TYtrOILBeR8JFdc0GJe3k9iEg1YBLgBSC7A3gTqAy0B0YB9e2xfEVPKApEKmsLhxo1alCjhppslC9fnqTkZH75fV9mYhlKxXNa7dNpZZPlbNj6K/UbNGDLmh84qRScWKY0S+Z+wlmtWmb2+fEn/6Np0pk0TayB2fMblUsBe3eGDUweiSdDpJ4RrVo0b5XXaxBiRm67GhvchyryPoig7RagKdAZmBrNYCWSkxORMsBHqKytLxoR+BHgKWA4cAlK4LxojUWt1itSpKamZsm3EAq1atXixx9/pFXzpuzZs4d1qRsZdOU1DLvhJn7fo4GV1/y8DsHmZOjcg7ET8p/xPpRnxNZtBRKxXuG2q4UKa/mwD2gebR8llZO7H+XYzjXG7BGRrqhBcF/0muwF3gcmoIH7+sUiHtbxiH379jFgwIDMfAvh2kyZMoXu3bvT/fwBVKt6Gr26d+WaK4eydMW33DpqNK88O4H0jAwWLFrM159/zElly3LuhQNp1bxZvsKVF6pnhIgjbLFDdRHpHGHbLWg49KhQ4oiciJyIWlXfYIzZA5kanG72+FXAYOBiY0yGiJyCekIMQSMklBgcOXKEAQMGMGTIkMx8C+HaDB06lFtuuQXzh0Y9ObvneTRsUJ8WTVM4/9KhACTUrEHnsztoTgZjYpKTobA9I8RtV2OFXrZEAoNKRKNCiSNyaBjmMoTY34tIW+A24BxjTAaAMeYPEXkXOJ8SROSMMQwfPpykpCRuuSV75I/gNiNGjGD//v2cBLw99X3iS8XTuFFDnnz2BZrYPKy9unfl8QnPkpaWRpnSpZm38Ctuvu7afM3T7xlRq2ZN3p4ylf++WlC3KXpOzsZIexDojWrzt6G7hQeMMZEZIRYvrEezcxU4SiKROwNYa4w57K8UkRrA60D/EA/dt8BfC2l+xwQWLlzIpEmTSElJoXlzFYeMGTOGQ4cOZXph9OzZkz179pCSkkLr1q1Zv349p516Kunp6VSsUJ5mHbqQWKc2Lzw1DoBKp5zCyBtGcFa33gjEJCdDoXpGeG5deT1N5Az0ha6KCtt/As5C/TN7i8jZxphdMZzp8YD5xpirC2Ogkkjk4glat02VNgX4pzEmlP1Bui0lBh07dgwp7wK46KKLwp7nbVfDYeglFzP0kottwufYoG/vnvTt3TNm/eWI6MxGnkMJ3D+MMc94lSLyBJqi72FUhOJQACiJUtT1QP2gLN3PAR8bY8KptJsAPxf4zByOcUietavWo6Ynmobv2aDDo9EUmJeLSM5ZuYsXUtHoPoWCkkjkPkddti4DEJEbUZeuf4VqbInhJWgOVoeSjOiikHSzn58YY476D9gI1AuBk4D8GTYeRzDGnOFP+1nQKHFEzga+fBV4yCaMHgZckYMD8Fjgd2BaIU3R4ZhF3jk5oKH9XBPmuBfLqiCzYpdolESZHMDdaPz4t4F2xph9wQ1EJA71ergG6GiMORLc5ptvvtkpIpElQnUoakQdWdbDN8tXfCzlKlUJqi4rIv5QNC8aY170/fZCpIQLqeLVn5Lf+TmERkklcodQA0MBPhKRF1CTklTUSLgDcCNQG40KvCRUJ8aY0wpltg7HBIwxvQugW0+TUSJjnlmZZVeU4/X/IawGPjPGbMjvGCWOyFmn/FdQP9VW6ulwAAAgAElEQVQ3UCf8q9DkNWXQh+1n1K/1BWPMjqKZqUMxgcephQt6VyGoXYmAJW7PoUoZyE7kxbabg4ZAXx/tWCWOyFn8zRiz335/GXjZEr+qwB5jzKGim5pDMcNq+xlO5ubFkw8nsyt2sBGAF6Lv2wY0huMasv4hnIkSwHOBr0SkpTFma4juch/PRQZ2cCg4WEPgn1FRyBl+DauIlEc9H+KA03x/vMUaIvIyunu6yeZyyKnt9cAzwOvRJrMpcdrVYx0i0khEkiNsGxNX9Fj145AdNorGJ2iU6RuCDj+Apt17o6QQOIvewPTcCByAMeY51P0tUj/XbHBE7hiCiLQDFgFXBxkrh8OJ9rwsBlsiUldEakYwXm0AY4yJIcHMddyCgPVaOVYJ9vXAduBpEXlfRB4Rkc9Qb4c1qBa/JOFU4Mc8tF9tz4kKjsgdIxCRRqjB8U/ATGNMWi7tWwGbRaS9jZYSZ+sboLZXY63sI6fxNorIfyBA6EQkWUTOi3INKbbPwknplRXDRSTes3c8loid5eZaA6+hmahuRX2onwbal0C/1W3kzfj5LAKxHfMOY4wrx0BBIxPvAy7z1dUH4sO0vxE4CuwB2ti6ONTeaj7q1fFvoFaY85sC39s+Jti6U4E0YDPQL4o1tAR22T7/FePrI+F+2zXPAXqg2dYqF/X9dCXHe/mEfUaeAMrl0K4cMN62fTLq8Yp6wa5k3tCn7c1sa38/iW5lGoZoG2c/b0Jt/vYHEbpKwEzb37+BhDBjNkO1XH5Cd73t83vgghzmK8GftrREBe05ErpgopWP6zbCXrsjqD/kLuBXVN7Vpajvqysh79kpwCr7jOwFZqPKhYfR4LXP2Lq9ts0q4JSoxyvqBbuSeePPAXajYdmftzf3FSDBI2q+tqV832+1L3gaytbH2XIq8J6P0NUJ1YclSgtsu4m2bpj9/QPwlxAELc7XT3U7x0TgZFt3lo/QPeJrGxfi/IpATaBmmOvSGg1ieisa59/rQ1Bbxt+AF9HAC9ehuXPvAb5DU03eUdT31pWQ97USGrBgH7rrCFX2orZ0lfI1VlEv1pXMm17OvsiHLXF4yxIOP0HogyYBmYmGcL/I1t9rCeR+oKWtOxXdkq5FOZxMji4E0TwfWGbHfdzWDbK/vwEu97WN932/HlgC7LDjvw9c65vrJpQAj/OdUzro/C/RLfcuNNxVTx8hu90SsaO+8jIq15qGbqtb27b3AUN8fSehAv0fgZuL+v66Eva5PxH1ePgbKrL5p/3eFTgxJmMU9SJdyXLD/21f5MPoNjLRd+we1FhyryUeB2zbb9Gt5Q4C7P8TqJZ2K+q+dtS2z7J1RXPOLraE0Dv/KCog7wh8iP6jGoJkInY+R+1cZgDLfedvscTR224YYEbQ+aPssY1oIMn1liBuA+709f8RcCGQgspnMlAO83egnq+/h4BRQWNUteN8CjQp6vvrStGUIp+AK/ZGQA3gC+AdVCZxEPgf0Bi42b7wk9DQ7FiOxiMgs22bqb6694HzgHqokmIxPhmdJSTpKMfYE03DeCW67TWoXO4LS3iMLY/Zsdta4jYBaGTrKlri6LWdgiYGetFX975t2xzl0J4Dkm1dHXS7uRol8vvQcPPJvmu0xDefH9Etvrft7g88ZL/7lRIploAOK+p77EoRvVtFPQFXMuVLcZb4VATKo/kADqLbuU2WG/EIQk/L0exD7a8SbX0rW2fIunWNB05HjVKPWgK0B9VI1vfN4+/2+Hb7+TjK0U3zEar7LeHYAbTynXs28Ked7yF7/iB7bKSP+E5Bt8drfevxZH1lgYGW0B4EhvrmvxTdEo+1xzJQDrCLbdMPNb2B7NvxJ1CCHUeMFB6uFPg70Rr9A061z9MBlIMfA1TMU19FvZiSWEK8hCfZT79CoYYldIctcbjN1gsBjul3oK6tPxENBHrYEkRv6+rJrEqhXOEKAtza7d58gFqosH6tj9AcRbeIpVGuyiN0y4A5vrkKGlj0KJrVzPv+DAH52m2+89cDi3xj+zmv2gS2uc/Z4/9DifK1qKxwNbrNPYoS4G6opviDoOvqEc+b7AsSl5f75EqRvR99COwkXkRNotKAcfb7OqB6xP0V9YJKWiG7IuF5S5j+A1xEVu1hdWCWJQxLUWH6RaiQ3gDn+9qWReVzR9Eox58TsKPrYNuVQU1TPK5qBVDDHquHcoH/Z3/fS0DG9irqfjTbR6h22/l4hOQW2/Zm29duS4wqENhSTvedvxOrhAgmPmioqzRgHipr24cqFsqjyb/XodtxT4P7Hhoia0rQ9Y23n4OAr/D9ibhy7Bb7XN4dVPcw8Lb9/irwWqT9OY+HQoSIiLEO2iJyDxq082qUUxqOytQe8Vnu/wb8157eDM0m9gSBkFCtRaSsbTsG3ZLuRYnnepTVrwC8a93EjqAEwQtb0wTlGLH9ZaDEAlR5sBO11bsS1f6+RIBIVQJusWODmmsAtEEVHr8DvxqNxOy1mWPrDar9fRjAGHNUROI8rw3b1247n7Ntf88YDRe+1F6vJuh29UZUuVAK5QIzk64aY7zkQx0JEGyHYx9norJpP/4P1bgCTETlvRHBEblChEcQRORmdCs6HbXSv5mAQP124H4foVuMGreCEpAawACsMB34q4iMQOVcn6LC+ZboAzEXVcf3MQE3sQtRAjINvf89re/rXmAl0E1E+qBc2iKUc6xvjPkN+BjlFufauQ4XkSG23xUo5zUYfQhPB9aKiEeQARqhGttrUe5suIhcY6/NUROI0DEQtZ1bixL3LSaQJnIvujU/xc7jeZQQ70ZNSrIQMhGpaufyujGx89F1KFBsQp9hP5oDv9jvO1GuPjIUNWtaHAs5CLfRf6l1KEFK8dqj8iePS7rQ1/5a9OX1OK0b0Zf8f7buiL35G1A3sDqo4sCgBK8eyuEkohrIr1F2/3IC5iX3owToGvv7VzvHD1EZm2eIe7MdqwrKgf5q249GjXCboRpbg08eaNfRByWiz9vf59k2+1BZSwLKRV6Ncmu/2fX8nx3jr6jccDhqZpIO9LB9eVv8yUC3oLq/oZxgy6J+LlyJ+P25ARWz/AsYinL8e4Eb7PGuwMqI+yvqBRXHQghFgu9Yb/vSXoYVuqPc0k/2hd+DakbPQremaah8bY49708Crlj7LUFJR7ehA4FH0W1qum1zEBXiHrYE8U/U5elLe75nB7cRNUv5wf72xvIb4m5HCWcFu5Yhluh5tn2LUa1rhq2bj7pcvWrbbQfO9F2Lvna9R1HOdAfqwvMNAc1rQ5TgemN4dnjT7XXwa3hfAwbY7wJcaucytKifCVfy/A6NQP8U99pn8nrfsYZYU6qI+irqxRS3gppx7EKjS4BPkWA/PTONjvb3RShHlGYJnWeUe9ASmxm+F/4OSww8ju83e45HlLzyoyV6HgE0BLS0ntLBM8Ld5iN0xhKd332/txEQ8BvUDKWrb72X2+OHgBdQ2VhbVKayz573B7r1bRzievWyY+5ETQZqEeSniHKhd9o+70W3nx1Rzm2vreuLcpwXoTLEZ+01vM5//V0peaXIJ1DcCmGig/iOn2+P/w3oZF9ug3I9N6JbQr8W8i57nqD2YpMt8RiOWvRfToDDecYSjT4EtqxbUPlfQ1Te5xGrV22/nhHuetvPRyiR/hwlsJtQ7m0nyglmoFrWLr41XW4Jyn7gBDvX0nbMnug2NmxkEDvfo6gJy8m5XF+/uUlN1DxkNVndvzaihtPnFvXz4ErRlyKfQHEshI4O4plRnEnALetXS3A+ABr4zr/VR+T2YbdkqLB9LbA4aLxBvpe8NwFj3X/Y8zehCoHzUfmeZ6x7qT2/LKoA2WjntQI4yx6rRMAQ9x/AUwQMcbv75jAQ31Y0imvWC+s9EUHb4LBLtVDu7hlUQ12OgMhAgtu7UrJKkU+gOBVyiA4S1M6zQfNkVz19xwaj2z9DwO5tL7oNPhXluFZajqk6kGTPu5SAp8JsdMt6Kirg34cqJuag28ZMY117rqBZ3B+y4y6w9cGGuGUskfYb4vYo6uvuu3ZtsVFPcIa/rthS5BMoTiX4xSJEdBDfsWcJcGtvWaI4xRIqj1h9ZYu3/T0X1XgeRd2YXkeF/RUsQRqMcm1HLYHta8fqR0A+9p0lVH5jXU9u2JCAEW5ZVOuaaYjrm3uwIW5MokXE4Pq3xm1RXQkqRT6B4lLIKnf7G2o4GxwdxPMl9QJMzvARuqOoqcMkS4Qm2boVPs5pKwEzj5329wtBBGg8ASXDPb76PqhCwTMZWUfAn9MjcvGoDO9zVKb2mSWElYLWutKOcwOWk3TFlWO1FPkEilux3NsBdGt4I+qc7kUHSUNlZh5RqY/KyLaicbROtxxUKdTI1dOwpqNaVC9Wm+cwn4HKyRLtOddbTm0PARndEHuOoHIvL1zTUdSVrAwBze9gW/80qtTYDXwatL7BluieX9TX2hVXIilFPoHiVAhEB3kLTRrsEbNTUEG9x7H1tPVjCZh4pKMa1psIxFrz7NQ8bu8Quh3900eojtp67/da1OVpCNmNdePRbbInC/wP4Y1whdCGuN+iHgYhI/m64sqxVop8AsWpEDDR6Onjjlr5CFEqga3rcz4C5o/Z5hGzb1BlwQ7gzSCithg1+t0TdN4ibLRgO3YoY91fUONKj4CGNMK154cyxN3ob+OKK8d6cb6rsUU1+7nbmEw/yb+iXFEZNMjkF2hEjxGoouFvKFH7E1U4YNu2RLerXYwxQ1GC420VG6KGwQtR49eWqHdEG2CoiHQFMMa8hW6DPQPi5aiLV0fU5ONPVPHwCWpC0t0Y84O3GGPMapQDvBtVMIy288ls4+BwzKOoqezxWAiTHo9AuKEHUNlaZVQ7uolAqr5ZBHxRt6N+qJ6ZRyq63fQ4thH47LxQM4+3CUTGzQw6SSBqbkTGurY+YiNcV1yJVfE933kq0Y4Xnz8SWfIgInEmEC4pEThsjPGiI0xHhf9DUFOP/ajt1mOoLOtpVPGwBSVSA4wxB0RkJmrn9jZKyDajHFcySsCM7b88auXvuX2tE5F4dJtaDzVD8cKdHxKRUsaYT40xk0TkIPCtMeaQtxZjzCwbcWSjMWZ/TC+Ug0N4bEaf60hRDdX2RwVH5PKAIAJ3LSqIP1VEOllCtxXdgt6BmpA8Z0890Rjztoi0QAlRLVTA74kLTkNlXwPR+HG10S1sT9R6f7+IeNrUhug203Nsn40qCu5AHdQT0S3s+UCGnfMcY8y7odZkjPk4/1fGwSFyGGPq5tZGRE5GxTPXoO+LyfmMnAd0JTIW2+8zOZpAdJCh/uOoG9QrBOzaDGpr9rz93E3AYNcLae4Z9Z5CwMwjA1gdYh4jCJh5PMBxZKzriiu5FTT6zksEFG1bUU+c06Pus6gXdbwVVJGQjhrhJvuImxAwGamE2sstJ6ux74/2nKEEtJ5Po1zcC6gcLw7d3npa1ycIb+bhjHVdOe4Lmrzp7wTcGA+jop8LiIF7XpEv8HgqBBQJ36ORc4PzMSSgMrGyqAyhM2of57lUTQCa23OGEDAB2W0JXzwq01tkj4U180ADVzpjXVeO22Lfl9cJxDVcj+bbjakNpseFOEQAETkLJUCPolvEo0bzE1yPclmJqOfBPPSmdUTlbRBQIOxANaGVgBYEQhOlEzDxWI+aliSgMdTS0dDjdwD7jTF7rHnKO8DFqPP8LFSp8Q9UY9vFBBQiDg7HHETkAlScshX4mzFmVoEMVNTU/HgpaEw0z1PgKV99uEzyBjXm9bK/TyFgFJyGbl1fQBURwUEnK/n6D2vmgTPWdeU4LmhCc0/+vBkNd1435uMU9UKPl4JqPj3H99XotjRUJvmrfQTnKDDO1vvt2D5HNadeyryQdmy+scPGWiNE1NyivlauuBJpQQ3jh6O5RzLQXctsNFlTfCzGcNvVINi0gcb79Ors4UGomUZZ+/k06i3Q2xjzjYhUQ29QI1Rb2g71XtiAbk/vJ2DH9iHwtDHmUzvGQNSObU0hLNPB4ZiDiDRFTUaGop4424E3gJeMMT9H3W9JJnLBBC3IDq46qgiIB34zxhyw9Y+jng1x6FZxhzGmvT12DWobtwclcG3QcOWHUcL2Gsp5zUFlb++j2avmFM6KHRyOfYhIWZShuAZNNA4wzxjTNfxZ4VHSjYErAnssgYs3NhlxKEWCiHxsjHneGHO7NVS8BjgDOF1E5hJIeLwLDW6ZiioBQG1+PgTSjTFrRGQ5unW9UIeTBR4RdXAo7hCRf0XYdC7KzfVDLRWiG6+kcnIi0gxYBlxh1JHdq78HTfyyBbXbSUCVDgATjTH/sO0eQzWgFVGlwQbU+f5xINUYkyEin6Ea1iXGmLN9Y6xEtaXrgc+MMT8W5FodHI4liEgGWd0VIzrNGBNVQJGSzMmdiQo6J4nIYWPMuyLSFjXHeIZAZvbywNnotvPvtu1txph/isj3KFGrhhrwvglkWM7wNNR6uzTwo80kfwT1US0PzDXGzCjMBTs4HCO4ujAHK7FEzhI1QQnaOyLyJ+qFcCIq7Fxj1AZuH8oyj0e1l7eIyHLL/b2F/huNR1256gILROQLVBt7CNUerUPt5Tw7Ni9enINDiYMx5o3CHK9EbleDFAyXolrSKugWMh7oZQmcALfZ4jfqXYVyfF9arm0IytFVJ6tR7zZ7XoKtj0ftgfoaF5PNwaFQUCI5OUvA4owxR41GBxF0u9kT5doaosa6d6MGirOAl9HcDdOAxqj3wQRUpvaWiMShjvu1USLXAVWDV0S3qC1Qn9I3jTEbCm2xDg7HMESkLrp97YiGETMoc7AQeMUYszHfY5RETs5DkEb1cjSqRyJKwLYAE1EzjwnGmO9tu4Oo7ZtBTUEeNcbM9fXxEirrq4zGmiu5F9jBIQeIyF2o7WhpsishBJVhP2CMGZOfcUoUJ+ffpgJ4BM5+n2Rjto215TvUs+FFY8z3IlIazZh1AuqkvxSN8pFrcEoHB4es8FkxHEJjKH6KMhaCine6oX7ZD1lm5MGoxyopjEaQHK41mk2rJkqsvjTGHLHHhqDhjaoCM+33eWjkkfqo32lblFtzRr0ODnmEiDRERTeefPqnMO2S0HewNuqPvTqqAYvad60wClkTP99OICepV14Gzva18afzO4z+w3jO72187VxwSldcyWNBZdlHgVYRtG2NMhTPRDtesd+uWpctj4MbRUCR8BJq2nEVcDNwioiUMcZ8brIrEgRNMNPLGJPq674+Gi/OM+p1XgsODrmjK7p7ytWMyhizVEQWAV2iHazYEzlj/w6sqchI1J7tSWNNOESkE/qvciFwsnVjnWtUvgZKDGuikUMyNT0iMhhn1OvgEA3qoH7ckWIJyoxEhWJP5ABsdJArCSgSfshjlqurUUPe4SLiD07pjHodHPKOsmi07EixDzXSjwolgsihe/qWwEPGmK/tVjTiLFci8h1qO/ci2Y16XfRdB4e8YReBROyRoBqwM9rBSgSRM8bstLGqvAz196GhkMYCk40xh4Gw0UGMMatFpDvOqNfBIRZYAXTKQ/suqIF9VCgxJiR+2OggtYB2xpjfffUuOoiDQwFDRK5CZeN9jTGzc2nbF00rcI0x5pWoxitpRE5EqqDRRZYbY7r76gejCWr+7hQJDg4FB2tYPx8N99/H7qTCtZuNKvjOMVEa2JdEIueyXDk4lCCUOCIHmRbXM4HTcdFBHByKNUokkQMQkUScIsHBodijxBI5BweHkoGoYqY7ODg4HC9wRM7BwaFYwxE5BweHYg1H5BwcHIo1HJFzcHAoVIjItSLyu4j0CqofKSJLbbk2ZuM57aqDg0NhQkSmoomeEowxGbZuKPA6Gnk7Hk0zcKExZnp+x3OcnIODQ2GjGbDAI3AW16PRuGujOR62ATfFYrBiTeREJFVEUotg3KtExFhH5MIa8347ZpfCGjPWEJEGIvKeiPxq17LH1h/3a8sLiuq5LURUB1K9HyJSAWgDvGOM2WOM2YU65SfFYrAcQy3Z7FV5wdXGmNein45DSYWIlEJzZtQHJqF5NQ4W6aQKCCIyF+hsjJGinksR4SCaUsBDB5Th+sJXtxNN65lv5BZP7oEQdTejCZMnAHuCjq2IxaSKAd5Dow5vK+qJHEc4HU3a/ZIxJljoPBF4G9hU6LMqGnTPvclxjY3AOb7fF9nPhb66qmhwzXwjRyJnjLk/uM5uwSoCTwUldXGwMMb8AfxR1PM4zlDTfmaLAGOM2Uk+IsMebzDGrCvqORQw/gs8LiIz0Mx5VwFf2PvsoSmaPCr/iCKdWCqa7ToxxLFBaI7SP4ADaK7Su4ATwvSTihLMicBWlI1dhYY9kjDj53mMEPUnAHeiCaTTgD/R+FaDwowpqBB0lZ3jVjvniqHGQG+aAa4K0VcC8DR6Aw8Cu4GvgXuD2nVFw62vsvM7gAYSGA2UDdHv/XbMLnm4lycB/0Rzz+5FY+n/aOdXLahtDeBZu97DwA5gGiHSyvnXb9cx1/b/J/ARkBTU3oQp9+e2NjR95DJ7fbajW92adkwT6X3xzWNuuOsKXAYsttcpNajfqWiw1QN2nQuBoUF9Jeaw1rm+dtmeqbw+t76xXrPf30b/KA7a+90vr+9+rApwsr2OGbbs8D9H6DZ1MXBdLMaLWfhzERmDEpudKKXeB/QBxgC9RKSHsQmcfSgD/A84xd6EMsAAdCvcEM1Qn98xgudZBo3+2xn4CX1xT0Ljy70jIs2NMXcHnfYscB3KZbyIvuQXAGcBpdGQ6bnCJrX+GL2J81AicRK6TbsfTZfo4Z9AI+BLlDCUBc627bqIyLkmq3YqTxCRSmgGsmZogp9X7LrOAIbZuf1m254OLECJx2fAZFQLNhA4T0QGGGM+DDFMP+AvaMy+F+w6+wJtRKSxCfxzP4C+iFeicpm5tn4uOUBEbkdD2P+Omh/8AfRACUysOelbbd8z0OtW0XfsefTPaB4qojgVXeckEWlojLnXttuDrvUqoC5ZxUGpOQ0e5XOLHedrlABPQp+9S4AP7DP0eQRrjymMMftFpANwLrqG+b5nAWPMbjSBe8wGzCsVTiWIkwPa27pNQHVffTz6UBjg7jD9LMDHhaE3YZ091ikGYwRzWXfZtjOBeF99Vd+cOvjqz7F1q4FTfPVl0IfahBjjKoI4Btt+g62/LMR1rR30ux4huFmUEBrgkqD6+8kDJ4f+SRj0BY0LOlYeqOj7/bFtOyqoXQc0Ht8uoFyI9acD3YPOecQeuyOovgs+7i23tdnrcwTlAmr76gUlwobYcnL7gRZhzjsjRF0Z4FM7x1pBx+YGz60AnttE7xoAo4P66uX1FcmzcryXvJ8Qmsi9ZOuuDdH+TJQlXR+mn3NCnOM9jK/GYIzgh2Utmme1UYh+htsxXvHV/cfWXRGi/dlETuQG2LoP8nXDlEvIMkdb772MXSLoo6q9Xr8AJ+fSNsH2uxEoHeL4pODr41v/myHan26PTQmq70LeiNw9tu6+EO3rogTW5HZfgo7nROSejOJe9Q/17BAdkcvrc5voPZtAqRDnbAR25udZLMiCJroZHYu+YmUn19J+fhZ8wBizBjUHOF1ETgk6nI5ux4Ix1362iMEYmRCR8qiJwi/GmJ9CNPH69o/rfV8Qov0idA2RoJ39nBVJYxE5WUTuFpElIvKHiBy1Jj0eW18rwnFDoQ2qsp9njNmfS1tv/fNNaFFAqGvmYWmIus32s1Kus4xsXtnuizFmo2+cWOHrcAdEpI6IPCsiP4lImrXpM6icDvJ3r6J9bj2sMKHFGpvJ/z3IF+wzniAidYMLmjFvtIgk+uqiQqxkcp58IpzJxDY0a3ZFspqd7AxzA34N6jc/Y+R1nqAywuBzfgtubIzJEJFI1dxen1tza2gTeHyGyvxWojkpdhCQ/Y1GhdDRIuK5EN0185DtPhhj0jXNBqUiGDuSeWW7L776xHyO4cevoSpFpB5KACuhSoBPUHlgBgE5Y37uFcT4HlikU0TOAPaavYnK3XKyFTSoLNFDVPONFZHzhLzVUXlaMGoEtfNQxWarDyZ01UO0j3aMcPMMhVB9/Gk/q5H1gnsGrKcSGbHwHrZI/tX/ghK4140xVwWNWQMlcvlBXuYSzTUrDPjvS6i8HKGSFx+1n9me+5x2ABYmTP0t6DOQzRDeZoC7Mpd+I8Gxeg+ixZPo8z0PlbGHYnSa2/JafgeLFSX3Er92CT4gIvVRuc4GY0zwv0o8KrwOhtePP6FstGNkwhizFyWQtUSkQYgmXe3nshDjdgzRvh2R/1Essp99Imhb335ODXGsc4Tj5YSv0Re+k4icnEvbzPWLSKi1hrpmhYGw98VubWqHOMfLsRvqWOso5xHNvfKc0iPiZqN8bo9ldADeNsZ0NcZcaYwZFlxQ7xeC6qJCrIicl/T1HhE5zau0N3GcHeflMOc+IiIn+M6pjAqVAV6N0RjBcxXUGDHzIbP5WO/1tfHwhv0cJSIVfe3LoKYrkWIGKgS+wP7DZ4GI+LmqVPvZJahNPeCxPIwZEsaYHajJTg1gnIhkeQ5EpJy3VmPMFmAOuvW6OahdW9R27HfUy6Mw8V90y3WjiGQSLZty8hFCb4eXosT9MhE5yXdOZdQUJRqk2s8u/kobRuivYc7xRBx18jBOXp/bYxmV0NzHhYKYbFeNMV+KyFjgDmCliEzBJo4FmqDC4cdDnLoNlVesFJHpqM3ZxejL95wxZl4MxgjGOHvOX4BvRWQmaqszENU6jjXGZAqzjTFfiMiLaI7WH2yYmCPA+ej24BcC26CcrtFhERmIymz+KyJ/Q7m7sqgjcncC92MG8DNwi4ikoFxLHdTu7CPy9nKEw9/R6zYCtbv7GLWTOx01MbiAgAJoBGp79riI9ESJhWcndxTdqu2NwZwihjFmnYjch/7RfCsi7xCwk6sMfItazfvP2SYibwGXAytE5COgAmrTNo/Qgvvc8BxwNfCufTa2ote1N/B/qE1aMD5Fr9P636YAAAW6SURBVN00+/wdADYaYyblME6enttjHELoLWrBIArVbirhPR4uRYnNXtSy+gdgFKEt9FMJeDw8iz4ch1CL+5w8HvI8Roj6ssDdqFD/gO1rATA4zJhxwEjUCPMQStietXPfi2qw/O2vIrzHQx30xdiAEpVdqHV3sA1abeAte10O2HXegRLCnEwduuThXp5sr51nQb8XNWp9Cqga1LYWalO30c57J7qlaBOi37Drt8dDzb8LeTAh8R27HP0TOIgqZ95EjZZXAntCtD8B/TPcYtfxM2qDFvV1Rbdfn6EcrfcsXRhuTSiXOQaV8R4JHjcWzy0+j4cwc55LDmYsBVmAK4CUXNo0JYTZVjSlyIJmeqFkjDGJRTKBGMDKR9ag8oVsW1CHooEN3fMb+ufTvqjn41C0KNbx5GIFEakeQm51EsrxQOHLoxwAETnNmtv46+KB8SjX4+7LMQwRaSgiz1hb0DUi8pWIPBlGuRL9OI6Tyx0i8igwGGXxt6Gq/O6oRncWcJ4pqgtZgiEiI4AHUf/nzagsrhPqAbMCdXM6UHQzdAgHERmCKkpC6QXSgRHGmFdDHMszYuagX8wxB3Vk74m+SOnoNvVpNOSUI3BFg8WoTKoTaqsGKut8GHjMEbhjEyKSjFpCbEBl3fNQu8dxKCPxFPCiiKwyxizO93ju/XRwcChMiMgrwFAg2Riz1tZlAA8aYx4Qkeqor+7HxpiL8zuek8k5ODgUNroCsz0CFwxjzK/AdDQARr7hiJyDg0Nhozq5GwN7MtZ8wxE5BweHwsZ+1PA/JzQgfPCFPMEROQcHh8LGL+TgtSOaevICNEBovuEUDw4ODoUKEXkedUerYWyMQqt4+ArN03EBysW1NsbkO+Od4+QcHBwKG++gxMyfx8GgKQ76AR8AbWNB4MBxcg4ODscAROQcNDHVGpN7tOq89e2InIODQ3GG83hwcHAoUOQnP4PRfB35G99xcg4ODgUJq1TIKZdDWBhj8q03cJycg4NDQeM1shO501Gf4z/QsO2/oXk5WqJxGucTlFMlWjhOzsHBoVAhIg3RVKQvAw/4FQ0258i9aCTu9saY1fkezxE5BweHwoSIvAdUNsaETcokIl8AvxtjLszveM5OzsHBobDRCc0ZkhMWAOfEYjBH5BwcHAobJ6B5OHJCDfKflBtwRM7BwaHwsQwYJCIhc92KSEs0y9nyUMfzCieTc3BwKFSIyLlo2gCD5s/9HPgVDcHUBRiCMmB9jDFz8j2eI3IODg6FDRG5GPg3cApZzUsMsAfN8fBuTMZyRM7BwaEoICLlgIvQpN4VUZu55cB7xph9MRvHETkHB4fiDKd4cHBwKNZwRM7BwaFQISLXisjvItIrqH6kiCy15dqYjee2qw4ODoUJEZkKdAASjDEZtm4o8DqQhvrUnwBcaIyZnt/xHCfn4OBQ2GgGLPAInMX1qBlJbSAB2AbcFIvBHJFzcHAobFQHUr0fIlIBaAO8Y4zZY4zZBcwAkmIxmCNyDg4OhY2DZLWN64DSoi98dTtxeVcdHByOU2wkq/P9RfbT77RfFdgVi8Fc0EwHB4fCxn+Bx0VkBhos8yrgC2PMTl+bpsDaWAzmiJyDg0Nh4wVgENDX/t4N3O4dFJHK6Hb2nVgM5kxIHBwcCh0iUgo4FzgJmB/ExcV2LEfkHBwcijOc4sHBwaHAISKNRKRsPvtoGc15jsg5ODgUBgYB60TkRmsXFzFEpI2ITAOWRjOw2646ODgUOESkDvA4cDFwCPgA+BhYBKw2PkLkMw4+G7gUaIR6QNxjjHk1z2M7Iufg4FBYEJEOwK3A+QSsO7xAmQdQA+CyvvpNaHDNp40xaVGN6Yicg4NDYUNEagEXAN2A5qjx74mop8MmYB4wB/ifySeRckTOwcGhWMMpHhwcHIo1HJFzcHAo1nBEzsHBoVjDETkHB4diDUfkHBwcijUckXNwcCjW+H8YUcUE9Gvo0AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "protein161.set_unfolding_path('cts1_2uy3_len_161')\n",
    "protein161.runHMM(s=.01,m=0.1, starting_state = 'abcdef')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use a kinetic threshold of 100 (note this is in units of millions of MC steps, so 100 means 100 million)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A value of T_A = 150 will be used to form clusters\n",
      "Doing loop clustering...\n",
      "We have obtained the following clusters: \n",
      " {0: ['abcdef'], 1: ['abcde', 'abcd', 'abd', 'ab', 'a'], 2: ['bcd', 'bd', 'd', 'b'], 3: ['∅']} \n",
      " Please make sure that unfolding rates fit well to Arrhenius equation before proceeding \n"
     ]
    }
   ],
   "source": [
    "protein161.form_clusters(T_A = 150)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pathway \t \t fraction of flux\n",
      "[0, 1, 2, 3] \t \t 0.704 \n",
      "\n",
      "[0, 2, 3] \t \t 0.172 \n",
      "\n",
      "[0, 1, 3] \t \t 0.121 \n",
      "\n"
     ]
    }
   ],
   "source": [
    "pathways, unique_pathways, fracs =  dbfold.folding_rates.get_unfolding_pathways(protein161.temp_unfolding_info['combined_trajs'], protein161.temp_unfolding_info['PDB files'], goal_state = -1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see how the resulting unfolding rates fit to the Arrhenius equation. Importantly,the Arrhenius temperatures are should be multipled by the Boltzmann constant (i.e. in units of energy) prior to inputting to this function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computing temperature 0.9\n",
      "Computing temperature 0.925\n",
      "Computing temperature 0.95\n",
      "Computing temperature 0.975\n",
      "Computing temperature 1\n",
      "0.9618645819520123\n",
      "0.996660195726096\n",
      "0.995671820761867\n"
     ]
    }
   ],
   "source": [
    "%matplotlib qt\n",
    "\n",
    "activation_energies, prefactors, mean_transition_rates, Ns=dbfold.folding_rates.Arrhenius_fit([ 0.9, 0.925, 0.95, 0.975, 1], [(0,1), (1,2), (2,3)],\n",
    "                                                                         protein161.temp_unfolding_info['combined_trajs'],\n",
    "                                                                         protein161.temp_unfolding_info['PDB files'], \n",
    "                                                                         'Pathway 1',\n",
    "                                                                         legend_loc='lower left', min_trans = 3, fontsize = 20,\n",
    "                                                                         labelsize = 20, legend_fontsize = 20, show_rsquared = False,\n",
    "                                                                         colors = ['firebrick', 'mediumblue', 'gold'], temp_norm = 1, \n",
    "                                                                         labels = ['f unfolding', 'a unfolding', 'b and d unfold'])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[119., 120., 120., 120., 120.],\n",
       "       [108., 109., 104., 102.,  96.],\n",
       "       [ 47.,  71.,  98.,  99., 109.],\n",
       "       [ 18.,  50.,  92., 109., 120.]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Ns\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we agree to this Arrhenius fit, then we will go ahead and compute folding rates and obtain an error distribution by bootstrapping. Note that if you set N_trials=1000 as below, this will take a few minutes to run.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Folding info already exists! Loading now...\n",
      "Folding info successfully loaded\n"
     ]
    }
   ],
   "source": [
    "protein161.obtain_folding_rates(Arrhenius_temps=[ 0.9, 0.925, 0.95, 0.975, 1], unfolding_transitions_of_interest= [(0,1), (1,2), (2,3)], N_trials=1000, min_trans=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us now plot some of these predicted folding rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib qt\n",
    "\n",
    "fontsize = 20\n",
    "labelsize = 20\n",
    "TM = 1\n",
    "\n",
    "folding_rates.plot_folding_rates('cts1_2uy3_len_161/Folding_info.dat', ['b and d fold', 'a fold', 'f fold'], [(3,2), (2,1), (1,0)], fontsize = fontsize, labelsize = labelsize,              \n",
    "                                                               legend_fontsize = fontsize, temp_norm = TM, legend = True)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1e-17, 1e-12)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "folding_rates.plot_folding_rates('cts1_2uy3_len_161/Folding_info.dat', '', [(2,1)], fontsize = fontsize, labelsize = labelsize,              \n",
    "                                                               legend_fontsize = fontsize, temp_norm = TM, legend = True, ax = ax, labels = ['length 161, [bcd, bd, d, b] -> [abcde, abcd, abd, ab]'])\n",
    "\n",
    "folding_rates.plot_folding_rates('cts1_2uy3_len_82/Folding_info.dat', '', [(3,2)], fontsize = fontsize, labelsize = labelsize,              \n",
    "                                                               legend_fontsize = fontsize, colors = ['green'], temp_norm = TM, legend = True, ax = ax, labels = ['length 82, [null] -> [a]'])\n",
    "\n",
    "folding_rates.plot_folding_rates('cts1_2uy3_len_50/Folding_info.dat', '', [(1,0)], fontsize = fontsize, labelsize = labelsize,              \n",
    "                                                               legend_fontsize = fontsize, colors = ['gold'], temp_norm = TM, legend = True, ax = ax, labels = ['length 50, [nul] -> [a]'])\n",
    "left, right = plt.xlim()\n",
    "plt.xlim(left = .5, right = .7)\n",
    "bottom, top = plt.ylim()\n",
    "plt.ylim(bottom = 10**(-17), top = 10**(-12))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Kinetic model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In general, the kinetic model class is designed for modelling co-translational folding, but it can also be used to solve master equation for a single protein (at its fully synthesized length). We will do that for protein 161. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "km = kinetic_model.Kinetic_model('1igd', [61], folding_rate_paths = ['1igd_data/Folding_info.dat'], omit_transitions = [(3,0)])\n",
    "#km.rates[0][3,0,:] = 0  #delete connections between 3 and 0 (abc and abcd) \n",
    "#since our method does not accurately predict this folding rate\n",
    "#km.rates[0][0,3,:] = 0\n",
    "km.diagonalize()\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now plot the solution to the master equation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "kinetic_model.plot_timecourse(km, 0.8, [5*10**(10)], [[4,3,2,1,0]], starting_state = -1, final_time = None, colors= ['k','gold', 'mediumblue', 'purple', 'firebrick'], \n",
    "                 linewidth = 8,colorchange = (None, None, None), scientific = True, ylim = (), ax = ax,ntimes = 200,fontsize=20, \n",
    "                 labelsize = 20, timenorm = 1, ylabel = True, xlabel= True)\n",
    "\n",
    "ax.legend(labels = ['$\\emptyset / b$', '(b)c', 'a(b)c', '(b)cd', 'a(b)cd'], fontsize = fontsize)\n",
    "ax.xaxis.get_offset_text().set_fontsize(20)\n",
    "\n",
    "ax.set_xlabel('Time (MC Steps)', fontsize = fontsize)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import importlib\n",
    "importlib.reload(dbfold.utils)\n",
    "importlib.reload(dbfold.folding_rates)\n",
    "importlib.reload(dbfold.compute_PMF)\n",
    "importlib.reload(dbfold.folding_rates)\n",
    "importlib.reload(dbfold.analyze_structures)\n",
    "importlib.reload(dbfold.nonnative_states)\n",
    "importlib.reload(dbfold.kinetic_model)\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}