<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 OCNC2017 STEPS Tutorial Execise
In this notebook we will try to create a STEPS simulation script from scratch by modifying the examples given in the tutorial. Please follow the tutor's instruction step by step.
Build a biochemical model
Here is the example of a 2nd order reaction $A+B\overset{k}{\rightarrow}C$, where the reaction constant $k$ is set to 200 /uM.s
Step1:  For complex model, we can break it down into elementary reactions, for example, the following model
$E+S\underset{k_{-1}}{\overset{k_{1}}{\rightleftarrows}}ES\overset{k_{2}}{\rightarrow}E+P$
is broken down into 3 reactions in STEPS
1
Step2:  Setup geometry
You can easily setup the geometry of a well-mixed model by providing the volume of the geometry as well as the voume system it associated with.
Step3:  Create a random number generator
You can use the follow code to create a random number generator for the simulation, currently available generators are "mt19937" and "r123".
Step4:  Execise 2
Step5:  Create and initialize a solver
For well-mixed simulation, we create a "wmdirect" solver and initialize it by adding molecules to the compartment.
Step6:  Run the solver and gather simulation data
After that we can run the solver until it reaches a specific time point, say 0.1 second. You can gather simulation data such as molecule counts using many STEPS APIs, for example
Step7:  In practice, it is often necessary to store simulation data in a numpy array or a file for plotting or further analysis. For example, here we record the number of molcules using numpy array.
Step8:  Let's check what is inside the array now
Step9:  Execise 3
Step10:  Visuzalize simulation data
Visuzliation is often needed to analyze the behavior of the simulation, here we use Matplotlib to plot the data.
Step11:  Execise 4
Step12:  Here is the complete script for our well-mixed kinase simulation
Step13:  From well-mixed simulation to spatial simulation
To convert a well-mixed simulation to a spatial one, here are the basic steps
Step14:  We now import a tetrahedral mesh using the steps.utilities.meshio module to replace the well-mixed geometry
Step15:  Finally, we replace the "Wmdirect" solver with the spatial "Tetexact" solver
Step16:  The "Wmdirect" solver and the "Tetexact" solver share most of the APIs, so we can reuse our old script for simulation control and plotting
Step17:  Execise 5
Step18:  Here is the modified script
<ASSISTANT_TASK:>
Python Code:
# Import biochemical model module
import steps.model as smod
# Create model container
mdl = smod.Model()
# Create chemical species
A = smod.Spec('A', mdl)
B = smod.Spec('B', mdl)
C = smod.Spec('C', mdl)
# Create reaction set container
vsys = smod.Volsys('vsys', mdl)
# Create reaction
# A + B - > C with rate 200 /uM.s
reac_f = smod.Reac('reac_f', vsys, lhs=[A,B], rhs = [C])
reac_f.setKcst(200e6)
Explanation: OCNC2017 STEPS Tutorial Execise
In this notebook we will try to create a STEPS simulation script from scratch by modifying the examples given in the tutorial. Please follow the tutor's instruction step by step.
Build a biochemical model
Here is the example of a 2nd order reaction $A+B\overset{k}{\rightarrow}C$, where the reaction constant $k$ is set to 200 /uM.s:
End of explanation
# Import biochemical model module
import steps.model as smod
# Create model container
execise_mdl = smod.Model()
# Create chemical species
MEKp = smod.Spec('MEKp', execise_mdl)
ERK = smod.Spec('ERK', execise_mdl)
MEKpERK = smod.Spec('MEKpERK', execise_mdl)
ERKp = smod.Spec('ERKp', execise_mdl)
# Create reaction set container (volume system)
execise_vsys = smod.Volsys('execise_vsys', execise_mdl)
# Create reactions (Do it yourself)
# MEKp + ERK -> MEKpERK, rate constant 16.2*10e6
# MEKpERK -> MEKp + ERK, rate constant 0.6
# MEKpERK -> MEKp + ERKp, rate constant 0.15
Explanation: For complex model, we can break it down into elementary reactions, for example, the following model
$E+S\underset{k_{-1}}{\overset{k_{1}}{\rightleftarrows}}ES\overset{k_{2}}{\rightarrow}E+P$
is broken down into 3 reactions in STEPS
1: $E+S\overset{k_{1}}{\rightarrow}ES$
2: $ES\overset{k_{-1}}{\rightarrow}E+S$
3: $ES\overset{k_{2}}{\rightarrow}E+P$
Execise 1: Create a kinase reaction model in STEPS
Modify the script below for this kinase reaction system:
$MEKp+ERK\underset{0.6}{\overset{16.2*10^{6}}{\rightleftarrows}}MEKpERK\overset{0.15}{\rightarrow}MEKp+ERKp$
Hint: Break it down in to these elementary reactions
1: $MEKp+ERK\overset{16.2*10^{6}}{\rightarrow}MEKpERK$
2: $MEKpERK\overset{0.6}{\rightarrow}MEKp+ERK$
3: $MEKpERK\overset{0.15}{\rightarrow}MEKp+ERKp$
End of explanation
# Import geometry module
import steps.geom as sgeom
# Create well-mixed geometry container
wmgeom = sgeom.Geom()
# Create cytosol compartment
cyt = sgeom.Comp('cyt', wmgeom)
# Give volume to cyt (1um^3)
cyt.setVol(1.0e-18)
# Assign reaction set to compartment
cyt.addVolsys('vsys')
Explanation: Setup geometry
You can easily setup the geometry of a well-mixed model by providing the volume of the geometry as well as the voume system it associated with.
End of explanation
# Import random number generator module
import steps.rng as srng
# Create random number generator, with buffer size as 256
r = srng.create('mt19937', 256)
# Initialise with some seed
r.initialize(899)
# Could use time to get random seed
#import time
#r.initialize(int(time.time()))
Explanation: Create a random number generator
You can use the follow code to create a random number generator for the simulation, currently available generators are "mt19937" and "r123".
End of explanation
# Import biochemical model module
import steps.model as smod
# Create model container
execise_mdl = smod.Model()
# Create chemical species
MEKp = smod.Spec('MEKp', execise_mdl)
ERK = smod.Spec('ERK', execise_mdl)
MEKpERK = smod.Spec('MEKpERK', execise_mdl)
ERKp = smod.Spec('ERKp', execise_mdl)
# Create reaction set container (volume system)
execise_vsys = smod.Volsys('execise_vsys', execise_mdl)
# Create reactions (Do it yourself)
# MEKp + ERK -> MEKpERK, rate constant 16.2*10e6
MEKp_ERK_to_MEKpERK = smod.Reac('MEKp_ERK_to_MEKpERK', execise_vsys, lhs=[MEKp,ERK], rhs = [MEKpERK])
MEKp_ERK_to_MEKpERK.setKcst(16.2e6)
# MEKpERK -> MEKp + ERK, rate constant 0.6
MEKpERK_to_MEKp_ERK = smod.Reac('MEKpERK_to_MEKp_ERK', execise_vsys, lhs = [MEKpERK], rhs=[MEKp,ERK])
MEKpERK_to_MEKp_ERK.setKcst(0.6)
# MEKpERK -> MEKp + ERKp, rate constant 0.15
MEKpERK_to_MEKp_ERKp = smod.Reac('MEKpERK_to_MEKp_ERKp', execise_vsys, lhs = [MEKpERK], rhs=[MEKp,ERKp])
MEKpERK_to_MEKp_ERKp.setKcst(0.15)
####### You script after execise 1 should look like above #######
# Create a compartment of  0.1um^3
# Associate the compartment with the volume system 'vsys'
# Create and initialize a 'r123' random number generator
Explanation: Execise 2: Create the geometry and random number generator for the kinase reaction model
Let's continue our kinase reaction simulation script, here are the tasks:
1. Create a compartment of $0.1um^{3}$ (note that STEPS uses S.I units)
2. Associate the compartment with the volume system we've previously created
3. Create a "r123" random number generator and initialize it with with some seed
End of explanation
# Import solver module
import steps.solver as ssolv
# Create Well-mixed Direct solver
sim_direct = ssolv.Wmdirect(mdl, wmgeom, r)
# Inject 10 ‘A’ molecules
sim_direct.setCompCount('cyt','A', 10)
# Set concentration of ‘B’ molecules
sim_direct.setCompConc('cyt', 'B', 0.0332e-6)
Explanation: Create and initialize a solver
For well-mixed simulation, we create a "wmdirect" solver and initialize it by adding molecules to the compartment.
End of explanation
# Run simulation for 0.1s
sim_direct.run(0.1)
# Return the number of A molecules
sim_direct.getCompCount('cyt', 'A')
Explanation: Run the solver and gather simulation data
After that we can run the solver until it reaches a specific time point, say 0.1 second. You can gather simulation data such as molecule counts using many STEPS APIs, for example
End of explanation
# Reset the solver and reinitizlize molecule counts
sim_direct.reset()
# Inject 10 ‘A’ molecules
sim_direct.setCompCount('cyt','A', 10)
# Set concentration of ‘B’ molecules
sim_direct.setCompConc('cyt', 'B', 0.0332e-6)
# Import numpy
import numpy as np
# Create time-point numpy array, starting at time 0, end at 0.5 second and record data every 0.001 second
tpnt = np.arange(0.0, 0.501, 0.001)
# Calculate number of time points
n_tpnts = len(tpnt)
# Create data array, initialised with zeros
res_direct = np.zeros([n_tpnts, 3])
# Run simulation and record data
for t in range(0, n_tpnts):
    sim_direct.run(tpnt[t])
    res_direct[t,0] = sim_direct.getCompCount('cyt','A')
    res_direct[t,1] = sim_direct.getCompCount('cyt','B')
    res_direct[t,2] = sim_direct.getCompCount('cyt','C')
