Search Results

Search found 2 results on 1 pages for 'user359597'.

Page 1/1 | 1 

  • Python/Biophysics- Trying to code a simple stochastic simulation!

    - by user359597
    Hey guys- I'm trying to figure out what to make of the following code- this is not the clear, intuitive python I've been learning. Was it written in C or something then wrapped in a python fxn? The code I wrote (not shown) is using the same math, but I couldn't figure out how to write a conditional loop. If anyone could explain/decipher/clean this up, I'd be really appreciative. I mean- is this 'good' python- or does it look funky? I'm brand new to this- but it's like the order of the fxns is messed up? I understand Gillespie's- I've successfully coded several simpler simulations. So in a nutshell- good code-(pythonic)? order? c? improvements? am i being an idiot? The code shown is the 'answer,' to the following question from a biophysics text (petri-net not shown and honestly not necessary to understand problem): "In a programming language of your choice, implement Gillespie’s First Reaction Algorithm to study the temporal behaviour of the reaction A---B in which the transition from A to B can only take place if another compound, C, is present, and where C dynamically interconverts with D, as modelled in the Petri-net below. Assume that there are 100 molecules of A, 1 of C, and no B or D present at the start of the reaction. Set kAB to 0.1 s-1 and both kCD and kDC to 1.0 s-1. Simulate the behaviour of the system over 100 s." def sim(): # Set the rate constants for all transitions kAB = 0.1 kCD = 1.0 kDC = 1.0 # Set up the initial state A = 100 B = 0 C = 1 D = 0 # Set the start and end times t = 0.0 tEnd = 100.0 print "Time\t", "Transition\t", "A\t", "B\t", "C\t", "D" # Compute the first interval transition, interval = transitionData(A, B, C, D, kAB, kCD, kDC) # Loop until the end time is exceded or no transition can fire any more while t <= tEnd and transition >= 0: print t, '\t', transition, '\t', A, '\t', B, '\t', C, '\t', D t += interval if transition == 0: A -= 1 B += 1 if transition == 1: C -= 1 D += 1 if transition == 2: C += 1 D -= 1 transition, interval = transitionData(A, B, C, D, kAB, kCD, kDC) def transitionData(A, B, C, D, kAB, kCD, kDC): """ Returns nTransition, the number of the firing transition (0: A->B, 1: C->D, 2: D->C), and interval, the interval between the time of the previous transition and that of the current one. """ RAB = kAB * A * C RCD = kCD * C RDC = kDC * D dt = [-1.0, -1.0, -1.0] if RAB > 0.0: dt[0] = -math.log(1.0 - random.random())/RAB if RCD > 0.0: dt[1] = -math.log(1.0 - random.random())/RCD if RDC > 0.0: dt[2] = -math.log(1.0 - random.random())/RDC interval = 1e36 transition = -1 for n in range(len(dt)): if dt[n] > 0.0 and dt[n] < interval: interval = dt[n] transition = n return transition, interval if __name__ == '__main__': sim()

    Read the article

  • Python/Biomolecular Physics- Trying to code a simple stochastic simulation of a system exhibiting co

    - by user359597
    *edited 6/17/10 I'm trying to understand how to improve my code (make it more pythonic). Also, I'm interested in writing more intuitive 'conditionals' that would describe scenarios that are commonplace in biochemistry. The conditional criteria in the below program is explained in Answer #2, but I am not satisfied with it- it is correct, but isn't obvious and isn't easy to implement for more complicated conditional scenarios. Ideas welcome. Comments/criticisms welcome. First posting experience @ stackoverflow- please comment on etiquette if needed. The code generates a list of values that are the solution to the following exercise: "In a programming language of your choice, implement Gillespie’s First Reaction Algorithm to study the temporal behaviour of the reaction A---B in which the transition from A to B can only take place if another compound, C, is present, and where C dynamically interconverts with D, as modelled in the Petri-net below. Assume that there are 100 molecules of A, 1 of C, and no B or D present at the start of the reaction. Set kAB to 0.1 s-1 and both kCD and kDC to 1.0 s-1. Simulate the behaviour of the system over 100 s." def sim(): # Set the rate constants for all transitions kAB = 0.1 kCD = 1.0 kDC = 1.0 # Set up the initial state A = 100 B = 0 C = 1 D = 0 # Set the start and end times t = 0.0 tEnd = 100.0 print "Time\t", "Transition\t", "A\t", "B\t", "C\t", "D" # Compute the first interval transition, interval = transitionData(A, B, C, D, kAB, kCD, kDC) # Loop until the end time is exceded or no transition can fire any more while t <= tEnd and transition >= 0: print t, '\t', transition, '\t', A, '\t', B, '\t', C, '\t', D t += interval if transition == 0: A -= 1 B += 1 if transition == 1: C -= 1 D += 1 if transition == 2: C += 1 D -= 1 transition, interval = transitionData(A, B, C, D, kAB, kCD, kDC) def transitionData(A, B, C, D, kAB, kCD, kDC): """ Returns nTransition, the number of the firing transition (0: A->B, 1: C->D, 2: D->C), and interval, the interval between the time of the previous transition and that of the current one. """ RAB = kAB * A * C RCD = kCD * C RDC = kDC * D dt = [-1.0, -1.0, -1.0] if RAB > 0.0: dt[0] = -math.log(1.0 - random.random())/RAB if RCD > 0.0: dt[1] = -math.log(1.0 - random.random())/RCD if RDC > 0.0: dt[2] = -math.log(1.0 - random.random())/RDC interval = 1e36 transition = -1 for n in range(len(dt)): if dt[n] > 0.0 and dt[n] < interval: interval = dt[n] transition = n return transition, interval if __name__ == '__main__': sim()

    Read the article

1