4230.sagews

## This sheet simulates two reversible first-order chemical reactions in series:
## A1 + A2 <-> A3 -> A4   (S+E ⇌ ES ➙ P+E)
## The initial value problem is solved by Runge-Kutta and the dynamics are graphed

reset()

## Initial Conditions
var('t A1 A2 A3 A4')
ivar = t
vars = [A1,A2,A3,A4]
ics = [0]           #starting at t = 0
ics.append(2.0)     #A1_0 = S0
ics.append(0.01)    #A2_0 = E0
ics.append(0.00)    #A3_0 = ES0
ics.append(0.00)    #A4_0 = P0

## Kinetic Parameters - adjust these!
k1 = .25
k_1 = .25
k2 = .5
k_2 = 0

## Equations 
rxns = []
rxns.append(diff(A1,t) == k_1*A3 - k1*A1*A2)
rxns.append(diff(A2,t) == k_1*A3 - k1*A1*A2 + k2*A3)
rxns.append(diff(A3,t) == k1*A1*A2 - k_1*A3 - k2*A3)
rxns.append(diff(A4,t) == k2*A3 - k_2*A4)

## Calculation Parameters
end_points = 3000
stepsize = 1
steps = end_points/stepsize

des = []
for i in range(0,len(rxns)):
    des.append(rxns[i].rhs())
sol = desolve_system_rk4(des,vars,ics,ivar,end_points=end_points,step=stepsize)

## Process the output into a form that can be graphed
sols=[]

for i in range(0,len(sol[0])-1):
    sols.append([])
for i in range(steps):
    for j in range(0,(len(sol[i])-1)):
        sols[j].append([sol[i][0],sol[i][j+1]])

################################################
####   Unnecessarily Fancy Plotting Stuff   ####
################################################
### 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)'

### USER-DEFINED PARAMETERS
title = "Enzymatic Reaction Kinetics"                     ## Graph Title
legend_labels = ['[S]','[E]','[ES]','[P]']               ## Legend Labels
colors = ['red','orange','yellow','green','blue']        ## Legend Colors
x_label = 'time'                                         ## X-axis label
y_label = 'concentration'                                ## Y-axis label

### Create a plot object
a = plot([])
a.axes_labels([x_label,y_label])
# a.set_axes_range(xmin=0,xmax=150,ymin=0,ymax=10)

### Add the desired lines to the plot
for i in range(0,len(sols)):
    a += list_plot(sols[i],color=colors[mod(i,len(colors))],legend_label=legend_labels[i])
    
### Set the plot parameters
a += text(title,(a.xmax()/1.8,a.ymax()),color='black',fontsize=15)
#a.axes_label_color('grey')
#a.set_legend_options(ncol=round(len(legend_labels)/2),borderaxespad=5,back_color='whitesmoke',fancybox=true)

#show(a)

b = plot([])

for i in range(1,3):
    b += list_plot(sols[i],color=colors[mod(i,len(colors))],legend_label=legend_labels[i])
    
title = "Free Enzyme and Enzyme-Substrate complex"
b += text(title,(b.xmax()/1.8,b.ymax()),color='black',fontsize=15)
b.axes_labels([x_label,y_label])

show(b)

Traceback (most recent call last): b += text(title,(b.xmax()/1.8,b.ymax()),color='black',fontsize=15) File "", line 1, in <module> File "/tmp/tmpYiprb8/___code___.py", line 6, in <module> b += list_plot(sols[i],color=colors[mod(i,len(colors))],legend_label=legend_labels[i]) IndexError: list index out of range

### Re-solve the problem using the Michaelis-Menten kinetic model
### (Briggs-Haldane derivation, applying PSSH to [ES])

## Initial Conditions
var('t S P')
ivar = t
vars = [S,P]
S0 = ics[1]
P0 = ics[4]
E0 = ics[2]

## Kinetic Parameters
Vmax =  (k2*E0)
Km   =  ((k_1 + k2) / k1)

## Equations 
rxns = []
rxns.append(diff(S,t) == -Vmax*S/(Km+S))
rxns.append(diff(P,t) == Vmax*S/(Km+S))

## Calculation Parameters
## Same as above!

des = []
for i in range(0,len(rxns)):
    des.append(rxns[i].rhs())
sol = desolve_system_rk4(des,vars,[0,S0,P0],ivar,end_points=end_points,step=stepsize)

## Process the output into a form that can be graphed
sols=[]

for i in range(0,len(sol[0])-1):
    sols.append([])
for i in range(steps):
    for j in range(0,(len(sol[i])-1)):
        sols[j].append([sol[i][0],sol[i][j+1]])

################################################
####   Unnecessarily Fancy Plotting Stuff   ####
################################################
### 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)'

### USER-DEFINED PARAMETERS
#title = "Enzymatic Reaction Kinetics"                    ## Graph Title
legend_labels = ['[S]MM','[P]MM']                         ## Legend Labels
colors = ['salmon','lightgreen']                          ## Legend Colors
#x_label = 'time'                                         ## X-axis label
#y_label = 'concentration'                                ## Y-axis label

### Create a plot object
#a = plot([])
#a.axes_labels([x_label,y_label])
#a.set_axes_range(xmin=0,xmax=150,ymin=0,ymax=10)

### Add the desired lines to the plot
for i in range(0,len(sols)):
    a += list_plot(sols[i],color=colors[mod(i,len(colors))],legend_label=legend_labels[i],size='4')
    
### Set the plot parameters
#a += text(title,(a.xmax()/1.8,a.ymax()),color='black',fontsize=15)
a.axes_label_color('grey')
a.set_legend_options(ncol=3,borderaxespad=5,back_color='whitesmoke',fancybox=true)

show(a)

show('Vmax = '+ Vmax.str())
show('Km = '+ Km.str())
#plot?

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\verb|Vmax|\phantom{\verb!x!}\verb|=|\phantom{\verb!x!}\verb|0.00500000000000000|$$

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\verb|Km|\phantom{\verb!x!}\verb|=|\phantom{\verb!x!}\verb|3.00000000000000|$$