Explanation: In practice, it is often necessary to store simulation data in a numpy array or a file for plotting or further analysis. For example, here we record the number of molcules using numpy array.
End of explanation
print(res_direct)
Explanation: Let's check what is inside the array now:
End of explanation
# Import biochemical model module
import steps.model as smod
# Create model container
execise_mdl = smod.Model()
# Create chemical species
MEKp = smod.Spec('MEKp', execise_mdl)
ERK = smod.Spec('ERK', execise_mdl)
MEKpERK = smod.Spec('MEKpERK', execise_mdl)
ERKp = smod.Spec('ERKp', execise_mdl)
# Create reaction set container (volume system)
execise_vsys = smod.Volsys('execise_vsys', execise_mdl)
# Create reactions (Do it yourself)
# MEKp + ERK -> MEKpERK, rate constant 16.2*10e6
MEKp_ERK_to_MEKpERK = smod.Reac('MEKp_ERK_to_MEKpERK', execise_vsys, lhs=[MEKp,ERK], rhs = [MEKpERK])
MEKp_ERK_to_MEKpERK.setKcst(16.2e6)
# MEKpERK -> MEKp + ERK, rate constant 0.6
MEKpERK_to_MEKp_ERK = smod.Reac('MEKpERK_to_MEKp_ERK', execise_vsys, lhs = [MEKpERK], rhs=[MEKp,ERK])
MEKpERK_to_MEKp_ERK.setKcst(0.6)
# MEKpERK -> MEKp + ERKp, rate constant 0.15
MEKpERK_to_MEKp_ERKp = smod.Reac('MEKpERK_to_MEKp_ERKp', execise_vsys, lhs = [MEKpERK], rhs=[MEKp,ERKp])
MEKpERK_to_MEKp_ERKp.setKcst(0.15)
####### You script after execise 1 should look like above #######
# Create a compartment of  0.1um^3
import steps.geom as sgeom
execise_wmgeom = sgeom.Geom()
execise_cyt = sgeom.Comp('execise_cyt', execise_wmgeom)
execise_cyt.setVol(0.1e-18)
# Associate the compartment with the volume system 'vsys'
execise_cyt.addVolsys('execise_vsys')
# Create and initialize a 'r123' random number generator
import steps.rng as srng
execise_r = srng.create('r123', 256)
execise_r.initialize(1)
####### You script after execise 2 should look like above #######
# Create a "wmdirect" solver and set the initial condition:
# MEKp = 1uM
# ERK = 1.5uM
# Run the simulation for 30 seconds, record concerntrations of each molecule every 0.01 seconds.
Explanation: Execise 3: Run your kinase model in STEPS
Here are the tasks:
1. Create a "wmdirect" solver and set the initial condition:
    * MEKp = 1uM
    * ERK = 1.5uM
