-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathIA_1D_MORPHOGENESIS_RATES_BINDING.py
190 lines (123 loc) · 6.86 KB
/
IA_1D_MORPHOGENESIS_RATES_BINDING.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
# Imtiaz Ali, Physics Graduate Student, Univeristy Of California, Merced
# Python 3.6.2 (v3.6.2:5fd33b5)
# 12/10/2016
# 1-D RANDOM WALK
# MORPHOGENESIS
# THIS CODE GETS THE DIFFUSION, SOURCE RATE AND BINDING OF CELLS
# PART 1 DIFFUSION -> APPLY MONTE CARLO TO GET DIFFUSION IN SOLUTION
# PART 2 SOURCE -> APPLY FIXED RATE FOR THE SOURCE
# PART 3 BINDING -> DO IRRIVERSIBLE BINDING - k_off DOESN'T EXIST, CELLS CANT LET GO
# PART 4 DIFFUSION AND BINDING
# ASSUMPTIONS:
## PARTICLES/MOLECULES ARE THE SAME SO SAME DIFFUSION CONSTANTS
## ALL CELLS ARE THE SAME SO SAME k_on
# FIRST BOX -> ONLY LET PARTICLE MOVE RIGHT TO SECOND BOX => SOURCE
# SECOND BOX -> NOT ALLOWED TO MOVE PARTICLE TO FIRST BOX, IF LEFT OCCURS, DONT MOVE IT
# LAST BOX -> ONLY LET PARTICLE STAY IN LAST BOX => SINK
# SECOND TO LAST BOX -> NOT ALLOWED TO RECIEVE PARTICLE FROM LAST BOX
##############################################################################
##############################################################################
# FUTURE WORK
# 1. LOOK AT CELLS THAT ARE NOT THE SAME -> DIFFERENT K_ON
# 2. ALLOW CELLS TO RELEASE PARTICLES
# 3. ALLOW MORE THAN 1 PARTICLE TO MOVE INSTEAD OF JUST ONE
# 4. LOOK AT OARTICLES THAT ARE NOT THE SAME -> DIFFERENT DIFFUSION CONSTANT
# 5. APPLY DIFFERENT SOURCE RATES -> RELATING TO BOUNDAY CONDITIONS
# 6. LOOK AT FLUCTUATIONS IN CONCENTRATION -> MAYBE MORE THAT JUST HEAD PLAYS A ROLE HERE
# 7. SEE HOW STEADY-STATE FLUCTUATES WITH SYSTEM PERTURBATIONS
##############################################################################
##############################################################################
#!/use/bin/python
from random import choice as ch # ALIASING MODULES
import numpy as np
import matplotlib.pyplot as plt
npart = 50000 # TIME-STEPS => t = 0 TO npart-1 => TOTAL STEP = npart
side = 10 # NUMBER OF BOXES, SHOULD BE AN ODD NUMBER
nparticles_total = 100000 # TOTAL NUMBER OF PARTICLES IN SYSTEM
nparticles = 10 # NUMBER OF PARTICLES TO START OFF WITH
nparticles_left = nparticles_total - nparticles # NUMBER OF PARTICLES LEFT TO ADD TO SOURCE
steps = [-1,1] # X STEP LENGTHS
grid = np.zeros((npart+1,side))
grid[0][0] = nparticles # START ALL PARTICLES IN FIRST BOX
t_step_in = 4 # TIME-STEP WE ADD CERTAIN NUMBER OF PARTICLES TO SOURCE
rate_in = 4 # RATE WE ADD PARTICLES TO SOURCE BOX
for t in range(npart): # START AT INITIAL TIME
# ADD PARTICLES TO SOURCE EVERY CERTAIN NUMBER OF TIME-STEPS
# THIS IS NOT CORRECT, NEED PROPER CONDITIONAL FROM BOUNDAY CONDITION PDE
if t != 0 and t%t_step_in == 0:
grid[t][0] = grid[t][0] + rate_in
nparticles_left = nparticles_left - rate_in
## print(grid[t][0:side])
grid_p = grid[t][0:side]
# LOOP OVER EACH BOX
for x in range(side):
# FIRST CHECK TO SEE IF BOX IS NOT EMPTY -> NOT EMPTY IMPLIES WE CAN MOVE PARTICLE
if grid_p[x] != 0: # NOT EMPTY
# BINDING PHASE, k_on -> MONTE CARLO STEP
k_on = round(np.random.rand(1)[0],4) # 1 NUMBER ARRAY -> [0] IS POINTER TO FIRST ARRAY ELEMENT
if k_on >= 0.00 and k_on < 0.50: # CELL DOESN'T LET ALL PARTICLES IN IT TO MOVE AT ALL
grid_p[x] += 0
else:
# RNG TO MOVE LEFT OR RIGHT -> MONTE CARLO STEP
r_move = round(np.random.rand(1)[0],4) # 1 NUMBER ARRAY -> [0] IS POINTER TO FIRST ARRAY ELEMENT
## print(r_move)
###################################################################################
if r_move >= 0.00 and r_move < 0.50: # MOVE RIGHT
if x == 0: # -> PARTIALLY CORRECT
if grid_p[x] > 0:
grid_p[x] += steps[0] # LOSE 1 PARTICLE
grid_p[x+1] += steps[1] # GAIN 1 PARTICLE
if x > 0 and x < side-1: # ISSUE MAYBE
grid_p[x] += steps[0] # LOSE 1 PARTICLE
grid_p[x+1] += steps[1] # GAIN 1 PARTICLE
if x == side-1: # LAST BOX IS A SINK - NOTHING LEAVES IT -> CORRECT
grid_p[x] += 0
###################################################################################
if r_move >= 0.50 and r_move < 1: # MOVE LEFT
# NOTE, IT DOESN'T ACCOUNT FOR THINGS ALREADY COMING IN AT TIME STEP BEFORE
if x == 0: # FIRST BOX CAN ONLY MOVE RIGHT -> CORRECT
grid_p[x] += 0
if x == 1: # SECOND BOX CAN ONLY MOVE RIGHT -> PARTIALLY CORRECT
grid_p[x] += 0
if x > 1 and x < side-1: # ISSUE MAYBE
grid_p[x] += steps[0] # LOSE 1 PARTICLE
grid_p[x-1] += steps[1] # GAIN 1 PARTICLE
if x == side-1: # LAST BOX IS A SINK -> NOTHING LEAVES IT -> CORRECT
grid_p[x] += 0
grid[t+1][0:side] = grid_p[0:side]
###################################################################################
# PLOT COUNT IN EACH BOX FOR SEPERATE TIME-STEPS
# PLOT ALL CELL BOXES AND WITHOUT SOURCE AND SINK -> SUBPLOT
##figure(num=None, figsize=(8, 6), dpi=80, facecolor='w', edgecolor='k')
fig, (ax1,ax2)=plt.subplots(1,2)
for i in range(npart//5000):
ax1.plot(grid[i*5000][0:side], label ='time-step = %d' %(i*5000))
ax2.plot(grid[i*5000][1:side-1], label ='time-step = %d' %(i*5000))
ax1.set_title('Concentration vs. Cell box')
ax2.set_title("Concentration vs. Cell box")
chartBox = ax2.get_position()
ax2.set_position([chartBox.x0, chartBox.y0, chartBox.width*0.6, chartBox.height])
ax2.legend(loc='upper center', bbox_to_anchor=(1.45, 0.8), shadow=True, ncol=1)
plt.show()
# PLOT ALL CELL BOXES
fig = plt.figure()
ax = plt.subplot(111)
for i in range(npart//5000):
ax.plot(grid[i*5000][0:side], label ='time-step = %d' %(i*5000))
plt.title('Concentration vs. Cell box')
plt.xlim(0,side)
chartBox = ax.get_position()
ax.set_position([chartBox.x0, chartBox.y0, chartBox.width*0.6, chartBox.height])
ax.legend(loc='upper center', bbox_to_anchor=(1.45, 0.8), shadow=True, ncol=1)
plt.show()
# PLOT ALL CELL BOXES EXCEPT SOURCE AND SINK
fig = plt.figure()
ax = plt.subplot(111)
for i in range(npart//5000):
ax.plot(grid[i*5000][1:side-1], label ='time-step = %d' %(i*5000))
plt.title('Concentration vs. Cell box')
##plt.xlim(1,side-1)
chartBox = ax.get_position()
ax.set_position([chartBox.x0, chartBox.y0, chartBox.width*0.6, chartBox.height])
ax.legend(loc='upper center', bbox_to_anchor=(1.45, 0.8), shadow=True, ncol=1)
plt.show()