Introduction to Python Programming.pptx

Introduction to Python Programming.pptx

• 1.
  +1 (315)557-6473 https://www.pythonhomeworkhelp.com/ support@pythonhomeworkhelp.com Kolb Rd., Tucson, Arizona, USA
• 2.
  Exercise 1: BiologicalComputation Part 1. Accompanying this assignment is some python code that implements Dan Gille- sipie’s “Direct Method” for computing the exact behavior of discrete chemical systems (gillespie.py). Modify the code so that it implements Dan’s “First Reaction Method Part 2. You should now have two discrete reaction event simulators. Use each simulator to compute the behavior of the system specified by the provided input files (system1.inp, species1.inp) over a period of 100 seconds, sampling the system state every 0.1 seconds. Report how long it takes each simulator to run and whether or not this time agrees with your expectation of the computational costs of each algorithm. Part 3. Produce a plot of the numbers of proteins A->D over time (100 seconds). Also, draw a simple sketch of the species and the reactions that define this system. What do you expect this system to do? Is the computed behavior what you expect? Python Homework Help
• 3.
  Solution: The important differencesbetween Gillespie's first reaction and direct methods: In the direct method, we consider the time and identity of the next reaction independently. First, we picked one random number (uniformly on (0,1)) to generate the time of the next reaction (with exponential distribution parametrized by a = sum(a_i) ): # how long until the next reaction r1 = random.random() tau = (1.0/(a+eps))*math.log(1.0/r1) Then we picked a second random number uniformly on (0,a) and identified which reaction's "bin" it fell into (where the width of each reactions's bin is equal to its current propensity). r2a = random.random()*a a_sum = 0 for i in a_i: if r2a < (a_i[i]+a_sum): mu = i break a_sum += a_i[i] Where a_sum is the "right edge" of the current bin you're considering. Python Homework Help
• 4.
  In the firstreaction method, one picks a random number (uniform on 0,1) and uses it to generate a a time (exponentially distributed with parameter a_i) for each reaction, and selects the reaction with the smallest time to fire at that time. Tau's for all reactions are recalculated at every step, never saved. Note that the behavior of the first and direct methods is equivalent. mintau = t_max eps = math.exp(-200) # which reaction will happen first? # caluculate each reaction's time based on a different random number for rxn in a_i.keys(): ai = a_i[rxn] ri = random.random() taui = (1.0/(ai+eps)) * math.log(1.0/ri) if taui < mintau: # "sort as you go" mu = rxn # the putative first rxn tau = taui # the putative first time mintau = tau # reset the min This method of sorting seemed fastest to us, but here's another way of tackling the problem: tau_i = {} # which reaction will happen first? # caluculate each reaction's time based on a different random number for rxn in a_i.keys(): Python Homework Help
• 5.
  ai = a_i[rxn] ri= random.random() taui = (1.0/(ai+eps)) * math.log(1.0/ri) tau_i[taui] = rxn tau = min(tau_i.keys()) mu = tau_i[tau] Coding # python implementation of "direct" gillespie method. from gillespie_helper import * import math import random import sys # # examples of the dictionaries: # consider RXN1: A + 2B -> C with propensity 0.01 # # substrate_dict['RXN1'] = {'A':1, 'B':2} # product_dict['RXN1'] = {'C':1} # c_dict['RXN'] = 0.01 # Python Homework Help
• 6.
  # counts_dict storesspecies' abundances. substrate_dict, product_dict, c_dict = ReadInputFile(sys.argv[1]) counts_dict = ReadCountsFile(sys.argv[2]) # lists storing time-steps and species concentrations at each step t_list = [0] s_list = [counts_dict.copy()] # how long a time period to simulate t_max = 1000 # calculate both a_i and a def GetA_i(substrate_dict,c_dict,counts_dict): a_i = {} a = 0 # here, we first calculate h and then take product with c for rxn in substrate_dict: a_i_rxn = 1 substrate_rxn = substrate_dict[rxn] Python Homework Help
• 7.
  # calculate h,the reaction propensity by successively taking # product of all species' counts involved in this reaction. # note that in most cases, substrate_rxn[species] is 1 and # NChooseK(n,1) just equals n. for species in substrate_rxn: a_i_rxn *= NChooseK(counts_dict[species],substrate_rxn[species]) a_i_rxn *= c_dict[rxn] a += a_i_rxn # store reaction probabilities in a dictionary a_i[rxn] = a_i_rxn return a_i, a # calculate the initial a_i a_i, a = GetA_i(substrate_dict,c_dict,counts_dict) # the meat of the algorithm; here, we step through time, separately # finding the time until the next reaction and the next occurring # reaction. eps = math.exp(-200) while t_list[-1] < t_max: Python Homework Help
• 8.
  # if thereare no species left, the simulation should end if a == 0: t.append(sys.maxint) break mintau = t_max # which reaction will happen first? # caluculate random time for each based on a for rxn in a_i.keys(): ai = a_i[rxn] ri = random.random() taui = (1.0/(ai+eps)) * math.log(1.0/ri) if taui < mintau: mu = rxn # the putative first rxn tau = taui # the putative first time mintau = tau # reset the min # update the substrate and product species counts involved for species in substrate_dict[mu]: counts_dict[species] -= substrate_dict[mu][species] for species in product_dict[mu]: counts_dict[species] += product_dict[mu][species] # record the events Python Homework Help
• 9.
  t_list.append(tau + t_list[-1]) s_list.append(counts_dict.copy()) #update a_i a_i, a = GetA_i(substrate_dict,c_dict,counts_dict) # end of while loop # write the results to an output file names = s_list[0].keys() output_file = open(sys.argv[3],"w") output_file.write("time ") for k in s_list[0].keys(): output_file.write("t" + k) output_file.write("n") # specify how many data points to use data_points = int(sys.argv[4]) interval = len(t_list)/data_points if interval < 1: interval = 1 Python Homework Help
• 10.
  for i inrange(0,len(t_list),interval): output_file.write(str(t_list[i])) for j in names: output_file.write("t" + str(s_list[i][j])) output_file.write("n") output_file.close() Python Homework Help