2. Run the simulation for 30 seconds, record concerntrations of each molecule every 0.01 seconds.
End of explanation
from pylab import *
%matplotlib inline
plot(tpnt, res_direct[:,0], label='A')
plot(tpnt, res_direct[:,1], label='B')
plot(tpnt, res_direct[:,2], label='C')
ylabel('Number of molecules')
xlabel('Time(sec)')
legend()
show()
Explanation: Visuzalize simulation data
Visuzliation is often needed to analyze the behavior of the simulation, here we use Matplotlib to plot the data.
End of explanation
# Import biochemical model module
import steps.model as smod
# Create model container
execise_mdl = smod.Model()
# Create chemical species
MEKp = smod.Spec('MEKp', execise_mdl)
ERK = smod.Spec('ERK', execise_mdl)
MEKpERK = smod.Spec('MEKpERK', execise_mdl)
ERKp = smod.Spec('ERKp', execise_mdl)
# Create reaction set container (volume system)
execise_vsys = smod.Volsys('execise_vsys', execise_mdl)
# Create reactions (Do it yourself)
# MEKp + ERK -> MEKpERK, rate constant 16.2*10e6
MEKp_ERK_to_MEKpERK = smod.Reac('MEKp_ERK_to_MEKpERK', execise_vsys, lhs=[MEKp,ERK], rhs = [MEKpERK])
MEKp_ERK_to_MEKpERK.setKcst(16.2e6)
# MEKpERK -> MEKp + ERK, rate constant 0.6
MEKpERK_to_MEKp_ERK = smod.Reac('MEKpERK_to_MEKp_ERK', execise_vsys, lhs = [MEKpERK], rhs=[MEKp,ERK])
MEKpERK_to_MEKp_ERK.setKcst(0.6)
# MEKpERK -> MEKp + ERKp, rate constant 0.15
MEKpERK_to_MEKp_ERKp = smod.Reac('MEKpERK_to_MEKp_ERKp', execise_vsys, lhs = [MEKpERK], rhs=[MEKp,ERKp])
MEKpERK_to_MEKp_ERKp.setKcst(0.15)
####### You script after execise 1 should look like above #######
# Create a compartment of  0.1um^3
import steps.geom as sgeom
execise_wmgeom = sgeom.Geom()
execise_cyt = sgeom.Comp('execise_cyt', execise_wmgeom)
execise_cyt.setVol(0.1e-18)
# Associate the compartment with the volume system 'vsys'
execise_cyt.addVolsys('execise_vsys')
# Create and initialize a 'r123' random number generator
import steps.rng as srng
execise_r = srng.create('r123', 256)
execise_r.initialize(143)
####### You script after execise 2 should look like above #######
# Create a "wmdirect" solver and set the initial condition:
# MEKp = 1uM
# ERK = 1.5uM
import steps.solver as ssolv
execise_sim = ssolv.Wmdirect(execise_mdl, execise_wmgeom, execise_r)
execise_sim.setCompConc('execise_cyt','MEKp', 1e-6)
execise_sim.setCompConc('execise_cyt','ERK', 1.5e-6)
# Run the simulation for 30 seconds, record concerntrations of each molecule every 0.01 seconds.
import numpy as np
execise_tpnts = np.arange(0.0, 30.01, 0.01)
n_tpnts = len(execise_tpnts)
execise_res = np.zeros([n_tpnts, 4])
# Run simulation and record data
for t in range(0, n_tpnts):
    execise_sim.run(execise_tpnts[t])
    execise_res[t,0] = execise_sim.getCompCount('execise_cyt','MEKp')
    execise_res[t,1] = execise_sim.getCompCount('execise_cyt','ERK')
    execise_res[t,2] = execise_sim.getCompCount('execise_cyt','MEKpERK')
    execise_res[t,3] = execise_sim.getCompCount('execise_cyt','ERKp')
