In [ ]:
# IMPORT PACKAGES
from PIL import Image, ImageDraw
import numpy as np
from os import path, makedirs
import time
from tifffile import imsave
%matplotlib inline
import matplotlib.path as pltpath
from matplotlib.pyplot import imshow
In [ ]:
# SUB-ROUTINES
def motor_boundary_check(x,y,params,boundary):
if ~boundary.contains_point(np.array([x,y])):
h = params['boundary_height'];
w = params['boundary_width'];
b = [np.less(y,-1*h),np.greater(y,h),np.less(x,-1*w),np.greater(x,w)];
# b[0] = y position below bottom boundary
# b[1] = y position above top boundary
# b[2] = x position left of left boundary
# b[3] = x position right of right boundary
if b[0] == True and b[2] == False and b[3] == False:
y = -1*h;
elif b[1] == True and b[2] == False and b[3] == False:
y = h;
elif b[2] == True and b[0] == False and b[1] == False:
x = -1*w;
elif b[3] == True and b[0] == False and b[1] == False:
x = w;
elif b[0] == True and b[2] == True:
x = -1*w;
y = -1*h;
elif b[0] == True and b[3] == True:
x = w;
y = -1*h;
elif b[1] == True and b[2] == True:
x = -1*w;
y = h;
elif b[1] == True and b[3] == True:
x = w;
y = h;
return x,y
def motorX_state_forward(j,J,params,boundary,check):
# check if motor ungergoes activation event
if J[0,j] == 0:
# define state transition
if params['activationRate']*params['tstep'] > np.random.uniform():
J[0,j] = 1
# define motor diffusion movement
J[3,j] = J[3,j] + np.sqrt(2*params['tstep']*params['diffusionRate'])*np.random.randn()
if check == 1:
J[3,j] = J[3,j] + (0.006 * np.random.randn() + 0.003)
J[4,j] = J[4,j] + np.sqrt(2*params['tstep']*params['diffusionRate'])*np.random.randn()
J[3,j],J[4,j] = motor_boundary_check(J[3,j],J[4,j],params,boundary)
# check if motor undergoes bundling event
elif J[0,j] == 1:
# condition 1: look for motor candidates to bundle with
candidate1 = np.where(J[0,:] == 1)[0]
# condition 2: look for motors within motorBundleRadius
candidate2 = np.where(np.sqrt((J[3,candidate1]-J[3,j])**2 + \
(J[4,candidate1]-J[4,j])**2)<params['motorBundleRadius'])[0]
# group the candidate lists together (returns indices of candidates)
bundle_candidates = [candidate1[i] for i in candidate2]
# check the probability of a bundling event
if (len(bundle_candidates)*params['tstep']*params['bundleRate']>np.random.uniform()
and len(bundle_candidates)>1):
# uniformly randomly select a motor to bundle with
J[1,j] = bundle_candidates[int(np.ceil(len(bundle_candidates)*np.random.uniform()))-1]
while J[1,j]==j:
J[1,j] = bundle_candidates[int(np.ceil(len(bundle_candidates)*np.random.uniform()))-1]
J[0,j] = 2
else:
# diffusion of motor occurs if no bundling event
J[3,j] = J[3,j] + np.sqrt(2*params['tstep']*params['diffusionRate'])*np.random.randn()
if check == 1:
J[3,j] = J[3,j] + (0.006 * np.random.randn() + 0.003)
J[4,j] = J[4,j] + np.sqrt(2*params['tstep']*params['diffusionRate'])*np.random.randn()
J[3,j],J[4,j] = motor_boundary_check(J[3,j],J[4,j],params,boundary)
return J[0,j],J[1,j],J[3,j],J[4,j]
def find_filaments_2D(Z,x,y,params):
# test cadidates
d0 = x - Z[0,:];
d1 = Z[3,:] - Z[0,:];
d2 = y - Z[1,:];
d3 = Z[4,:] - Z[1,:];
d = (d0*d1 + d2*d3) / \
(d1**2 + d3**2)
T1 = (d>0) & (d<1)
var = ((d0*d3)-(d2*d1))**2 / \
(d1**2+d3**2)
dpp = d0**2 + d2**2
dmm = (x-Z[3,:])**2 + (y-Z[4,:])**2
T4 = np.where((var < params['r']**2) & (T1 == True) | \
((dpp < params['r']**2) ^ (dmm < params['r']**2)) & (T1 == False))[0]
return len(T4),T4
def update_motor_position(x1,y1,x2,y2,xm,ym):
# find closest position on line to motor as a fraction of the filament length
var = ((xm-x1)*(x2-x1) + (ym-y1)*(y2-y1)) / ((x2-x1)**2 + (y2-y1)**2);
if var>1:
var = 1;
elif var<0:
var = 0;
# calculate new x and y positions for motor
new_xm = x1 + var*(x2-x1);
new_ym = y1 + var*(y2-y1);
return new_xm,new_ym
def motorX_attach_detach(lfil,rfil,lx,ly,rx,ry,Z,params,boundary):
# toggle switches; disallows right motor from performing event if switched (left motor event occured)
fil_switch = 0;
diffuse_event = 0;
detach_event = 0;
# check left motor events
if lfil<0: # left motor not bound to filament
# find number of filament binding cadidates and their indices for left motor
l_num_fils,l_fil_indices = find_filaments_2D(Z,lx,ly,params);
if l_num_fils*params['tstep']*params['p1']>np.random.uniform():
# select a random filament for left motor to bind to
lfil = l_fil_indices[int(np.ceil(l_num_fils*np.random.uniform()))-1];
# set new position of newly attached left motor
lx,ly = update_motor_position(Z[0,lfil],Z[1,lfil],Z[3,lfil],Z[4,lfil],lx,ly);
# switch fil_switch toggle
fil_switch = 1;
else:
diffuse_event = 1;
lx = rx;
ly = ry;
else:
# check if both motors are attached to the same filament
# to see if left motor switches binding to a different filament
if lfil == rfil and np.random.uniform()>0.5:
fil_switch = 1;
# find number of filament binding cadidates and their indices for left motor
l_num_fils,l_fil_indices = find_filaments_2D(Z,lx,ly,params);
# probability of binding event
if l_num_fils*params['tstep']*params['p1']>np.random.uniform():
# select a random filament for left motor to bind to
lfil = l_fil_indices[int(np.ceil(l_num_fils*np.random.uniform()))-1];
# set new position of newly attached left motor
lx,ly = update_motor_position(Z[0,lfil],Z[1,lfil],Z[3,lfil],Z[4,lfil],lx,ly);
# check for detachment/unbinding event of left motor
if np.random.uniform()>0.5:
detach_event = 1;
# probability of left motor detachment from filament
if params['p0']*params['tstep']>np.random.uniform():
lfil = -1;
lx = rx;
ly = ry;
# check right motor events
if rfil<0: # right motor not bound to filament
rx = lx;
ry = ly;
# find number of filament binding cadidates and their indices for right motor
r_num_fils,r_fil_indices = find_filaments_2D(Z,rx,ry,params);
# probability of binding event
if r_num_fils*params['tstep']*params['p1']>np.random.uniform():
# select a random filament for right motor to bind to
rfil = r_fil_indices[int(np.ceil(r_num_fils*np.random.uniform()))-1];
# set new position of newly attached right motor
rx,ry = update_motor_position(Z[0,rfil],Z[1,rfil],Z[3,rfil],Z[4,rfil],rx,ry);
if lfil<0:
lx = rx;
ly = ry;
# check if both motors in bundle unbound to a filament; diffusion event occurs
elif diffuse_event == 1:
lx = lx + np.sqrt(2*params['tstep']*params['diffusionRate'])*np.random.randn();
ly = ly + np.sqrt(2*params['tstep']*params['diffusionRate'])*np.random.randn();
# check if new positions are within boundary
lx,ly = motor_boundary_check(lx,ly,params,boundary);
rx = lx;
ry = ly;
else: # right motor bound to filament
# check if right motor and left motor are bound to same filament
# and if left motor did not have a filament switch event
if rfil==lfil and fil_switch == 0:
# find number of filament binding cadidates and their indices for right motor
r_num_fils,r_fil_indices = find_filaments_2D(Z,rx,ry,params);
# probability of binding event
if r_num_fils*params['tstep']*params['p1']>np.random.uniform():
# select a random filament for right motor to bind to
rfil = r_fil_indices[int(np.ceil(r_num_fils*np.random.uniform()))-1];
# set new position of newly attached right motor
rx,ry = update_motor_position(Z[0,rfil],Z[1,rfil],Z[3,rfil],Z[4,rfil],rx,ry);
# check if both detachment event for left motor and filament switch of
# either motor did not occur
if detach_event == 0 and fil_switch == 0:
# probability of right motor detachment from filament
if params['p0']*params['tstep']>np.random.uniform():
rfil = -1;
rx = lx;
ry = ly;
return lfil,rfil,lx,ly,rx,ry
def bound_motorX_movement(lfil,rfil,lx,ly,rx,ry,Z,params):
# left motor not bound, right motor bound to filament; right motor walks
if lfil==-1 and rfil!=-1:
cb = np.sqrt((Z[3,rfil] - rx)**2 + (Z[4,rfil] - ry)**2);
# motor travels past end of filament, falls off
if cb>params['fil_length']:
rx = Z[0,rfil];
ry = Z[1,rfil];
rfil = -1;
# motor continues to travel along filament
else:
s = np.sqrt((Z[0,rfil] - rx)**2 + (Z[1,rfil] - ry)**2) - (params['v']*params['tstep']);
rx = Z[0,rfil] + (s*np.cos(Z[2,rfil]));
ry = Z[1,rfil] + (s*np.sin(Z[2,rfil]));
# update location of bundled motor
lx = rx;
ly = ry;
# right motor not bound, left motor bound to filament; left motor walks
elif lfil!=-1 and rfil==-1:
cb = np.sqrt((Z[3,lfil] - lx)**2 + (Z[4,lfil] - ly)**2);
# motor travels past end of filament, falls off
if cb>params['fil_length']:
lx = Z[0,lfil];
ly = Z[1,lfil];
lfil = -1;
# motor continues to travel along filament
else:
s = np.sqrt((Z[0,lfil] - lx)**2 + (Z[1,lfil] - ly)**2) - (params['v']*params['tstep']);
lx = Z[0,lfil] + (s*np.cos(Z[2,lfil]));
ly = Z[1,lfil] + (s*np.sin(Z[2,lfil]));
# update location of bundled motor
rx = lx;
ry = ly;
# both motors bound to separate filaments
elif lfil!=-1 and rfil!=-1 and lfil!=rfil:
# calculate travel distance for left motor
cb = np.sqrt((Z[3,lfil] - lx)**2 + (Z[4,lfil] - ly)**2);
# motor travels past end of filament, falls off
if cb>params['fil_length']:
lx = Z[0,lfil];
ly = Z[1,lfil];
lfil = -1;
# motor continues to travel along filament
else:
s = np.sqrt((Z[0,lfil] - lx)**2 + (Z[1,lfil] - ly)**2) - (params['v']*params['tstep']);
lx = Z[0,lfil] + (s*np.cos(Z[2,lfil]));
ly = Z[1,lfil] + (s*np.sin(Z[2,lfil]));
# calculate travel distance for right motor
cb = np.sqrt((Z[3,rfil] - rx)**2 + (Z[4,rfil] - ry)**2);
# motor travels past end of filament, falls off
if cb>params['fil_length']:
rx = Z[0,rfil];
ry = Z[1,rfil];
rfil = -1;
# motor continues to travel along filament
else:
s = np.sqrt((Z[0,rfil] - rx)**2 + (Z[1,rfil] - ry)**2) - (params['v']*params['tstep']);
rx = Z[0,rfil] + (s*np.cos(Z[2,rfil]));
ry = Z[1,rfil] + (s*np.sin(Z[2,rfil]));
# calculate new distance between bundled motors
motor_stretch = np.sqrt((lx-rx)**2 + (ly-ry)**2);
# check if bundled motors have passed stretch threshold
# if so, both motors "fall off" and unbind
if motor_stretch>params['r']:
lfil = -1;
rfil = -1;
# equal chance of one motor moving to position of other motor
if np.random.uniform()>0.5:
lx = rx;
ly = ry;
else:
rx = lx;
ry = ly;
# both motors bound to the same filament
elif lfil==rfil and lfil!=-1:
# calculate travel distance for left motor
cb = np.sqrt((Z[3,lfil] - lx)**2 + (Z[4,lfil] - ly)**2);
# motor travels past end of filament, falls off
if cb>params['fil_length']:
lx = Z[0,lfil];
ly = Z[1,lfil];
lfil = -1;
rfil = -1;
# motor continues to travel along filament
else:
s = np.sqrt((Z[0,lfil] - lx)**2 + (Z[1,lfil] - ly)**2) - (params['v']*params['tstep']);
lx = Z[0,lfil] + (s*np.cos(Z[2,lfil]));
ly = Z[1,lfil] + (s*np.sin(Z[2,lfil]));
# update location of bundled motor
rx = lx;
ry = ly;
return lfil,rfil,lx,ly,rx,ry
def depolymerization(A,J,fil,params,boundary):
# remove all bound motors from filament
for j in range(0,params['motorsX']):
# check if motor is active and bundled
if J[0,j] == 2:
# if motor j is on filament, but motor bundled to j is not
if J[2,j]==fil and J[2,J[1,j].astype(int)]!=fil:
# mark motor j as unbound and move it to position of bundled motor
J[2,j] = -1;
J[3,j] = J[3,J[1,j].astype(int)];
J[4,j] = J[4,J[1,j].astype(int)];
# if motor bundled to j is on filament, but motor j is not
elif J[2,j]!=fil and J[2,J[1,j].astype(int)]==fil:
# mark bundled motor as unbound and move it to position of motor j
J[2,J[1,j].astype(int)] = -1;
J[3,J[1,j].astype(int)] = J[3,j];
J[4,J[1,j].astype(int)] = J[4,j];
# if both bundled motors are on filament
elif J[2,j]==fil and J[2,J[1,j].astype(int)]==fil:
# mark both bundled motors as unbound
J[2,j] = -1;
J[2,J[1,j].astype(int)] = -1;
# equal probability of one motor to move to position of the other
if np.random.uniform()>0.5:
J[3,j] = J[3,J[1,j].astype(int)];
J[4,j] = J[4,J[1,j].astype(int)];
else:
J[3,J[1,j].astype(int)] = J[3,j];
J[4,J[1,j].astype(int)] = J[4,j];
for a in range(0,params['motorsI']):
# check if motor a is activate and bound to IgG and if motor a is on filament
if A[2,a]==fil:
A[2,a] = -1
if A[3,a]==2:
A[0,a] = A[4,a]
A[1,a] = A[5,a]
# reset position of filament (to keep number of filaments constant
# by assuming for each filament deploymerized, one is polymerized)
x_plus_end = 2*params['boundary_width']*np.random.uniform()-params['boundary_width'];
y_plus_end = 2*params['boundary_height']*np.random.uniform()-params['boundary_height'];
while 1:
angle = 2*np.pi*np.random.uniform();
x_minus_end = x_plus_end+params['fil_length']*np.cos(angle);
y_minus_end = y_plus_end+params['fil_length']*np.sin(angle);
# check if minus ends of filament are within boundary
if boundary.contains_point(np.array([x_minus_end,y_minus_end])):
break
fil_params = [x_plus_end,y_plus_end,angle,x_minus_end,y_minus_end];
return A,J,x_plus_end,y_plus_end,angle,x_minus_end,y_minus_end,0
def transplant_plus(x,y,ctr,rad):
if not (ctr[0]==x and ctr[1]==y): # point not at center
if x == ctr[0]: # infinite slope
if y>ctr[1]:
xnew = x
ynew = rad
elif y<ctr[1]:
xnew = x
ynew = -rad
elif y == ctr[1]: # zero slope
if x>ctr[0]:
xnew = rad
ynew = y
elif x<ctr[0]:
xnew = -rad
ynew = y
else: # not infinite or zero slope
y_dist = (y-ctr[1])
x_dist = (x-ctr[0])
slope = y_dist/x_dist
theta = np.arctan(slope)
if x_dist>0 and y_dist>0: # point in quadrant I
xnew = rad * np.cos(theta)
ynew = rad * np.sin(theta)
elif x_dist<0 and y_dist>0: # point in quadrant II
xnew = -rad * np.cos(theta)
ynew = -rad * np.sin(theta)
elif x_dist<0 and y_dist<0: # point in quadrant III
xnew = -rad * np.cos(theta)
ynew = -rad * np.sin(theta)
elif x_dist>0 and y_dist<0: # point in quadrant IV
xnew = rad * np.cos(theta)
ynew = rad * np.sin(theta)
return xnew+ctr[0],ynew+ctr[1]
else:
# point at center of circle, choose random point on circumference
theta = np.random.rand()*2*np.pi
xnew = rad*np.cos(theta)
ynew = rad*np.sin(theta)
def filament_boundary_check(Z,f,params,boundary):
xp = Z[0,f];
yp = Z[1,f];
angle = Z[2,f];
xm = Z[3,f];
ym = Z[4,f];
h = params['boundary_height'];
w = params['boundary_width'];
# check if plus end of filament is within boundary
if ~boundary.contains_point(np.array([xp,yp])):
b = [np.less(yp,-1*h),np.greater(yp,h),np.less(xp,-1*w),np.greater(xp,w)];
# b[0] = yp position below bottom boundary
# b[1] = yp position above top boundary
# b[2] = xp position left of left boundary
# b[3] = xp position right of right boundary
if b[0] == True and b[2] == False and b[3] == False:
yp = -1*h;
ym = yp + params['fil_length']*np.sin(angle);
elif b[1] == True and b[2] == False and b[3] == False:
yp = h;
ym = yp + params['fil_length']*np.sin(angle);
elif b[2] == True and b[0] == False and b[1] == False:
xp = -1*w;
xm = xp + params['fil_length']*np.cos(angle);
elif b[3] == True and b[0] == False and b[1] == False:
xp = w;
xm = xp + params['fil_length']*np.cos(angle);
elif b[0] == True and b[2] == True:
xp = -1*w;
yp = -1*h;
xm = xp + params['fil_length']*np.cos(angle);
ym = yp + params['fil_length']*np.sin(angle);
elif b[0] == True and b[3] == True:
xp = w;
yp = -1*h;
xm = xp + params['fil_length']*np.cos(angle);
ym = yp + params['fil_length']*np.sin(angle);
elif b[1] == True and b[2] == True:
xp = -1*w;
yp = h;
xm = xp + params['fil_length']*np.cos(angle);
ym = yp + params['fil_length']*np.sin(angle);
elif b[1] == True and b[3] == True:
xp = w;
yp = h;
xm = xp + params['fil_length']*np.cos(angle);
ym = yp + params['fil_length']*np.sin(angle);
# check if minus end of filament is within boundary
if ~boundary.contains_point(np.array([xm,ym])):
b = [np.less(ym,-1*h),np.greater(ym,h),np.less(xm,-1*w),np.greater(xm,w)];
# b[0] = ym position below bottom boundary
# b[1] = ym position above top boundary
# b[2] = xm position left of left boundary
# b[3] = xm position right of right boundary
if b[0] == True and b[2] == False and b[3] == False:
ym = -1*h;
yp = ym - params['fil_length']*np.sin(angle);
elif b[1] == True and b[2] == False and b[3] == False:
ym = h;
yp = ym - params['fil_length']*np.sin(angle);
elif b[2] == True and b[0] == False and b[1] == False:
xm = -1*w;
xp = xm - params['fil_length']*np.cos(angle);
elif b[3] == True and b[0] == False and b[1] == False:
xm = w;
xp = xm - params['fil_length']*np.cos(angle);
elif b[0] == True and b[2] == True:
xm = -1*w;
ym = -1*h;
xp = xm - params['fil_length']*np.cos(angle);
yp = ym - params['fil_length']*np.sin(angle);
elif b[0] == True and b[3] == True:
xm = w;
ym = -1*h;
xp = xm - params['fil_length']*np.cos(angle);
yp = ym - params['fil_length']*np.sin(angle);
elif b[1] == True and b[2] == True:
xm = -1*w;
ym = h;
xp = xm - params['fil_length']*np.cos(angle);
yp = ym - params['fil_length']*np.sin(angle);
elif b[1] == True and b[3] == True:
xm = w;
ym = h;
xp = xm - params['fil_length']*np.cos(angle);
yp = ym - params['fil_length']*np.sin(angle);
return [xp,yp,angle,xm,ym]
def force_on_filament(A,J,xp,yp,angle,xm,ym,igg_bound,fil,params,check):
# frictional coefficients
eta = 0.6; # pN*s/um^2 or Pa*s
di = 0.008; # diameter of actin filament in um
length = params['fil_length'];
p = length/di;
gparl = -0.2;
gperp = 0.84;
grot = -0.662;
l_pr = (2*np.pi*eta*length)/(np.log(p) + gparl);
l_pp = (4*np.pi*eta*length)/(np.log(p) + gperp);
l_r = ((1/3)*np.pi*eta*(length**3))/(np.log(p) + grot);
# thermal diffusion of filament paramters
KT = 0.005;
Dpr = KT/l_pr;
Dpp = KT/l_pp;
Dr = KT/l_r;
# pre_allocate force vector component arrays and torque arrays
force_vectorsX = np.zeros((2,params['motorsX']));
torquesX = np.zeros(params['motorsX']);
force_vectorsI = np.zeros((2,params['motorsI']));
torquesI = np.zeros(params['motorsI']);
# rotation matrix and inverse
rot = np.array([[np.cos(angle),np.sin(angle)],[-1*np.sin(angle),np.cos(angle)]]);
rot_inv = np.array([[np.cos(angle),-1*np.sin(angle)],[np.sin(angle),np.cos(angle)]]);
# center of mass position of filament
COMx = (xp + xm)/2;
COMy = (yp + ym)/2;
COM = np.array([[COMx],[COMy]]);
# center of mass position in rotated coordinate system; "matmul = matrix multiplication"
COM_rot = np.matmul(rot,COM);
# Calculate forces from myosin X
# every motor on the filament (group1)
v1 = np.where(J[2,:]==fil)[0];
# bundle index for each motor in v1 (group2)
v2 = J[1,v1].astype(int)
# bound filament index for v2 motors (group3)
v3 = J[2,v2];
# motors of group3 that are not bound to a filament
r1 = np.where(v3==-1)[0];
# motor bundles where both legs are bound to same filament
r2 = np.where(fil==v3)[0];
# group1 with exclusions (r1 and r2)
g1 = np.delete(v1,np.union1d(r1,r2));
# position of first motor
x1X = J[3,g1];
y1X = J[4,g1];
# position of second motor
x2X = J[3,J[1,g1].astype(int)];
y2X = J[4,J[1,g1].astype(int)];
# compute x and y components of force vector
F_x_X = params['k']*(x2X - x1X);
F_y_X = params['k']*(y2X - y1X);
# compute length of motor from center of mass of filament
leverX = np.sqrt((x1X - COMx)**2 + (y1X - COMy)**2);
# check which side of COM the motor is
for motX in range(0,len(leverX)):
if np.sqrt((xp - x1X[motX])**2 + (yp - y1X[motX])**2)<length/2:
leverX[motX] = -1*leverX[motX];
# rotate x and y components into parallel and perpendicular components
# then insert parallel and perpindicular components of force vector into array
force_vectorsX[0,g1] = rot[0,0]*F_x_X + rot[0,1]*F_y_X;
force_vectorsX[1,g1] = rot[1,0]*F_x_X + rot[1,1]*F_y_X;
torquesX[g1] = leverX*force_vectorsX[1,g1];
# Calculate forces from myosin I
# every myosin I motor bound to both the filament and the IgG domain
a1 = np.where(A[3,:]==2)[0]
a2 = np.where(A[2,:]==fil)[0]
g2 = np.intersect1d(a1,a2)
# positions of myosin I motors' head and tail domains
x1I = A[0,g2] # head
y1I = A[1,g2] # head
x2I = A[4,g2] # tail
y2I = A[5,g2] # tail
# force vectors from myosinI motors
F_x_I = params['k_I']*(x2I - x1I)
F_y_I = params['k_I']*(y2I - y1I)
# compute length of myosin I motor from center of mass of filament
leverI = np.sqrt((x1I - COMx)**2 + (y1I - COMy)**2);
# check which side of COM the myosin I motor is
for motI in range(0,len(leverI)):
if np.sqrt((xp - x1I[motI])**2 + (yp - y1I[motI])**2)>length/2:
leverI[motI] = -1*leverI[motI];
# rotate x and y components into parallel and perpendicular components
# then insert parallel and perpindicular components of force vector into array
force_vectorsI[0,g2] = rot[0,0]*F_x_I + rot[0,1]*F_y_I;
force_vectorsI[1,g2] = rot[1,0]*F_x_I + rot[1,1]*F_y_I;
torquesI[g2] = leverI*force_vectorsI[1,g2];
# pre-allocate array for new COM positions in rotated coordinate system
COM_rot_new = np.zeros([2,1]);
if igg_bound != 1:
# new parallel component of COM of filament
COM_rot_new[0] = COM_rot[0] + \
params['tstep']*np.sum(force_vectorsX[0,:])/l_pr + \
params['tstep']*np.sum(force_vectorsI[0,:])/l_pr + \
np.sqrt(2*params['tstep']*Dpr)*np.random.randn()
# new perpendicular component of COM of filament
COM_rot_new[1] = COM_rot[1] + \
params['tstep']*np.sum(force_vectorsX[1,:])/l_pp + \
params['tstep']*np.sum(force_vectorsI[1,:])/l_pp + \
np.sqrt(2*params['tstep']*Dpp)*np.random.randn()
# rotate COM positions in rotated cooridnate system back to standard coordinate system
COM_new = np.matmul(rot_inv,COM_rot_new);
# new angular component of filament
angle_new = angle + \
params['tstep']*np.sum(torquesX)/l_r + \
params['tstep']*np.sum(torquesI)/l_r + \
np.sqrt(2*params['tstep']*Dr)*np.random.randn()
# update positions of plus end of filament
xp_new = COM_new[0] - (length/2)*np.cos(angle_new)
if check == 1:
xp_new = xp_new + (0.004 * np.random.randn() + 0.002)
yp_new = COM_new[1] - (length/2)*np.sin(angle_new)
else:
##**** HOW TO INCORPORATE LEVER PARAMTER FOR ACTIN W/ BOUND MOTOR AND ON IGG DOMAIN? ****##
angle_new = angle + \
params['tstep']*np.sum(torquesX)/l_r + \
params['tstep']*np.sum(torquesI)/l_r + \
np.sqrt(2*params['tstep']*Dr)*np.random.randn()
xp_new = xp
yp_new = yp
return xp_new,yp_new,angle_new
# image builders
def image_build(A,J,Z,params):
r = params['res']
s = params['s']
h = params['boundary_height']
w = params['boundary_width']
# configure data for image construction
Zpx = (Z[0,:]+(w))*((s[0]/w/2)*r)
Zpy = (Z[1,:]+(h))*((s[1]/h/2)*r)
Zmx = (Z[3,:]+(w))*((s[0]/w/2)*r)
Zmy = (Z[4,:]+(h))*((s[1]/h/2)*r)
Jx1 = (J[3,:]+(w))*((s[0]/w/2)*r)
Jy1 = (J[4,:]+(h))*((s[1]/h/2)*r)
Jx2 = (J[3,J[1,:].astype(int)]+(w))*((s[0]/w/2)*r)
Jy2 = (J[4,J[1,:].astype(int)]+(h))*((s[1]/h/2)*r)
Ax1 = (A[0,:]+(w))*((s[0]/w/2)*r)
Ay1 = (A[1,:]+(h))*((s[1]/h/2)*r)
Ax2 = (A[4,:]+(w))*((s[0]/w/2)*r)
Ay2 = (A[5,:]+(h))*((s[1]/h/2)*r)
# create new PIL image object
im = Image.new('RGB', s)
draw = ImageDraw.Draw(im)
# draw IgG domain
ctrx = (params['igg_center'][0]+(w))*((s[0]/w/2)*r)
ctry = (params['igg_center'][1]+(h))*((s[1]/h/2)*r)
radx = (params['igg_radius']*((s[0]/w/2)*r))
rady = (params['igg_radius']*((s[1]/h/2)*r))
draw.ellipse([ctrx-radx,ctry-rady,ctrx+radx,ctry+rady],fill=(80,80,80)) # dark gray
# plot filament positions
igg_off = np.where(Z[5,:]==0)[0]
[draw.line(((Zpx[i],Zpy[i]),(Zmx[i],Zmy[i])),fill=(255,0,0)) for i in igg_off] # red
# igg_on = np.where(Z[5,:]==1)[0]
igg_on = A[2,np.where(A[2,:]!=-1)[0]]
[draw.line(((Zpx[i],Zpy[i]),(Zmx[i],Zmy[i])),fill=(255,0,255)) for i in igg_on] # magenta
# plot myosin I motors that are bound to the IgG domain
igg_on_I = np.where(A[3,:]==2)[0]
[draw.line(((Ax1[i],Ay1[i]),(Ax2[i],Ay2[i])),fill=(0,255,255)) for i in igg_on_I] # cyan
# plot non-IgG-bound myosin I motor positions
igg_off_I = np.where(A[3,:]!=2)[0]
[draw.point([(Ax1[i],Ay1[i])],fill=(255,122,0)) for i in igg_off_I] # orange
# find all bundled pairs of myosin X motors
motorsX = (np.where(J[1,:]!=-1)[0])
# plot motorX bundles
[draw.line(((Jx1[i],Jy1[i]),(Jx2[i],Jy2[i])),fill=(0,255,0)) for i in motorsX] # green
# plot myosin X motor positions
[draw.point([(Jx1[i],Jy1[i])],fill=(255,255,255)) for i in range(params['motorsX'])] # white
return np.asarray(im);
def image_build_z(A,J,Z,params):
r = params['res']
s = params['s']
h = params['boundary_height']
w = params['boundary_width']
# configure data for image construction
Zpx = (Z[0,:]+(w))*((s[0]/w/2)*r)
Zpy = (Z[1,:]+(h))*((s[1]/h/2)*r)
Zmx = (Z[3,:]+(w))*((s[0]/w/2)*r)
Zmy = (Z[4,:]+(h))*((s[1]/h/2)*r)
Jx1 = (J[3,:]+(w))*((s[0]/w/2)*r)
Jy1 = (J[4,:]+(h))*((s[1]/h/2)*r)
Jx2 = (J[3,J[1,:].astype(int)]+(w))*((s[0]/w/2)*r)
Jy2 = (J[4,J[1,:].astype(int)]+(h))*((s[1]/h/2)*r)
Ax1 = (A[0,:]+(w))*((s[0]/w/2)*r)
Ay1 = (A[1,:]+(h))*((s[1]/h/2)*r)
Ax2 = (A[4,:]+(w))*((s[0]/w/2)*r)
Ay2 = (A[5,:]+(h))*((s[1]/h/2)*r)
# create new PIL image object
im = Image.new('RGB', s)
draw = ImageDraw.Draw(im)
# draw IgG domain
ctrx = (params['igg_center'][0]+(w))*((s[0]/w/2)*r)
ctry = (params['igg_center'][1]+(h))*((s[1]/h/2)*r)
radx = (params['igg_radius']*((s[0]/w/2)*r))
rady = (params['igg_radius']*((s[1]/h/2)*r))
draw.ellipse([ctrx-radx,ctry-rady,ctrx+radx,ctry+rady],fill=(80,80,80)) # dark gray
# plot filament positions
[draw.line(((Zpx[i],Zpy[i]),(Zmx[i],Zmy[i])),fill=(255,0,0)) for i in range(params['filaments'])] # red
# plot myosin X motor positions
[draw.point([(Jx1[i],Jy1[i])],fill=(255,255,255)) for i in range(params['motorsX'])] # white
return np.asarray(im);
# initializers
def initialize_positions_box(params,boundary,A,J,Z,aa,jj,zz):
boundary_height = params['boundary_height'];
boundary_width = params['boundary_width'];
boundary = pltpath.Path([(-1*boundary_width,-1*boundary_height),
(-1*boundary_width,boundary_height),
(boundary_width,boundary_height),
(boundary_width,-1*boundary_height)],
readonly=True)
if zz == 1:
# assign random positions to ends of filaments
for z in range(0,params['filaments']):
# randomly assign positions to plus ends of filaments
Z[0,z] = 2*(boundary_width/params['init_zone_factor'])*np.random.uniform()-boundary_width;
Z[1,z] = 2*boundary_height*np.random.uniform()-boundary_height;
while 1:
angle = 2*np.pi*np.random.uniform();
x_minus_end = Z[0,z]+params['fil_length']*np.cos(angle);
y_minus_end = Z[1,z]+params['fil_length']*np.sin(angle);
# check if minus ends of filament are within boundary
if boundary.contains_point(np.array([x_minus_end,y_minus_end])):
Z[2,z] = angle;
Z[3,z] = x_minus_end;
Z[4,z] = y_minus_end;
break
if jj == 1:
# assign random positions to myosinX motors
for j in range(0,params['motorsX']):
J[3,j] = 2*(boundary_width/params['init_zone_factor'])*np.random.uniform()-boundary_width;
J[4,j] = 2*boundary_height*np.random.uniform()-boundary_height;
if aa == 1:
# assign random positions to myosinI motors
for a in range(0,params['motorsI']):
A[0,a] = 2*(boundary_width/params['init_zone_factor'])*np.random.uniform()-boundary_width;
A[1,a] = 2*boundary_height*np.random.uniform()-boundary_height;
return A,J,Z
def initialize_positions_ring(params,boundary,loc,A,J,Z,aa,jj,zz):
min_r_z = loc[0]
max_r_z = loc[1]
min_r_j = loc[2]
max_r_j = loc[3]
min_r_a = loc[4]
max_r_a = loc[5]
if zz == 2:
# assign random positions to ends of filaments
for z in range(0,params['filaments']):
while 1:
# randomly assign positions to plus ends of filaments
Z[0,z] = 2*(boundary_width/params['init_zone_factor'])*np.random.uniform()-boundary_width;
Z[1,z] = 2*boundary_height*np.random.uniform()-boundary_height;
if np.sqrt((Z[0,z]-params['igg_center'][0])**2 + (Z[1,z]-params['igg_center'][1])**2)>params['igg_radius']*min_r_z and \
np.sqrt((Z[0,z]-params['igg_center'][0])**2 + (Z[1,z]-params['igg_center'][1])**2)<params['igg_radius']*max_r_z:
break
while 1:
angle = 2*np.pi*np.random.uniform();
x_minus_end = Z[0,z]+params['fil_length']*np.cos(angle);
y_minus_end = Z[1,z]+params['fil_length']*np.sin(angle);
# check if minus ends of filament are within boundary
if boundary.contains_point(np.array([x_minus_end,y_minus_end])) and \
np.sqrt((x_minus_end-params['igg_center'][0])**2 + (y_minus_end-params['igg_center'][1])**2)>params['igg_radius']*(min_r_z) and \
np.sqrt((x_minus_end-params['igg_center'][0])**2 + (y_minus_end-params['igg_center'][1])**2)<params['igg_radius']*(max_r_z):
Z[2,z] = angle;
Z[3,z] = x_minus_end;
Z[4,z] = y_minus_end;
break
if jj == 2:
# assign random positions to myosin X motors
for j in range(0,params['motorsX']):
while 1:
J[3,j] = 2*(boundary_width/params['init_zone_factor'])*np.random.uniform()-boundary_width;
J[4,j] = 2*boundary_height*np.random.uniform()-boundary_height;
if np.sqrt((J[3,j]-params['igg_center'][0])**2 + (J[4,j]-params['igg_center'][1])**2)>params['igg_radius']*min_r_j and \
np.sqrt((J[3,j]-params['igg_center'][0])**2 + (J[4,j]-params['igg_center'][1])**2)<params['igg_radius']*max_r_j:
break
if aa == 2:
# assign random positions to myosin I motors
for a in range(0,params['motorsI']):
while 1:
A[0,a] = 2*(boundary_width/params['init_zone_factor'])*np.random.uniform()-boundary_width;
A[1,a] = 2*boundary_height*np.random.uniform()-boundary_height;
if np.sqrt((A[0,a]-params['igg_center'][0])**2 + (A[1,a]-params['igg_center'][1])**2)>params['igg_radius']*min_r_a and \
np.sqrt((A[0,a]-params['igg_center'][0])**2 + (A[1,a]-params['igg_center'][1])**2)<params['igg_radius']*max_r_a:
break
return A,J,Z
In [ ]:
# DEFINE PARAMETERS HERE
# define dictionary to store all parameters
params = {
# image parameters
'save_image':1, #0=No, 1=Yes
'res':1, #set resolution (r=1 == s[0] x s[1] pixels)
's':(300,300), # size of image in pixels
'image_iter':1,
'image_save_path':"./sim_ims/",
'tiff_save_name':"sim_.tif",
# simulation parameters
'time':1000, #iterations of simulation
'tstep':0.05, #time step between each iteration
'filaments':0, #number of filaments
'motorsX':0, #number of myosin X motors
'motorsI':0, #number of myosin I motors
# filament parameters
'fil_length':0.3, #filament length um
'p0':0.1, #myosinX detachment rate
'p1':1.0, #myosinX attachement rate
'p2':0.1, #myosinI detachment rate (5.6)
'p3':1.0, #myosinI attachement rate (0.7)
'p4':0.0, #depolymerization (turnover) rate
'igg_attach_rate':0.0, #IgG binding rate for plus ends of filaments to IgG domain
'igg_detach_rate':0.0, #IgG detachment rate for plus ends of filaments from IgG domain
# myosinX motor parameters
'r':0.2, #radius of search for motor attachment to filament & max stretch of bundle um (.066)
'v':0.2, #velocity of myosinX motor um/s (.2)
'k':2.0, #stiffness of dimer 'spring' nN/um or pN/nm (2.0)
'diffusionRate':0.0029,
'unbundleRate':0.2,
'bundleRate':0.8,
'activationRate':3.0,
'inactivationRate':1.0,
'motorBundleRadius':0.2, #same as r according to literature for myosinII um
# myosinI motor parameters
'r_I':0.1, #radius of parameter to find filament once bound to IgG domain um
'v_I':0.15, #velocity of myoIB motor um/s
'k_I':2.0, #stiffness of neck 'spring' nN/um or pN/nm
'myoI_igg_on':2.0, #binding rate to IgG domain (membrane)
'myoI_igg_off':0.1, #unbinding rate from IgG domain (membrane)
'myoI_activationRate':3.0,
'myoI_inactivationRate':1.0,
# edge boundary parameters
'boundary_width':7, #um 9
'boundary_height':7, #um 9
# IgG domain parameters
'igg_radius':5.0, #um 7.5
'igg_center':[0.0,0.0], #um from center at [0,0]
# initial zoning parameter
'init_zone_factor':1.0
}
# # check if save path exists; if not, create
# if params['save_image'] == 1:
# if not path.exists(params['image_save_path']):
# makedirs(params['image_save_path'])
# J -> myosinX data:
# J(0,:) = progression of myosinX state
# 0-inactive
# 1-active
# 2-bundled & active
# J(1,:) = index of motor bundled to (-1 means no bundle)
# J(2,:) = index of filament bound to (-1 means no binding to filament)
# J(3,:) = x position of motor head
# J(4,:) = y position of motor head
# A -> myosinI data:
# A(0,:) = x position of motor head
# A(1,:) = y position of motor head
# A(2,:) = index of filament bound to (-1 means no binding to filament)
# A(3,:) = progression of myosinI state
# 0-inactive
# 1-active
# 2-bound to IgG domain & active
# A(4,:) = x position of motor tail bound to IgG domain
# A(5,:) = y position of motor tail bound to IgG domain
# Z -> filament data:
# Z(0,:) = x position of filament +end
# Z(1,:) = y position of filament +end
# Z(2,:) = angle of orientation (unit circle, CCW)
# Z(3,:) = Z(1,:) - (fil_length * cos(Z(3,:))) x position of filament -end
# Z(4,:) = Z(2,:) - (fil_length * sin(Z(3,:))) y position of filament -end
# Z(5,:) = indicator of binding to IgG domain
In [ ]:
# INITIALIZATION
from PIL import Image, ImageDraw
from os import path, makedirs
def init(params):
# fix saved image size to match specified resolution
params['s'] = np.asarray((params['s'][0]*params['res'],params['s'][1]*params['res']))
# pre-allocate memory for output arrays
A = np.zeros((6,params['motorsI'])) #A == 1st leter, myosinI
A[[2,4,5],:] = -1
J = np.zeros((5,params['motorsX'])) #J == 10th letter, myosinX
J[[1,2],:] = -1
Z = np.zeros((6,params['filaments']))
# define boundary
boundary_height = params['boundary_height'];
boundary_width = params['boundary_width'];
boundary = pltpath.Path([(-1*boundary_width,-1*boundary_height),
(-1*boundary_width,boundary_height),
(boundary_width,boundary_height),
(boundary_width,-1*boundary_height)],
readonly=True)
# initialize positions
# 1 = box; 2 = ring
aa = 2 # initliaize boundary shape for myosin I
jj = 2 # initliaize boundary shape for myosin X
zz = 2 # initliaize boundary shape for actin filaments
if aa == 1 or jj == 1 or zz == 1:
A,J,Z = initialize_positions_box(params,boundary,A,J,Z,aa,jj,zz)
if aa == 2 or jj == 2 or zz == 2:
loc = np.zeros(6)
# min radius for actin
loc[0] = 0.98
# max radius for actin
loc[1] = 1.02
# min radius for myosin X
loc[2] = 1.0
# max radius for myosin X
loc[3] = 1.02
# min radius for myosin I
loc[4] = 0.99
# max radius for myosin I
loc[5] = 1.005
A,J,Z = initialize_positions_ring(params,boundary,loc,A,J,Z,aa,jj,zz)
strZ = str(params['filaments'])
strJ = str(params['motorsX'])
if len(strJ) == 1:
strJ = '000'+strJ
if len(strJ) == 3:
strJ = '0'+strJ
strA = str(params['motorsI'])
if len(strA) == 1:
strA = '000'+strA
if len(strA) == 3:
strA = '0'+strA
# write image data to file
if params['save_image']==1:
order = len(str(params['time']))
lead_zeros = order - 1
built_name = ''
for i in range(lead_zeros):
built_name+='0'
# name = params['image_save_path']+'sim_test'+built_name+'0'+'.png'
name = params['image_save_path']+'sim_z='+strZ+'_j='+strJ+'_a='+ \
strA+'_k=2.0_tstep=0.05_t='+built_name+'0.png'
im_png = Image.fromarray(image_build(A,J,Z,params))
im_png.save(name)
# imshow(im_png)
return A,J,Z,boundary
In [ ]:
# SIMULATION HEAD
from PIL import Image, ImageDraw
from tqdm import tqdm
def sim(A,J,Z,params,boundary):
advec_check = 0 #1 = adevection induced, 0 = no advection induced
strZ = str(params['filaments'])
strJ = str(params['motorsX'])
if len(strJ) == 1:
strJ = '000'+strJ
if len(strJ) == 3:
strJ = '0'+strJ
strA = str(params['motorsI'])
if len(strA) == 1:
strA = '000'+strA
if len(strA) == 3:
strA = '0'+strA
for t in range(0,params['time']):
# determine myosinX motor state
# myosinX motor state forward 0->1->2
# myosinX motor state backward 2->1->0
tempJ = J
# apply forward change of states on myosinX motors
for j in range(0,params['motorsX']):
if tempJ[0,j] != 2:
tempJ[0,j],tempJ[1,j],tempJ[3,j],tempJ[4,j] = motorX_state_forward(j,tempJ,params,boundary,advec_check);
if tempJ[0,j] == 2 and tempJ[0,tempJ[1,j].astype(int)] != 2:
tempJ[0,tempJ[1,j].astype(int)] = 2
tempJ[1,tempJ[1,j].astype(int)] = j
# move newly paired myosinX motor bundles to the same location, halfway between original locations
if tempJ[0,j] == 2 and J[0,j] != 2 and tempJ[2,j] == -1 and tempJ[2,tempJ[1,j].astype(int)] == -1 and \
j<tempJ[1,j]:
tempJ[3,j] = (tempJ[3,j] + tempJ[3,tempJ[1,j].astype(int)])/2
tempJ[4,j] = (tempJ[4,j] + tempJ[4,tempJ[1,j].astype(int)])/2
tempJ[3,tempJ[1,j].astype(int)] = tempJ[3,j]
tempJ[4,tempJ[1,j].astype(int)] = tempJ[4,j]
# apply backward change of states on myosinX motors
for j in range(0,params['motorsX']):
if tempJ[0,j] == 1 and params['inactivationRate']*params['tstep']>np.random.uniform():
tempJ[0,j] = 0
elif tempJ[0,j] == 2 and tempJ[2,j] == -1 and tempJ[2,tempJ[1,j].astype(int)] == -1 and \
params['unbundleRate']*params['tstep']>np.random.uniform():
tempJ[0,tempJ[1,j].astype(int)] = 1
tempJ[1,tempJ[1,j].astype(int)] = -1
tempJ[0,j] = 1
tempJ[1,j] = -1
J = tempJ
# determine myosinI state
# myosinI motor state forward 0->1->2
# myosinI motor state backward 2->1->0
# apply forward change of states on myosinI
for a in range(0,params['motorsI']):
if A[3,a] == 0 and params['myoI_activationRate']*params['tstep']>np.random.uniform() and \
np.sqrt((A[0,a]-params['igg_center'][0])**2 + (A[1,a]-params['igg_center'][1])**2)<params['igg_radius']:
A[3,a] = 1
elif A[3,a] == 1 and params['myoI_igg_on']*params['tstep']>np.random.uniform() and \
np.sqrt((A[0,a]-params['igg_center'][0])**2 + (A[1,a]-params['igg_center'][1])**2)<params['igg_radius']*(1+params['r_I']):
A[3,a] = 2
# if already bound to a filament
if A[2,a] >- 1:
# attach tail domain of myosin to closest point on IgG domain's surface
A[4,a],A[5,a] = transplant_plus(A[0,a],A[1,a],np.asarray(params['igg_center']),params['igg_radius'])
# not yet attached to a filament
else:
# move motor to closest point on IgG domain's surface
A[0,a],A[1,a] = transplant_plus(A[0,a],A[1,a],np.asarray(params['igg_center']),params['igg_radius'])
A[4,a] = A[0,a]
A[5,a] = A[1,a]
if A[2,a] == -1 and A[3,a] != 2:
# diffusion
A[0,a] = A[0,a] + np.sqrt(2*params['tstep']*params['diffusionRate'])*np.random.randn()
if advec_check == 1:
A[0,a] = A[0,a] + (0.006 * np.random.randn() + 0.003)
A[1,a] = A[1,a] + np.sqrt(2*params['tstep']*params['diffusionRate'])*np.random.randn()
A[0,a],A[1,a] = motor_boundary_check(A[0,a],A[1,a],params,boundary)
# apply backward change of states to myosinI
for a in range(0,params['motorsI']):
if A[3,a] == 1 and params['myoI_inactivationRate']*params['tstep']>np.random.uniform():
A[2,a] = -1
A[3,a] = 0
elif A[3,a] == 2 and A[2,a] == -1 and params['myoI_igg_off']*params['tstep']>np.random.uniform():
A[3,a] = 1
A[[4,5],a] = -1
# check if myosinX motors attach or detach from filaments
# call sub-routine if myosinX motor is bundled and has not been edited yet
cands = np.where(J[1,:]>J[1,J[1,:].astype(int)])[0]
if len(cands>0):
for j in cands:
J[2,j],J[2,J[1,j].astype(int)],J[3,j],J[4,j],J[3,J[1,j].astype(int)],J[4,J[1,j].astype(int)] = \
motorX_attach_detach(J[2,j].astype(int),J[2,J[1,j].astype(int)].astype(int),
J[3,j],J[4,j],
J[3,J[1,j].astype(int)],J[4,J[1,j].astype(int)],
Z,params,boundary)
# determine movement of myosinX motors that are bound to filaments
cands = np.where( ((J[2,:]!=-1) | (J[2,J[1,:].astype(int)]!=-1)) & (J[1,:]>J[1,J[1,:].astype(int)]) )[0]
# for j in range(0,params['motorsX']):
for j in cands:
# call sub-routine if myosinX motor is bundled, attached to a filament and has not been edited yet
J[2,j],J[2,J[1,j].astype(int)],J[3,j],J[4,j],J[3,J[1,j].astype(int)],J[4,J[1,j].astype(int)] = \
bound_motorX_movement(J[2,j].astype(int),J[2,J[1,j].astype(int)].astype(int),
J[3,j],J[4,j],
J[3,J[1,j].astype(int)],J[4,J[1,j].astype(int)],
Z,params)
# check if any myosinI binds/unbinds to actin filaments bound to the IgG domain
# find myosinI motors that are active
myoI_cands = np.where(A[3,:]>0)[0]
for a in myoI_cands:
# if myoI not bound to a filament, find a filament to bind to
if A[2,a] == -1:
# find filaments that are bound to IgG domain
cand_filsI = np.where((np.sqrt((Z[0,:]-A[0,a])**2 + (Z[1,:]-A[1,a])**2)<params['r_I']))[0] #(Z[5,:]==1) &
# check if binding event happens
if len(cand_filsI)*params['p3']*params['tstep']>np.random.uniform():
fil = cand_filsI[int(np.ceil(len(cand_filsI)*np.random.uniform()))-1]
A[0,a],A[1,a] = update_motor_position(Z[0,fil],Z[1,fil],Z[3,fil],Z[4,fil],A[0,a],A[1,a])
A[2,a] = fil
# if myoI bound to a filament, check if it switches to a different filament or if it unbinds
else:
# check if myosinI switches filaments it is bound to
# cand_filsI = np.where((Z[5,:]==1) & (np.sqrt((Z[0,:]-A[0,a])**2 + (Z[1,:]-A[1,a])**2)<params['r_I']))[0]
# if len(cand_filsI)*params['p3']*params['tstep']>np.random.uniform():
# A[2,a] = cand_filsI[int(np.ceil(len(cand_filsI)*np.random.uniform()))-1]
# check if myosinI unbinds to actin filament
if params['p2']*params['tstep']>np.random.uniform():
A[2,a] = -1
if A[3,a]==2:
A[0,a] = A[4,a]
A[1,a] = A[5,a]
# check if filament plus ends attach to IgG domain
# find all filaments where plus end of filament is within IgG domain
cands = np.where(np.sqrt((Z[0,:]-params['igg_center'][0])**2 + (Z[1,:]-params['igg_center'][1])**2)<params['igg_radius'])[0]
for f in cands:
# lock plus end of filament to IgG domain
if Z[5,f] == 0 and params['igg_attach_rate']>np.random.uniform():
Z[5,f] = 1
# # attach plus end of filament to closest point on IgG domain's surface
# Z[0,f],Z[1,f] = transplant_plus(Z[0,f],Z[1,f],np.asarray(params['igg_center']),params['igg_radius'])
# # update minus end of filament, retain angle orientation
# Z[3,f] = Z[0,f] - params['fil_length']*np.cos(Z[2,f]);
# Z[4,f] = Z[1,f] - params['fil_length']*np.sin(Z[2,f]);
# unlock plus end of filament from IgG domain
elif Z[5,f] == 1 and params['igg_detach_rate']>np.random.uniform():
Z[5,f] = 0
# check if any filaments undergo depolymerization
for f in range(0,params['filaments']):
# check if filament depolymerizes
if params['p4']*params['tstep']>np.random.uniform():
A,J,Z[0,f],Z[1,f],Z[2,f],Z[3,f],Z[4,f],Z[5,f] = depolymerization(A,J,f,params,boundary)
# iteratively determine filament movement
for f in range(0,params['filaments']):
# determine force and movement on filament
Z[0,f],Z[1,f],Z[2,f] = force_on_filament(A,J,Z[0,f],Z[1,f],Z[2,f],Z[3,f],Z[4,f],Z[5,f],f,params,advec_check)
Z[3,f] = Z[0,f] + params['fil_length']*np.cos(Z[2,f])
Z[4,f] = Z[1,f] + params['fil_length']*np.sin(Z[2,f])
# check if filament hit boundaries
Z[:5,f] = filament_boundary_check(Z[:5],f,params,boundary)
# move motors attached to filament accordingly
# find myosin X motors on filament
motsX1 = np.where((J[0,:]==2) & (J[0,J[1,:].astype(int)]==2) & (J[2,:]!=J[2,J[1,:].astype(int)]) & (J[2,:]!=-1) & (J[2,J[1,:].astype(int)]!=-1) & (J[2,:]==f))[0]
# distance between motor and plus end of filament
s = np.sqrt((Z[0,f] - J[3,motsX1])**2 + (Z[1,f] - J[4,motsX1])**2)
# x and y movement of motor to be on filament
J[3,motsX1] = Z[0,f] + s*np.cos(Z[2,f])
J[4,motsX1] = Z[1,f] + s*np.sin(Z[2,f])
# find myosin X motors with bundled motor on filament
motsX2 = np.where((J[0,:]==2) & (J[0,J[1,:].astype(int)]==2) & (J[2,:]!=J[2,J[1,:].astype(int)]) & (J[2,:]!=-1) & (J[2,J[1,:].astype(int)]!=-1) & (J[2,J[1,:].astype(int)]==f))[0]
# distance between motor and plus end of filament
s = np.sqrt((Z[0,f] - J[3,J[1,motsX2].astype(int)])**2 + (Z[1,f] - J[4,J[1,motsX2].astype(int)])**2)
# x and y movement of motor to be on filament
J[3,J[1,motsX2].astype(int)] = Z[0,f] + s*np.cos(Z[2,f])
J[4,J[1,motsX2].astype(int)] = Z[1,f] + s*np.sin(Z[2,f])
# find myosin I motors on filament
motsI = np.where(A[2,:]==f)
# distance between motor and plus end of filament
s = np.sqrt((Z[0,f] - A[0,motsI])**2 + (Z[1,f] - A[1,motsI])**2)
# x and y movement of motor to be on filament
A[0,motsI] = Z[0,f] + s*np.cos(Z[2,f])
A[1,motsI] = Z[1,f] + s*np.sin(Z[2,f])
# check if any bundled myosin X motors have surpassed stretch radius limit
# every motor on a filament (group1)
v1 = np.where(J[2,:]!=-1)[0]
# bundle index for each myoX motor in v1 (group2)
v2 = J[1,v1].astype(int)
# bound filament index for v1 myoX motors (group3)
v3 = J[2,v1]
# bound filament index for v2 myoX motors (group4)
v4 = J[2,v2]
# myoX motors of group4 that are not bound to a filament
r1 = np.where(v4==-1)[0]
# myoX motor bundles where both legs are bound to same filament
r2 = np.where(v3==v4)[0]
# group1 with exclusions (r1 and r2)
g1 = np.delete(v1,np.union1d(r1,r2))
# find stretch distances between candidate myosin X motors
motor_stretch = np.sqrt((J[3,g1] - J[3,J[1,g1].astype(int)])**2 + (J[4,g1] - J[4,J[1,g1].astype(int)])**2)
# find indices of stretch that surpass bundled myosin X stretch limit
overstretched = np.where(motor_stretch>params['r'])[0]
# for bundled pairs passing radius, randomly select a leg to unbind from filament and move to other leg's position
for motor in overstretched:
mot = g1[motor]
if np.random.uniform()<0.5:
J[2,mot] = -1
J[3,mot] = J[3,J[1,mot].astype(int)]
J[4,mot] = J[4,J[1,mot].astype(int)]
else:
J[2,J[1,mot].astype(int)] = -1
J[3,J[1,mot].astype(int)] = J[3,mot]
J[4,J[1,mot].astype(int)] = J[4,mot]
# check if any myosin I motors have surpassed stretch radius limit
# every motor bound to IgG domain
a1 = np.where(A[3,:]==2)[0]
# every motor bound to a filament
a2 = np.where(A[2,:]>-1)[0]
# group to check stretch radius
g2 = np.intersect1d(a1,a2)
# find stretch distances between candidate myosin I motors
motor_stretch = np.sqrt((A[0,g2] - A[4,g2])**2 + (A[1,g2] - A[5,g2])**2)
# find indices of stretch that surpass bundled myosin I stretch limit
overstretched = np.where(motor_stretch>params['r_I'])[0]
# if overstretched, randomly choose unbinding by either the tail or head domain of myosin I motor
for motor in overstretched:
mot = g2[motor]
# unbind filament (head)
if np.random.uniform()<0.5:
A[0,mot] = A[4,mot]
A[1,mot] = A[5,mot]
A[2,mot] = -1
# unbind IgG domain (tail)
else:
A[3,mot] = 1
A[[4,5],mot] = -1
# check how many units have reached opposite side of computational domain
if advec_check == 1:
advection_group = len(np.where(J[3,:]>20.0)[0])
# if a fraction of units is over range, turn off advection
if advection_group > params['motorsX']*.2:
advec_check = 0
# write image data to '.tiff' file
if params['save_image']==1:
if (t+params['image_iter']+1)%params['image_iter'] == 0:
time_point = str(int(((t+1)/params['image_iter'])))
order = len(str(params['time']))
lead_zeros = order - len(str(time_point))
built_name = ''
for i in range(lead_zeros):
built_name+='0'
# name = params['image_save_path']+'sim_test'+built_name+str(time_point)+'.png'
name = params['image_save_path']+'sim_z='+strZ+'_j='+strJ+'_a='+ \
strA+'_k=2.0_tstep=0.05_t='+built_name+str(time_point)+'.png'
im_png = Image.fromarray(image_build(A,J,Z,params))
im_png.save(name)
#image of just actin
# name = params['image_save_path']+'sim_test_z'+built_name+str(time_point)+'.png'
# im_png = Image.fromarray(image_build_z(A,J,Z,params))
# im_png.save(name)
# im_tiff = image_build(J,Z,params)
# imsave(params['image_save_path']+params['tiff_save_name'],im_tiff,append=True,compress=1) # for TIFF file
# finish and pass data back
return A,J,Z
In [ ]:
# RUN THE SIMULATION
from tqdm import tqdm
numA = [0,500,1000,1500,2000,2500,3000,6000]
numJ = [0,500,1000,1500,2000,2500,3000,6000]
numZ = [1000,2000,3000]
for nZ in tqdm(numZ):
params['filaments'] = nZ
for nJ in tqdm(numJ):
params['motorsX'] = nJ
for nA in tqdm(numA):
params['motorsI'] = nA
# initialize simulation
A,J,Z,boundary = init(params)
# run simulation
A,J,Z = sim(A,J,Z,params,boundary)
# import cProfile
# cProfile.run('sim(A,J,Z,params,boundary)')
In [ ]:
# VIEW IMAGES
import glob
from IPython.html.widgets import interact
from IPython.display import Image, display
cache = {}
def get_image(key):
if not key in cache:
cache[key] = open(key,'rb').read()
return cache[key]
images = glob.glob(params['image_save_path']+'*.png')
for image in images:
get_image(image)
def pltimg(Time=0):
order = len(str(params['time']))
lead_zeros = order - len(str(Time))
built_name = ''
for i in range(lead_zeros):
built_name+='0'
filename = params['image_save_path']+'sim_test'+built_name+str(Time)+'.png'
display(Image(data=cache[filename]))
interact(pltimg, Time=(0,params['time'],params['image_iter']))
In [ ]:
# CONVERT '.PNG' IMAGES INTO A '.TIFF' STACK
path = './../../Desktop/analyze_new/'
# path = params['image_save_path']
from scipy import misc
from tqdm import tqdm
for nZ in numZ:
for nJ in numJ:
for nA in numA:
first = 0
strZ = str(nZ)
strJ = str(nJ)
if len(strJ) == 1:
strJ = '000'+strJ
if len(strJ) == 3:
strJ = '0'+strJ
strA = str(nA)
if len(strA) == 1:
strA = '000'+strA
if len(strA) == 3:
strA = '0'+strA
filename = 'sim_z='+strZ+'_j='+strJ+'_a='+strA+'_k=2.0_tstep=0.05_t=*.png'
tiff_name = 'sim_z='+strZ+'_j='+strJ+'_a='+strA+'_k=2.0_tstep=0.05.tif'
for image_path in tqdm(glob.glob("./sim_ims/"+filename)):
im = misc.imread(image_path)
if first == 0:
imsave(path+tiff_name,im,compress=1)
first = 1
else:
imsave(path+tiff_name,im,append=True,compress=1)
In [ ]:
(1.375*120 + 2.221*77 + 3.702*47)/244