simulate a nanodisc

nanodisc.py

import numpy as np
import matplotlib.pyplot as plt
import mrcfile
import random

#################################################################################################

nanodisc_diameter = 170	# Angstrom, change this value for nanodisc of different sizes

#################################################################################################
# no need to change parameters below

msp_thickness	= 10	# Angstrom
core_diameter = 40 	# Angstrom

max_shift = nanodisc_diameter/2.0 - msp_thickness - core_diameter/2.	# Angstrom

# in theory, the nanodisc can be centered around the core -> min_shift = 0
# however, it seems that we need set min_shift > 0 to make the average to appear smaller
# it means that the nanodisc might be always off-centered
min_shift = max(0, max_shift - 20)	# Angstrom


def buildNanoDisc(n, nanodisc_diameter, msp_thickness):
	# http://mathworld.wolfram.com/Torus.html
	m = np.zeros((n, n, n), dtype=float)

	z, y, x = np.mgrid[-n/2:n/2, -n/2:n/2, -n/2:n/2] * apix
	x2 = x*x
	y2 = y*y
	z2 = z*z

	r2 = np.square((nanodisc_diameter-msp_thickness)/2.-np.sqrt(x2+y2))+z2

	m = np.exp( -r2/(msp_thickness*msp_thickness/4.))
	return m

n = 160	# 3D map size in pixels
apix = 2.0	# Angstrom per pixel

nanodisc = buildNanoDisc(n, nanodisc_diameter, msp_thickness)
with mrcfile.new('nanodisc.%gA.mrc' % (nanodisc_diameter), overwrite=True) as mrc:
	mrc.voxel_size = apix
	mrc.set_data(nanodisc.astype(np.float32))

avg = nanodisc*0
K=5000
i=0
while i<K:
	dx = random.uniform(-max_shift, max_shift)
	dy = random.uniform(-max_shift, max_shift)
	r2 = dx*dx+dy*dy
	if r2>max_shift*max_shift: continue
	elif r2<min_shift*min_shift: continue
		
	i+=1
	print("%d/%d: min,max shift=%g-%g Angstrom\t\tshift=%g (dx=%g\tdy=%g) Angstrom" % (i, K, min_shift, max_shift, np.sqrt(r2), dx, dy))

	dxi = int(dx/apix + 0.5)
	dyi = int(dy/apix + 0.5)
	avg += np.roll(nanodisc, shift=(dyi, dxi), axis=(1,2))
	
avg /= K
with mrcfile.new('nanodisc.%gA.avg.mrc' % (nanodisc_diameter), overwrite=True) as mrc:
	mrc.voxel_size = apix
	mrc.set_data(avg.astype(np.float32))

plt.figure(0)
plt.imshow(nanodisc[n/2,:,:])	# torus itself
plt.gray()

plt.figure(1)
plt.imshow(avg[n/2,:,:])	# average of many nanodisc instances randomly shifted in X-Y plane

plt.gray()
plt.show()