####### You script after execise 3 should look like above #######
# Plot execise_res
Explanation: Execise 4: Plot the results of the kinase simulation
Let's now plot the result of our execise.
End of explanation
# Import biochemical model module
import steps.model as smod
# Create model container
execise_mdl = smod.Model()
# Create chemical species
MEKp = smod.Spec('MEKp', execise_mdl)
ERK = smod.Spec('ERK', execise_mdl)
MEKpERK = smod.Spec('MEKpERK', execise_mdl)
ERKp = smod.Spec('ERKp', execise_mdl)
# Create reaction set container (volume system)
execise_vsys = smod.Volsys('execise_vsys', execise_mdl)
# Create reactions (Do it yourself)
# MEKp + ERK -> MEKpERK, rate constant 16.2*10e6
MEKp_ERK_to_MEKpERK = smod.Reac('MEKp_ERK_to_MEKpERK', execise_vsys, lhs=[MEKp,ERK], rhs = [MEKpERK])
MEKp_ERK_to_MEKpERK.setKcst(16.2e6)
# MEKpERK -> MEKp + ERK, rate constant 0.6
MEKpERK_to_MEKp_ERK = smod.Reac('MEKpERK_to_MEKp_ERK', execise_vsys, lhs = [MEKpERK], rhs=[MEKp,ERK])
MEKpERK_to_MEKp_ERK.setKcst(0.6)
# MEKpERK -> MEKp + ERKp, rate constant 0.15
MEKpERK_to_MEKp_ERKp = smod.Reac('MEKpERK_to_MEKp_ERKp', execise_vsys, lhs = [MEKpERK], rhs=[MEKp,ERKp])
MEKpERK_to_MEKp_ERKp.setKcst(0.15)
####### You script after execise 1 should look like above #######
# Create a compartment of  0.1um^3
import steps.geom as sgeom
execise_wmgeom = sgeom.Geom()
execise_cyt = sgeom.Comp('execise_cyt', execise_wmgeom)
execise_cyt.setVol(0.1e-18)
# Associate the compartment with the volume system 'vsys'
execise_cyt.addVolsys('execise_vsys')
# Create and initialize a 'r123' random number generator
import steps.rng as srng
execise_r = srng.create('r123', 256)
execise_r.initialize(143)
####### You script after execise 2 should look like above #######
# Create a "wmdirect" solver and set the initial condition:
# MEKp = 1uM
# ERK = 1.5uM
import steps.solver as ssolv
execise_sim = ssolv.Wmdirect(execise_mdl, execise_wmgeom, execise_r)
execise_sim.setCompConc('execise_cyt','MEKp', 1e-6)
execise_sim.setCompConc('execise_cyt','ERK', 1.5e-6)
# Run the simulation for 30 seconds, record concerntrations of each molecule every 0.01 seconds.
import numpy as np
execise_tpnts = np.arange(0.0, 30.01, 0.01)
n_tpnts = len(execise_tpnts)
execise_res = np.zeros([n_tpnts, 4])
# Run simulation and record data
for t in range(0, n_tpnts):
    execise_sim.run(execise_tpnts[t])
    execise_res[t,0] = execise_sim.getCompCount('execise_cyt','MEKp')
    execise_res[t,1] = execise_sim.getCompCount('execise_cyt','ERK')
    execise_res[t,2] = execise_sim.getCompCount('execise_cyt','MEKpERK')
    execise_res[t,3] = execise_sim.getCompCount('execise_cyt','ERKp')
