#!/usr/bin/python import sys import random import math import matplotlib matplotlib.use('TKAgg') import numpy as np from matplotlib import pyplot as plt from matplotlib import animation # We start with the construction of the objects. # Every object has a mass_center. # Each of these has x, y, and z coordinates and d (diameter). # Arguably, we could attach the linkers to points on the periphery of # the particles and incorporate rotation of these bodies. # Defining parameters. membrane_elast = 15. Temp_scaling_factor = 100 epsilon = .3 # Lenard Jones Epsilon for Particle attractions. # Bond Lengths and Constants in nanometer and piconewton per nanometer: cas2tmem = 131. _cas2tmem = 30. tmem2tf = 65. _tmem2tf = 20. cpf2pmem = 98. _cpf2pmem = 30. pmem2prot = 45. _pmem2prot = 30. # Construction of the Cas9 component. cas9 = { 'x':-100., 'y':131, 'z':0., 'old_x':-100., 'old_y':100., 'old_z':0., 'r':90., # nanometer 'm':163., # Kilodalton 'bonds': { 'tc_membrane': { 'len':cas2tmem, 'con':_cas2tmem } } } # Construction of the trans-membrane domain. TC_membrane = { 'x':-100., 'y':0., 'z':0., 'old_x':-100., 'old_y':0., 'old_z':0., 'r':20., # nanometer 'm':70., # kilodalton 'bonds': { 'cas9': { 'len':cas2tmem, 'con':_cas2tmem }, 'TF': { 'len':tmem2tf, 'con':_tmem2tf } } } # Construction of the transcription factor. TF = { 'x':-100., 'y':-55., 'z':0., 'old_x':-100., 'old_y':-60., 'old_z':0., 'r':12., 'm':46., 'bonds': { 'tc_membrane': { 'len':tmem2tf, 'con':_tmem2tf } } } # Construction of the Cpf1 component. cpf1 = { 'x':100., 'y':95., 'z':0., 'old_x':100., 'old_y':130., 'old_z':0., 'r':60., 'm':126., # kilodalton 'bonds': { 'pc_membrane': { 'len':cpf2pmem, 'con':_cpf2pmem } } } # Of the trans-membrane domain. PC_membrane = { 'x':100., 'y':0., 'z':0., 'old_x':100., 'old_y':0., 'old_z':0., 'r':20., 'm':70., 'bonds': { 'cpf1': { 'len':cpf2pmem, 'con':_cpf2pmem }, 'protease': { 'len':pmem2prot, 'con':_pmem2prot } } } # And the protease component. protease = { 'x':100., 'y':-51., 'z':0., 'old_x':100., 'old_y':-50., 'old_z':0., 'r':10., 'm':27., 'bonds': { 'pc_membrane': { 'len':pmem2prot, 'con':_pmem2prot } } } # All stuffed together in one object. ProteaseChain = { 'cpf1': cpf1, 'pc_membrane': PC_membrane, 'protease': protease } # Combination of the above objects into one group. TargetChain = { 'cas9': cas9, 'tc_membrane': TC_membrane, 'TF': TF } # Which we of-course store in an object dictionary. objects = { 'TargetChain': TargetChain, 'ProteaseChain': ProteaseChain } # Create a dictionary that contains the forces that each particle undergoes, but start without # any actual forces on anything: directions = ['x','y','z'] forces = {} for chain in objects: forces[chain] = {} for particle in objects[chain]: forces[chain][particle] = {} for direction in directions: forces[chain][particle][direction] = 0. # Now, request seed values. random.seed(input("Please provide a random number generator seed: \n")) timestep_max = input("Please provide the maximal number of steps you want to make: \n") timestep_size = input("Please tell me how long one timestep should be: \n") output_doc = open('bouncy.out', 'w') output_doc.close() # Now, we actually start with the true fuctions of updating the coordinates. # At the end of each time-step, we increase the number of timesteps we have # done so far. # def VelocityVerlet(objects, timestep_size, forces): # This is basically just an ordinary velocity verlet thing, using only the input of forces. for chain in objects.keys(): for particle in objects[chain].keys(): # Abbreviate the particle you look at because it's easier to write. point = objects[chain][particle] # Calculate the velocities for the verlet integration, because I can: x_velocity = (point['x']-point['old_x']) y_velocity = (point['y']-point['old_y']) z_velocity = (point['z']-point['old_z']) # Abbreviate the full notation because that's easier to write. force = forces[chain][particle] # Calculate the accelerations for the verlet integration, because, you know: x_acceleration = force['x']/point['m'] y_acceleration = force['y']/point['m'] z_acceleration = force['z']/point['m'] # Perform the actual velocity verlet calculation, resulting in new coordinates. x_new = point['x']+x_velocity*timestep_size+(x_acceleration*timestep_size**2)*0.5 y_new = point['y']+y_velocity*timestep_size+(y_acceleration*timestep_size**2)*0.5 z_new = point['z']+z_velocity*timestep_size+(z_acceleration*timestep_size**2)*0.5 # Update the coordinates of each particle. point['old_x'] = point['x'] point['x'] = x_new point['old_y'] = point['y'] point['y'] = y_new point['old_z'] = point['z'] point['z'] = z_new def ClearForces(forces): # This clears the full forces dictionary, used at the beginning of each time-step to overwrite # the whole thing in one go. for chain in forces.keys(): for particle in forces[chain].keys(): for direction in forces[chain][particle].keys(): forces[chain][particle][direction] = 0.0 def MembraneRepulsion(forces, objects): # This calculates the proximity to the membrane for non-membrane particles and determines # the repulsive force they undergo from it. for chain in forces.keys(): for particle in forces[chain].keys(): if (particle != 'tc_membrane') and (particle != 'pc_membrane'): forces[chain][particle]['y'] += 0.01/(objects[chain][particle]['y']**2) # Coordinate is negative -> force is negative. Always pushes away from asymptote. def MembraneAttraction(forces, objects): # This calculates the proximity to the membrane middle-region for membrane particles and # determines the elastic force forcing them back among the 'lipid layer'. for chain in forces.keys(): for particle in forces[chain].keys(): if (particle == 'tc_membrane') or (particle == 'pc_membrane'): forces[chain][particle]['y'] += - membrane_elast * objects[chain][particle]['y'] def BondForce(forces, objects): # This calculates the distances between two particles with pythagorean math, then finds the force # and reverse engineers it into sub-vectors using another set of pythagorean equations. for chain in objects.keys(): for particle1 in objects[chain].keys(): for particle2 in objects[chain][particle1]['bonds'].keys(): # For each defined bond, calculate the effect of the bonded particle on the other. # Calculate the distances per normal. dist_x = objects[chain][particle2]['x']-objects[chain][particle1]['x'] dist_y = objects[chain][particle2]['y']-objects[chain][particle1]['y'] dist_z = objects[chain][particle2]['z']-objects[chain][particle1]['z'] # Use those distances to find in-planar distances. dist_xyz = math.sqrt(dist_z**2 + dist_y**2 + dist_x**2) # Use the distances and compare them to the required length to get a force. bond = objects[chain][particle1]['bonds'][particle2] if dist_xyz > bond['len']: force_xyz = -bond['con']*(dist_xyz-bond['len']) else: force_xyz = 0. # Reverse vectorize the force according to the distance vectorization. force_z = force_xyz*(dist_z/dist_xyz) force_x = force_xyz*(dist_x/dist_xyz) force_y = force_xyz*(dist_y/dist_xyz) # Implement these forces into the forces dictionary. forces[chain][particle2]['x'] += force_x forces[chain][particle2]['y'] += force_y forces[chain][particle2]['z'] += force_z def ParticleRepulsion(forces, objects): for chain1 in objects.keys(): for particle1 in objects[chain1].keys(): for chain2 in objects.keys(): for particle2 in objects[chain2].keys(): if particle1 != particle2: dist_x = objects[chain2][particle2]['x']-objects[chain1][particle1]['x'] dist_y = objects[chain2][particle2]['y']-objects[chain1][particle1]['y'] dist_z = objects[chain2][particle2]['z']-objects[chain1][particle1]['z'] # Use those distances to find in-planar distances. sigma = objects[chain1][particle1]['r']+objects[chain2][particle2]['r'] dist_xyz = math.sqrt(dist_z**2 + dist_y**2 + dist_x**2) # WolframAlpha derivative of LJ potential. (Force from potential is its derivative) force_xyz = -6.*4.*epsilon*sigma**6*((dist_xyz**6-2.*sigma**6)/(dist_xyz**13)) # Yes, this is pure mathemagic. # Reverse vectorize the force according to the distance vectorization. force_z = force_xyz*(dist_z/dist_xyz) force_x = force_xyz*(dist_x/dist_xyz) force_y = force_xyz*(dist_y/dist_xyz) # Implement these forces into the forces dictionary. forces[chain1][particle1]['x'] += force_x forces[chain1][particle1]['y'] += force_y forces[chain1][particle1]['z'] += force_z def TemperatureRandom(forces): # Add a random term to your forces that is scaled by the temperature. The temperature scaling # factor calculation should be performed seperately, or can be defined as a parameter. for chain in forces.keys(): for particle in forces[chain].keys(): for direction in forces[chain][particle].keys(): forces[chain][particle][direction] += Temp_scaling_factor * random.uniform(-1.,1.) forces[chain][particle]['z'] = 0. # # This is the actual simulation part. For each timestep, it recalculates all the forces and then # applies them using the functions defined above. Each function adds something to the forces dict # and these are then applied to the particles to recalculate their positions with a VelocityVerlet # scheme. # def step(timestep): if timestep != timestep_max: ClearForces(forces) MembraneRepulsion(forces, objects) MembraneAttraction(forces, objects) BondForce(forces, objects) ParticleRepulsion(forces, objects) TemperatureRandom(forces) VelocityVerlet(objects, timestep_size, forces) output_doc = open('bouncy.out', 'a') for chain in objects.keys(): for particle in objects[chain].keys(): output_doc.write('%i %s %f %f %f %f %f\n' %( timestep, particle, objects[chain][particle]['x'], objects[chain][particle]['y'], objects[chain][particle]['z'], objects[chain][particle]['m'], objects[chain][particle]['r']) ) output_doc.close() # create a figure with axes for each particle fig = plt.figure(figsize=(8,8)) ax1 = fig.add_subplot(111, aspect='equal') ax2 = fig.add_subplot(111, aspect='equal') ax3 = fig.add_subplot(111, aspect='equal') ax4 = fig.add_subplot(111, aspect='equal') ax5 = fig.add_subplot(111, aspect='equal') ax6 = fig.add_subplot(111, aspect='equal') ax7 = fig.add_subplot(111, aspect='equal') plt.axhline(y=0, color="red", alpha=0.5) ax1.set_xlim([-150, 150]) ax1.set_ylim([-100, 200]) ax7.transData.transform([(0,1),(1,0)]) # assign axes to particle names cpf1, = ax1.plot([], [], 'bo', ms=60) pc_membrane, = ax2.plot([], [], 'bo', ms=20) protease, = ax3.plot([], [], 'bo', ms=10) TF, = ax4.plot([], [], 'bo', ms=12) cas9, = ax5.plot([], [], 'bo', ms=90) tc_membrane, = ax6.plot([], [], 'bo', ms=20) # define the function which first calls step function to update coordinates # and subsequently adds new coordinates to respective axes # return axis after each step def animate(i): timestep = i+1 step(timestep) cpf1.set_data(objects['ProteaseChain']['cpf1']['x'], objects['ProteaseChain']['cpf1']['y']) pc_membrane.set_data(objects['ProteaseChain']['pc_membrane']['x'], objects['ProteaseChain']['pc_membrane']['y']) protease.set_data(objects['ProteaseChain']['protease']['x'], objects['ProteaseChain']['protease']['y']) TF.set_data(objects['TargetChain']['TF']['x'], objects['TargetChain']['TF']['y']) cas9.set_data(objects['TargetChain']['cas9']['x'], objects['TargetChain']['cas9']['y']) tc_membrane.set_data(objects['TargetChain']['tc_membrane']['x'], objects['TargetChain']['tc_membrane']['y']) return cpf1, pc_membrane, protease, TF, cas9, tc_membrane # call the animation ani = animation.FuncAnimation(fig, animate, frames=timestep_max, interval=10, blit=False) plt.tight_layout() plt.show() ''' # can also make animation from a data file live update graphs from data file --> https://www.youtube.com/watch?v=7w8jk0r4lxA def animate(i): pullData = open('bouncy.out','r').read() dataArray = pullData.split('\n') xar = [] yar = [] for eachLine in dataArray: timestep = eachLine.split(' ')[1] print timestep '''