### Check out Sage Math! www.sagemath.org (I'll help you set it up)
### STOCHASTIC
### Initialize some stuff
var('t dt tf Tprint')
var('at m i')
X = [100,100,100,100,100] ## holds [species], ICs (n0,n1,n2,n3,n4)
k = [0.1,0.1,0.1,0.1,0.2,0.2,0.2,0.2,10.0] ## holds rate constants (4x k+,4x k-)
a = [0,0,0,0,0,0,0,0,0] ## holds calc'd rxn rates (once they ARE calculated)
r = [0,0,0,0,0,0,0,0,0] ## holds the random numbers (once they ARE created)
sol = [] ## holds the solution
t = 0 ## time, start
tf = 200 ## time, end
dt = 0.01 ## time, step (observational frequency)
m = len(a) ## the number of reactions under consideration
Tprint = t
### Loop 0: continues until target time ###
while (t<tf):
at = 0 ## 'overall reaction rate'; takes into account all reactions
## calculate the current reaction rates
## forward reactions:
a[0] = k[0]*X[0]
a[1] = k[1]*X[1]
a[2] = k[2]*X[2]
a[3] = k[3]*X[3]
## reverse reactions:
a[4] = k[4]*X[1]
a[5] = k[5]*X[2]
a[6] = k[6]*X[3]
a[7] = k[7]*X[4]
## irreversible activation:
a[8] = k[8]*X[4]
for j in range(0,m):
at = at + a[j] ## calculate the 'overall reaction rate'
r[j] = random() ## generate the needed random numbers
t = t+(1/at)*ln(1/r[0]) ## next rxn occurs in time increment 'tau'
while(t>=Tprint): ## until the next reaction occurs..
sol.append([Tprint,X[0],X[1],X[2],X[3],X[4]]) ## the concentrations stay the same
Tprint=Tprint+dt ## as time increments
i = 1 ## (which reaction is being investigated)
while (i <= m): ## and when a reaction DOES occur...
sum = 0
for l in range(0,i):
sum = sum + a[l] ## sum together the rates (for normalization)
if (sum > r[1]*at): ## randomly select a value (0, total rate)
if (i == 1): ## if value falls within rate #1, rxn #1 occurs
X[0] = X[0] - 1
X[1] = X[1] + 1
X[2] = X[2]
X[3] = X[3]
X[4] = X[4]
break
elif (i == 2): ## if value falls within rate #2, rxn #2 occurs
X[0] = X[0]
X[1] = X[1] - 1
X[2] = X[2] + 1
X[3] = X[3]
X[4] = X[4]
break
elif (i == 3): ## if value falls within rate #3, rxn #3 occurs
X[0] = X[0]
X[1] = X[1]
X[2] = X[2] - 1
X[3] = X[3] + 1
X[4] = X[4]
break
elif (i == 4): ## if value falls within rate #4, rxn #4 occurs
X[0] = X[0]
X[1] = X[1]
X[2] = X[2]
X[3] = X[3] - 1
X[4] = X[4] + 1
break
elif (i == 5): ## if value falls within rate #5, rxn #5 occurs
X[0] = X[0] + 1
X[1] = X[1] - 1
X[2] = X[2]
X[3] = X[3]
X[4] = X[4]
break
elif (i == 6): ## if value falls within rate #6, rxn #6 occurs
X[0] = X[0]
X[1] = X[1] + 1
X[2] = X[2] - 1
X[3] = X[3]
X[4] = X[4]
break
elif (i == 7): ## if value falls within rate #7, rxn #7 occurs
X[0] = X[0]
X[1] = X[1]
X[2] = X[2] + 1
X[3] = X[3] - 1
X[4] = X[4]
break
elif (i == 8): ## if value falls within rate #8, rxn #8 occurs
X[0] = X[0]
X[1] = X[1]
X[2] = X[2]
X[3] = X[3] + 1
X[4] = X[4] - 1
break
elif (i == 9): ## if value falls within rate #9, rxn #9 occurs
X[0] = X[0]
X[1] = X[1]
X[2] = X[2]
X[3] = X[3]
X[4] = X[4] - 1
break
else:
show("something went wrong...")
i = i + 1
#show(sol)
### Simple Plotting
sol0 = []
sol1 = []
sol2 = []
sol3 = []
sol4 = []
sol5 = []
for q in range(0,len(sol)):
sol0.append(sol[q][0])
sol1.append(sol[q][1])
sol2.append(sol[q][2])
sol3.append(sol[q][3])
sol4.append(sol[q][4])
sol5.append(sol[q][5])
a = list_plot([])
a += list_plot(zip(sol0,sol1), color='indianred')
a += list_plot(zip(sol0,sol2), color='lightsalmon')
a += list_plot(zip(sol0,sol3), color='khaki')
a += list_plot(zip(sol0,sol4), color='lightgreen')
a += list_plot(zip(sol0,sol5), color='lightblue')
#a.set_axes_range(0,100,0,200) ## (xmin,xmax,ymin,ymax)
#show(a)
|
|
### DETERMINISTIC
## This sheet simulates the dynamics of T-Cell Receptor activation
var('k1 k2 k3 k4 k5 k6 k7 k8 k9') ## kinetic parameters
var('N00 N10 N20 N30 N40') ## initial conditions
var('N0 N1 N2 N3 N4 t') ## variables
var('r1 r2 r3 r4 r5') ## reactions
## INITIAL CONDITIONS
N00 = 100 ##
N10 = 100 ##
N20 = 100 ##
N30 = 100 ##
N40 = 100 ##
## CONSTANTS
k1 = k[0]
k2 = k[1]
k3 = k[2]
k4 = k[3]
k5 = k[4]
k6 = k[5]
k7 = k[6]
k8 = k[7]
k9 = k[8]
## CALCULATION PARAMETERS
end_points = 200
stepsize = 0.1
steps = end_points/stepsize
ics = [0,N00,N10,N20,N30,N40]
dvars = [N0,N1,N2,N3,N4]
## EQUATIONS
r1 = (diff(N0,t) == k5*N1-k1*N0)
r2 = (diff(N1,t) == k1*N0+k6*N2-k2*N1-k5*N1)
r3 = (diff(N2,t) == k2*N1+k7*N3-k3*N2-k6*N2)
r4 = (diff(N3,t) == k3*N2+k8*N4-k4*N3-k7*N3)
r5 = (diff(N4,t) == k4*N3-k8*N4-k9*N4)
## NUMERICAL SOLUTION OF A SERIES OF DIFFERENTIAL EQUATIONS
des = [r1.rhs(), r2.rhs(), r3.rhs(), r4.rhs(), r5.rhs()]
sol = desolve_system_rk4(des,dvars,ics,ivar=t,end_points=end_points,step=stepsize)
## Process the output into a form that can be graphed
sols_1=[]
sols_2=[]
sols_3=[]
sols_4=[]
sols_5=[]
for i in range(steps):
sols_1.append([sol[i][0],sol[i][1]])
sols_2.append([sol[i][0],sol[i][2]])
sols_3.append([sol[i][0],sol[i][3]])
sols_4.append([sol[i][0],sol[i][4]])
sols_5.append([sol[i][0],sol[i][5]])
################################################
#### Unnecessarily Fancy Plotting Stuff ####
################################################
## Create a plot object
#a = plot([])
## Set the plot parameters
title = "TCR Activation" ## Graph Title
xmin = 0 ## X minimum
xmax = end_points ## X maximum
ymin = 0 ## Y minimum
ymax = 280 ## Y maximum
## Add a title to the plot
a += text(title,(xmax/1.8,ymax),color='black',fontsize=15)
## Add the desired lines to the plot
a += list_plot(sols_1,color='red',legend_label='N0')
a += list_plot(sols_2,color='orange',legend_label='N1')
a += list_plot(sols_3,color='yellow',legend_label='N2')
a += list_plot(sols_4,color='green',legend_label='N3')
a += list_plot(sols_5,color='blue',legend_label='N4')
## For more information on plots in general, evaluate 'plot?'
## For a list of legend options, evaluate 'a.set_legend_options?'
## For a list of Sage predefined colors, evaluate 'sorted(colors)'
## View predefined colors: http://en.wikipedia.org/wiki/X11_color_names
## Additional plot options
#a.set_axes_range(xmin,xmax,ymin,ymax)
a.axes_labels(['time','number'])
a.axes_label_color('grey')
a.set_legend_options(ncol=5,borderaxespad=5,back_color='whitesmoke',fancybox=true)
show(a)
|
|