####### You script after execise 3 should look like above #######
# Plot execise_res
from pylab import *
plot(execise_tpnts, execise_res[:,0], label='MEKp')
plot(execise_tpnts, execise_res[:,1], label='ERK')
plot(execise_tpnts, execise_res[:,2], label='MEKpERK')
plot(execise_tpnts, execise_res[:,3], label='ERKp')
ylabel('Number of molecules')
xlabel('Time(sec)')
legend()
show()
####### You script after execise 4 should look like above #######
Explanation: Here is the complete script for our well-mixed kinase simulation:
End of explanation
# Import biochemical model module
import steps.model as smod
# Create model container
mdl = smod.Model()
# Create chemical species
A = smod.Spec('A', mdl)
B = smod.Spec('B', mdl)
C = smod.Spec('C', mdl)
# Create reaction set container
vsys = smod.Volsys('vsys', mdl)
# Create reaction
# A + B - > C with rate 200 /uM.s
reac_f = smod.Reac('reac_f', vsys, lhs=[A,B], rhs = [C])
reac_f.setKcst(200e6)
###### Above is the previous well-mixed biochemical model
# We add diffusion rules for species A, B and C
diff_a = smod.Diff('diff_a', vsys, A)
diff_a.setDcst(0.02e-9)
diff_b = smod.Diff('diff_b', vsys, B)
diff_b.setDcst(0.02e-9)
diff_c = smod.Diff('diff_c', vsys, C)
diff_c.setDcst(0.02e-9)
Explanation: From well-mixed simulation to spatial simulation
To convert a well-mixed simulation to a spatial one, here are the basic steps:
1. Add diffusion rules for every diffusive species in the biochemical model.
2. Change the well-mixed geoemtry to a tetrahedral mesh.
3. Change the solver to "Tetexact".
First, let's see how to add diffusion rules in our example well-mixed model:
End of explanation
'''
# Import geometry module
import steps.geom as sgeom
# Create well-mixed geometry container
wmgeom = sgeom.Geom()
# Create cytosol compartment
cyt = sgeom.Comp('cyt', wmgeom)
# Give volume to cyt (1um^3)
cyt.setVol(1.0e-18)
# Assign reaction set to compartment
cyt.addVolsys('vsys')
'''
##### above is the old well-mixed geometry ##########
import steps.geom as sgeom
import steps.utilities.meshio as meshio
# Import the mesh
mesh = meshio.importAbaqus('meshes/1x1x1_cube.inp', 1.0e-6)[0]
# Create mesh-based compartment
cyt = sgeom.TmComp('cyt', mesh, range(mesh.ntets))
# Add volume system to the compartment
cyt.addVolsys('vsys')
Explanation: We now import a tetrahedral mesh using the steps.utilities.meshio module to replace the well-mixed geometry:
End of explanation
# Import solver module
import steps.solver as ssolv
'''
# Create Well-mixed Direct solver
sim_direct = ssolv.Wmdirect(mdl, wmgeom, r)
'''
##### above is the old well-mixed Wmdirect solver ##########
# Create a spatial Tetexact solver
sim_tetexact = ssolv.Tetexact(mdl, mesh, r)
Explanation: Finally, we replace the "Wmdirect" solver with the spatial "Tetexact" solver:
End of explanation
# Inject 10 ‘A’ molecules
sim_tetexact.setCompCount('cyt','A', 10)
# Set concentration of ‘B’ molecules
sim_tetexact.setCompConc('cyt', 'B', 0.0332e-6)
# Import numpy
import numpy as np
# Create time-point numpy array, starting at time 0, end at 0.5 second and record data every 0.001 second
tpnt = np.arange(0.0, 0.501, 0.001)
# Calculate number of time points
n_tpnts = len(tpnt)
# Create data array, initialised with zeros
res_tetexact = np.zeros([n_tpnts, 3])
# Run simulation and record data
for t in range(0, n_tpnts):
    sim_tetexact.run(tpnt[t])
    res_tetexact[t,0] = sim_tetexact.getCompCount('cyt','A')
    res_tetexact[t,1] = sim_tetexact.getCompCount('cyt','B')
    res_tetexact[t,2] = sim_tetexact.getCompCount('cyt','C')
    
