-
Notifications
You must be signed in to change notification settings - Fork 0
/
0-2_fates_functions.py
168 lines (144 loc) · 6.2 KB
/
0-2_fates_functions.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
import math
import numpy as np
import scipy.special as sp
import random as rand
import pandas as pd
from numba import jit,int64,float64
import time
import scipy.optimize as sciopt
#Parameters - self-activation a, cross-repression r, decay k, noise alpha:
a = np.array([1,1])
r = np.array([1,1])
k = np.array([1,1])
alpha = np.array([0.05,0.05])
#Hill coefficients:
n = np.array([[4,4],[4,4]])
theta = np.array([[0.5,0.5],[0.5,0.5]])
#The total time and dt::
T = 20
dt = 0.005
#Other time-related variables related to these two:
Nt = int(T/dt)
sqrt_dt = np.sqrt(dt)
TimeRange = np.arange(0,T,dt)
# #This is to calculate the dynamic threshold moving average:
# timeBase = 2*int(Nt/100)
# #Time to plot utility against:
# time_trunc = TimeRange[int((timeBase/2)-1):-int(timeBase/2)]
#Let's call num_traj the number of trajectories:
num_traj = 10000
#The initial conditions:
x0 = np.zeros((2,num_traj))
#Uniform around mean:
init_bias = 0
x0[0,:] = init_bias*np.ones(num_traj)
#The threshold above which we say x_i is high:
thresh = np.array([0.8,0.8])
#The stability threshold; a trajectory has to stay above this:
stability_thresh = 0.9
#Self-activation:
@jit(nopython=True)
def activation(x,a,n,theta):
if (x>=0):
return (a*x**n)/(x**n + theta**n)
else:
return 0
#Cross-inhibition
@jit(nopython=True)
def repression(x,r,n,theta):
if (x>0):
return (r*theta**n)/(x**n + theta**n)
else:
return 0
#Bias (for external signals):
@jit(nopython=True)
def ext_bias(x1,x2,t):
# return (x1-x2,x2-x1)
temp = 0
return (temp,0)
#This solves the system fwd using simple Euler-Maruyama:
@jit(nopython=True)
def Solver(initial,a,r,k,n,theta,alpha):
final = np.empty((2,num_traj,Nt))
final[:,:,0] = initial
#Solving the system forward in time:
for i in range(num_traj):
for t in range(1,Nt):
#Equation for first species:
temp = np.sqrt(np.maximum(final[0,i,t-1],0.01))
noise = rand.normalvariate(0,alpha[0]*temp)
final[0,i,t] = final[0,i,t-1] + dt*(activation(final[0,i,t-1],a[0],n[0,0],theta[0,0]) \
+ repression(final[1,i,t-1],r[0],n[1,0],theta[1,0]) \
- k[0]*final[0,i,t-1] + ext_bias(final[0,i,t-1],final[1,i,t-1],t)[0]) \
+ sqrt_dt*noise
#Equation for second:
temp = np.sqrt(np.maximum(final[1,i,t-1],0.01))
noise = rand.normalvariate(0,alpha[1]*temp)
final[1,i,t] = final[1,i,t-1] + dt*(activation(final[1,i,t-1],a[1],n[1,1],theta[1,1]) \
+ repression(final[0,i,t-1],r[1],n[0,1],theta[0,1]) \
- k[1]*final[1,i,t-1] + ext_bias(final[0,i,t-1],final[1,i,t-1],t)[1]) \
+ sqrt_dt*noise
return final
#Classifier:
# @jit(nopython=True)
def fate_classifier(traj):
#Stability factor of trajectories:
cross_flags = np.zeros((2,num_traj))
cross_times = np.ones((2,num_traj))*Nt
for axis_idx in range(2):
#Axis crossings:
for traj_idx in range(num_traj):
if (np.size(np.where(traj[axis_idx,traj_idx]>thresh[axis_idx])[0]) != 0):
cross_flags[axis_idx,traj_idx] = 1
cross_times[axis_idx,traj_idx] = np.where(traj[axis_idx,traj_idx]>thresh[axis_idx])[0][0]
#Stability factor: after the threshold is crossed, how much time the traj spends above it:
stability_factrs = np.zeros((2,num_traj))
for axis_idx in range(2):
for traj_idx in range(num_traj):
# if (cross_flags[axis_idx,traj_idx]==0):
# stability_factrs[axis_idx,traj_idx] = -1
if (cross_flags[axis_idx,traj_idx]==1):
stability_factrs[axis_idx,traj_idx] = np.sum(traj[axis_idx,traj_idx,int(cross_times[axis_idx,traj_idx]):]>thresh[axis_idx])\
/len(traj[axis_idx,traj_idx,int(cross_times[axis_idx,traj_idx]):])
#Stability threshold - trajectories that spend more than this above the concentration threshold are considered committed:
# stability_thresh = 0.8
#Classifying fates:
fates = np.zeros((2,num_traj))
for axis_idx in range(2):
for traj_idx in range(num_traj):
if (cross_times[axis_idx,traj_idx]<=int(Nt/2) and stability_factrs[axis_idx,traj_idx]>=stability_thresh):
fates[axis_idx,traj_idx]=1
return stability_factrs,fates
def fate_fractions(fates):
#Initializing:
fate_frax = np.zeros(4)
fate_frax[0] = np.sum((fates[0]==0) & (fates[1]==0))/num_traj
fate_frax[1] = np.sum((fates[0]==1) & (fates[1]==0))/num_traj
fate_frax[2] = np.sum((fates[0]==0) & (fates[1]==1))/num_traj
fate_frax[3] = np.sum((fates[0]==1) & (fates[1]==1))/num_traj
return fate_frax
def traj_moments(traj,fates):
#Finding the fate fractions first:
fate_frax = fate_fractions(fates)
#Flags, if fate_frac is zero for a fate then no avg or std is calculated:
fate_flags = np.array([False for i in range(4)])
#Flag = True if that fate exists in the population:
for fate_idx in range(4):
if (fate_frax[fate_idx] != 0):
fate_flags[fate_idx] = True
#Average trajectories and standard deviation around them:
avg_traj = np.zeros((4,2,Nt))
std_traj = np.zeros((4,2,Nt))
if (fate_flags[0]==True):
avg_traj[0,:,:] = np.mean(traj[:,(fates[0]==0) & (fates[1]==0),:],axis=1)
std_traj[0,:,:] = np.std(traj[:,(fates[0]==0) & (fates[1]==0),:],axis=1)
if (fate_flags[1]==True):
avg_traj[1,:,:] = np.mean(traj[:,(fates[0]==1) & (fates[1]==0),:],axis=1)
std_traj[1,:,:] = np.std(traj[:,(fates[0]==1) & (fates[1]==0),:],axis=1)
if (fate_flags[2]==True):
avg_traj[2,:,:] = np.mean(traj[:,(fates[0]==0) & (fates[1]==1),:],axis=1)
std_traj[2,:,:] = np.std(traj[:,(fates[0]==0) & (fates[1]==1),:],axis=1)
if (fate_flags[3]==True):
avg_traj[3,:,:] = np.mean(traj[:,(fates[0]==1) & (fates[1]==1),:],axis=1)
std_traj[3,:,:] = np.std(traj[:,(fates[0]==1) & (fates[1]==1),:],axis=1)
return (avg_traj,std_traj)