3607.sagews

### Simulates a system in which T-cell activation is irreversible
### Check out Sage Math! www.sagemath.org (I'll help you set it up)
### STOCHASTIC - simulated with the Gillespie algorithm

### 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([],figsize=[6,4])

## 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)