monte carlo simulation of protein structure and grid

Question

I am working Monte Carlo simulation script over protein structure. I have never done before Monte Carlo scripting. I will extent this program at large scale. According to protein xyz coordinates I have to define the box size. This box will be divided into a grid of size 0.5 A. Based on distance and angle criteria I have to assign the point based on Boltzmann probability distribution.

My program should be move in each direction by taking grid of 0.5 A and generate the random point and check the condition of distance and angle. If satisfy the condition put point there otherwise discard that point based on Boltzmann probability distribution.

Here is my code for generation of random points

from __future__ import division    
import math as mean    
from numpy import *   
import numpy as np   
from string import *    
from random import *    

def euDist(cd1, cd2):# calculate distance
    d2 = ((cd1[0]-cd2[0])**2 + (cd1[1]-cd2[1])**2 + (cd1[2]-cd2[2])**2)
    d1 = d2 ** 0.5
    return round(d1, 2)

def euvector(c2,c1):# generate vector
    x_vec = (c2[0] - c1[0])
    y_vec = (c2[1] - c1[1])
    z_vec = (c2[2] - c1[2])
    return (x_vec, y_vec, z_vec)


 for arang in range(1000):  # generate random point
        arang = arang + 1
        x,y,z = uacoord
        #print x,y,z

        x1,y1,z1 = (uniform(x-3.5,x+3.5), uniform(y-3.5,y+3.5), uniform(z-3.5,z+5))
        pacord = [x1,y1,z1]                 # random point coordinates
        print pacord

I am completely struck to generate the box size from the xyz coordinates of protein structure and how to define the grid of size 0.5. How to check every point in the box.
Any help will be appreciable.

Heteropolar

So, you have a bunch of points in 3D and you want to generate a box which contains them all? – Wholesome

As an aside, I sincerely hope you know what you're doing: protein folding typically has mutliple minima, and naive MC simulations tend to get stuck. – Wholesome

defining the box and grid size i will generate point and then will do further calculation. How to define the box and grid do not know??????? how my loop will move from each direction????? – Heteropolar

Out of curiosity why are you not using one of the current, free, and GPU enabled MC programs? The ones that I have used will randomize starting points in any fashion you wish. – Petersburg

Have you looked at github.com/ndexter/Aeolotopic-Monte-Carlo-Simulation/blob/… ? Although it is in C, you will get some idea about implementation and it would help you in writing python code, I believe. – Glochidium

What did you end up doing? – Clown

Answer 1

Did you consider using PyRosetta? It is much easier to use since a lot of the functions you need are already built in. You can visualize your output in real time in PyMol. I wrote a similar script in PyRosetta, fairly easy to write and modify and it did what it was supposed to do.

Reckon

I am aware of PyRosetta. I am working on a big project and developing a program. So I can't use such available program. If possible can you provide your script to me or give me some idea how to generate box and define the grid. Thank you. – Heteropolar

Answer 2

I love your topic, question, and approach. I don't know how long it'll stick around here.

A 3D mesh in a rectangular space is common in finite element analysis and all other techniques for analyzing physics problems. Have a look at how meshes are generated.

There are usually two parts: a mesh of (x,y,z) points in space and boxes that connect them. A simple 2D mesh with two elements would look like this:

D               E               F
o (1, 0) ------ o (1, 1) ------ o (1, 2)
+               +               +
+               +               +
+               +               +
o (0, 0) -------+ (0, 1) -------+ (0, 2)
A               B               C

There are six points in this mesh:

A (0, 0)
B (0, 1)
C (0, 2)
D (1, 0)
E (1, 1)
F (1, 2)

and two boxes:

1 - (A, B, E, D)
2 - (B, C, F, E)

So one possible approach would be to iterate over all the boxes, check position at the centroid, and adjust accordingly.

I'd externalize the mesh definition from your code and read it in from a file. That way you can handle different meshes with the same code.

If I understand you correctly, the mesh will remain fixed in space. You're trying to randomize the motion of the protein itself. It's like a fluids problem, an Eulerian approach: you're tracking the motion of material against a fixed control volume. Except it's a protein instead of a fluid.

So you'll also have a separate definition of the initial condition for the protein in space. You plan to increment it in time to see how the protein folds according to your rules.

Am I close?

Olivero

Thank duffymo. I am not looking the motion of protein but try to assign a bunch of points which in in close proximity of the protein and form a interaction with protein atoms. I did it with random point generation but with that program every time i got new result. So i thought this approach. – Heteropolar

Answer 3

My code was written for a protein folding application, but the overall idea is the same. It starts with a certain temperature and decreases the temperature step wise, the protein (or in your case the 'point) is moved randomly, if the energy score of the new position/structure is lower than the previous one it is accepted, if not the pose will be evaluated according to the Metropolis distribution. You have to define your scorefunction(), a function to define a random start position and a mover which moves your point from its original position. The code below is modified from my original script, just to give you a rough idea.

kT_lower=0.1
kT_upper=100
ktemp=kT_upper

max_iterations=15000

i=-1

#create random start point
pose=create_pose()

#evaluate start point
starting_score=scorefunction(pose)

while i<max_iterations:
    i=i+1
    new_pose=random_move(pose)
    if scorefunction(new_pose)<scorefunction(pose):
        pose=new_pose
    else:
        #linear decrease of kT
        ktemp=kT_upper-i*(kT_upper-kT_lower)/max_iterations
        #exponentatial decrease of kT
        #ktemp=math.exp(float(i)/float(max_iterations)*float(-5))*kT_upper+kT_lower

        try:
            p=math.exp(DeltaE/ktemp)
        except OverflowError:
            p=-1

        if random.random()<p:
            pose=new_pose
            print str(i)+'; accept new pose, metropolis'
        else:
            print str(i)+'; reject new pose!'

Reckon

Looks and sounds very reminscent of optimization using simulated annealing. – Olivero

It is simulated annealing, or at least how I understood the concept of simulated annealing. – Reckon

Hi Ashafix. I am stuck again. Can you provide me your email address so that i can take your help without spam this site. I have same query. I got the idea but found difficulty in implementation. – Heteropolar

Hi Awanit, I'd be happy to help. Just check my profile and send me a message via other profile page. – Reckon