from pylab import *
plot(tpnt, res_tetexact[:,0], label='A')
plot(tpnt, res_tetexact[:,1], label='B')
plot(tpnt, res_tetexact[:,2], label='C')
ylabel('Number of molecules')
xlabel('Time(sec)')
legend()
show()
Explanation: The "Wmdirect" solver and the "Tetexact" solver share most of the APIs, so we can reuse our old script for simulation control and plotting:
End of explanation
# Import biochemical model module
import steps.model as smod
# Create model container
execise_mdl = smod.Model()
# Create chemical species
MEKp = smod.Spec('MEKp', execise_mdl)
ERK = smod.Spec('ERK', execise_mdl)
MEKpERK = smod.Spec('MEKpERK', execise_mdl)
ERKp = smod.Spec('ERKp', execise_mdl)
# Create reaction set container (volume system)
execise_vsys = smod.Volsys('execise_vsys', execise_mdl)
# Create reactions (Do it yourself)
# MEKp + ERK -> MEKpERK, rate constant 16.2*10e6
MEKp_ERK_to_MEKpERK = smod.Reac('MEKp_ERK_to_MEKpERK', execise_vsys, lhs=[MEKp,ERK], rhs = [MEKpERK])
MEKp_ERK_to_MEKpERK.setKcst(16.2e6)
# MEKpERK -> MEKp + ERK, rate constant 0.6
MEKpERK_to_MEKp_ERK = smod.Reac('MEKpERK_to_MEKp_ERK', execise_vsys, lhs = [MEKpERK], rhs=[MEKp,ERK])
MEKpERK_to_MEKp_ERK.setKcst(0.6)
# MEKpERK -> MEKp + ERKp, rate constant 0.15
MEKpERK_to_MEKp_ERKp = smod.Reac('MEKpERK_to_MEKp_ERKp', execise_vsys, lhs = [MEKpERK], rhs=[MEKp,ERKp])
MEKpERK_to_MEKp_ERKp.setKcst(0.15)
########### execise 5.1: Add diffusion constants
#    * MEKp = 30e-12 m^2/s
#    * ERK = 30e-12 m^2/s
#    * MEKpERK = 10e-12 m^2/s
####### You script after execise 1 should look like above #######
########### execise 5.2: Replace the geometry to use mesh 'meshes/sp_0.1v_1046.inp'
# Create a compartment of  0.1um^3
import steps.geom as sgeom
execise_wmgeom = sgeom.Geom()
execise_cyt = sgeom.Comp('execise_cyt', execise_wmgeom)
execise_cyt.setVol(0.1e-18)
# Associate the compartment with the volume system 'vsys'
execise_cyt.addVolsys('execise_vsys')
# Create and initialize a 'r123' random number generator
import steps.rng as srng
execise_r = srng.create('r123', 256)
execise_r.initialize(143)
####### You script after execise 2 should look like above #######
# Create a "wmdirect" solver and set the initial condition:
# MEKp = 1uM
# ERK = 1.5uM
import steps.solver as ssolv
########### execise 5.3: Change the solver to Tetexact
execise_sim = ssolv.Wmdirect(execise_mdl, execise_wmgeom, execise_r)
execise_sim.setCompConc('execise_cyt','MEKp', 1e-6)
execise_sim.setCompConc('execise_cyt','ERK', 1.5e-6)
# Run the simulation for 30 seconds, record concerntrations of each molecule every 0.01 seconds.
import numpy as np
execise_tpnts = np.arange(0.0, 30.01, 0.01)
n_tpnts = len(execise_tpnts)
execise_res = np.zeros([n_tpnts, 4])
# Run simulation and record data
for t in range(0, n_tpnts):
    execise_sim.run(execise_tpnts[t])
    execise_res[t,0] = execise_sim.getCompCount('execise_cyt','MEKp')
    execise_res[t,1] = execise_sim.getCompCount('execise_cyt','ERK')
    execise_res[t,2] = execise_sim.getCompCount('execise_cyt','MEKpERK')
    execise_res[t,3] = execise_sim.getCompCount('execise_cyt','ERKp')
####### You script after execise 3 should look like above #######
# Plot execise_res
from pylab import *
plot(execise_tpnts, execise_res[:,0], label='MEKp')
plot(execise_tpnts, execise_res[:,1], label='ERK')
plot(execise_tpnts, execise_res[:,2], label='MEKpERK')
plot(execise_tpnts, execise_res[:,3], label='ERKp')
ylabel('Number of molecules')
xlabel('Time(sec)')
legend()
show()
####### You script after execise 4 should look like above #######
Explanation: Execise 5: Modify your well-mixed kinase simulation to a spatial one
Let's now convert the below well-mixed kinase model to a spatial one, here are the tasks:
Add diffusion constants:
MEKp = 30e-12 $m^2/s$
ERK = 30e-12 $m^2/s$
MEKpERK = 10e-12 $m^2/s$
Replace the geometry to use mesh 'meshes/sp_0.1v_1046.inp'
Change the solver to Tetexact
Run the simulation again
End of explanation
# Import biochemical model module
import steps.model as smod
# Create model container
execise_mdl = smod.Model()
# Create chemical species
MEKp = smod.Spec('MEKp', execise_mdl)
ERK = smod.Spec('ERK', execise_mdl)
MEKpERK = smod.Spec('MEKpERK', execise_mdl)
ERKp = smod.Spec('ERKp', execise_mdl)
# Create reaction set container (volume system)
execise_vsys = smod.Volsys('execise_vsys', execise_mdl)
# Create reactions (Do it yourself)
# MEKp + ERK -> MEKpERK, rate constant 16.2*10e6
MEKp_ERK_to_MEKpERK = smod.Reac('MEKp_ERK_to_MEKpERK', execise_vsys, lhs=[MEKp,ERK], rhs = [MEKpERK])
MEKp_ERK_to_MEKpERK.setKcst(16.2e6)
# MEKpERK -> MEKp + ERK, rate constant 0.6
MEKpERK_to_MEKp_ERK = smod.Reac('MEKpERK_to_MEKp_ERK', execise_vsys, lhs = [MEKpERK], rhs=[MEKp,ERK])
MEKpERK_to_MEKp_ERK.setKcst(0.6)
# MEKpERK -> MEKp + ERKp, rate constant 0.15
MEKpERK_to_MEKp_ERKp = smod.Reac('MEKpERK_to_MEKp_ERKp', execise_vsys, lhs = [MEKpERK], rhs=[MEKp,ERKp])
MEKpERK_to_MEKp_ERKp.setKcst(0.15)
########### execise 5.1: Add diffusion constants
#    * MEKp = 30e-12 m^2/s
#    * ERK = 30e-12 m^2/s
#    * MEKpERK = 10e-12 m^2/s
diff_MEKp = smod.Diff('diff_MEKp', execise_vsys, MEKp)
diff_MEKp.setDcst(30e-12)
diff_ERK = smod.Diff('diff_ERK', execise_vsys, ERK)
diff_ERK.setDcst(30e-12)
diff_MEKpERK = smod.Diff('diff_MEKpERK', execise_vsys, MEKpERK)
diff_MEKpERK.setDcst(10e-12)
####### You script after execise 1 should look like above #######
########### execise 5.2: Replace the geometry to use mesh 'meshes/sp_0.1v_1046.inp'
import steps.geom as sgeom
import steps.utilities.meshio as meshio
mesh = meshio.importAbaqus('meshes/sp_0.1v_1046.inp', 1.0e-6)[0]
execise_cyt = sgeom.TmComp('execise_cyt', mesh, range(mesh.ntets))
execise_cyt.addVolsys('execise_vsys')
# Create and initialize a 'r123' random number generator
import steps.rng as srng
execise_r = srng.create('r123', 256)
execise_r.initialize(143)
####### You script after execise 2 should look like above #######
# Create a "wmdirect" solver and set the initial condition:
# MEKp = 1uM
# ERK = 1.5uM
import steps.solver as ssolv
########### execise 5.3: Change the solver to Tetexact
execise_sim = ssolv.Tetexact(execise_mdl, mesh, execise_r)
execise_sim.setCompConc('execise_cyt','MEKp', 1e-6)
execise_sim.setCompConc('execise_cyt','ERK', 1.5e-6)
# Run the simulation for 30 seconds, record concerntrations of each molecule every 0.01 seconds.
import numpy as np
execise_tpnts = np.arange(0.0, 30.01, 0.01)
n_tpnts = len(execise_tpnts)
execise_res = np.zeros([n_tpnts, 4])
# Run simulation and record data
for t in range(0, n_tpnts):
    execise_sim.run(execise_tpnts[t])
    execise_res[t,0] = execise_sim.getCompCount('execise_cyt','MEKp')
    execise_res[t,1] = execise_sim.getCompCount('execise_cyt','ERK')
    execise_res[t,2] = execise_sim.getCompCount('execise_cyt','MEKpERK')
    execise_res[t,3] = execise_sim.getCompCount('execise_cyt','ERKp')
####### You script after execise 3 should look like above #######
# Plot execise_res
from pylab import *
plot(execise_tpnts, execise_res[:,0], label='MEKp')
plot(execise_tpnts, execise_res[:,1], label='ERK')
plot(execise_tpnts, execise_res[:,2], label='MEKpERK')
plot(execise_tpnts, execise_res[:,3], label='ERKp')
ylabel('Number of molecules')
xlabel('Time(sec)')
legend()
show()
####### You script after execise 4 should look like above #######
Explanation: Here is the modified script
End of explanation
<END_TASK>