How to Write a Process¶
Why Write Processes?¶
Vivarium comes with a number of processes ready for you to use, but combining processes to form composites will only take you so far. For many models, the existing processes will be insufficient, so you will need to create your own.
Note
Processes are the building blocks of Vivarium models, so creating a process to model a phenomenon you know well lets other modelers build on your expertise. Processes are a great way to share your knowledge with the modeling community!
Let’s Write a Process!¶
Suppose we want to model how hexokinase helps maintain intracellular glucose concentrations. Hexokinase catalyzes glucose phosphorylation, which we can describe by the following reaction:
In the reaction above and throughout this tutorial, we will use the following abbreviations:
GLC: D-Glucose
ATP: Adenosine triphosphate
G6P: α-D-Glucose-6-phosphate
ADP: Adenosine diphosphate
v: Rate of the forward reaction
HK: Hexokinase
Tip
Once you have worked through the example in this tutorial, try modeling a different reaction that you work with!
Conceptualize the Model¶
Make Assumptions¶
A critical component of modeling is deciding what parts to describe mechanistically and what to approximate. For our scenario, we could model this reaction using molecular dynamics, but this would be computationally intensive and would not scale to larger simulations. Instead, we will assume Michaelis-Menten, sequential bisubstrate kinetics:
Let’s also assume that diffusion is much faster than the reaction so that the concentrations of our enzyme and substrates are constant throughout the cell.
Note
If you actually wanted to model a reaction like this, the convenience kinetics process can be configured to model any Michaelis-Menten enzyme kinetics.
Translate Model into Updates¶
Processes in Vivarium work by repeatedly changing the state of the model. A process makes these changes by computing an update based on the model’s current state and a timestep \(t\).
For the current example, each update will include the following:
A decrease in GLC: \(-v t\)
A decrease in ATP: \(-v t\)
An increase in G6P: \(v t\)
An increase in ADP: \(v t\)
Note
In this example and most of the time in Vivarium, we work in terms of concentrations. We also normally use units of mM, but you can use different units so long as you are consistent.
Determine the Ports¶
We partition the overall simulation state into stores, which can be shared among processes. Each process declares ports, each of which will receive a store. When creating a process, you need to decide what ports to declare.
When someone else uses your process, they will create a composite of it and other processes. These processes will interact by sharing stores. While any number of your process’s ports may be linked to the same store, a port cannot be split between stores. This means that you should put in separate ports any variables that a user might want in separate stores.
For example, ATP and ADP are turned over rapidly in the cell, so a user might want to isolate those variables from others that get updated more slowly. We will therefore create two ports:
nucleoside_phosphates
: This port will store theATP
andADP
variables.cytoplasm
: This port will store theGLC
,G6P
, andHK
variables.
Implement the Model¶
To implement the model, create a new Python file named
glucose_phosphorylation.py
in the vivarium/processes/
directory.
Then we create a new class that inherits from
vivarium.core.process.Process
:
from vivarium.core.process import Process
class GlucosePhosphorylation(Process):
pass
The Constructor¶
In the constructor we can configure the process. We accept
configurations as a dictionary called initial_parameters
. For
example, we can let the user configure the kinetic parameters
\(k_{cat}\), \(K_{GLC}\), and \(K_{ATP}\). We can also
provide default values for these parameters.
The configurations (with any missing parameters filled in with defaults) and ports are passed to the superclass constructor to instantiate the process.
import copy
from vivarium.core.process import Process
class GlucosePhosphorylation(Process):
defaults = {
'k_cat': 2e-3,
'K_ATP': 5e-2,
'K_GLC': 4e-2,
}
def __init__(self, initial_parameters=None):
if initial_parameters == None:
initial_parameters = {}
parameters = copy.deepcopy(GlucosePhosphorylation.defaults)
parameters.update(initial_parameters)
super().__init__(parameters)
It turns out that the vivarium.core.process.Process
constructor handles merging initial_parameters
with the defaults
automatically, so we don’t actually have to include a constructor at
all!
Even though we’re just getting started on our process, let’s try it out!
At the bottom of the glucose_phosphorylation.py
file, instantiate
the process and take a look at some of its attributes:
if __name__ == '__main__':
parameters = {
'k_cat': 1.5,
}
my_process = GlucosePhosphorylation(parameters)
print(my_process.parameters['k_cat'])
print(my_process.parameters['K_ATP'])
Then run your code by executing the whole file:
$ python glucose_phosphorylation.py
1.5
0.05
Notice that the k_cat
parameter updated to the value we supplied and
that k_ATP
took on the default value.
Wait! Where did the parameters
attribute come from? We never
created that attribute, but
vivarium.core.process.Process
made it from the
parameters
argument we passed to its constructor. We’ll take
advantage of this in the next step.
Generating Updates¶
Now we can write the next_update
method, which generates updates for
each port based on a provided simulation state and timestep.
Warning
The states
parameter passed into the update function is
a view of the overall state, so it must not be changed.
def next_update(self, timestep, states):
# Get concentrations from state
cytoplasm = states['cytoplasm']
nucleoside_phosphates = states['nucleoside_phosphates']
hk = cytoplasm['HK']
glc = cytoplasm['GLC']
atp = nucleoside_phosphates['ATP']
# Get kinetic parameters
k_cat = self.parameters['k_cat']
k_atp = self.parameters['K_ATP']
k_glc = self.parameters['K_GLC']
# Compute reaction rate with michaelis-menten equation
rate = k_cat * hk * glc * atp / (
k_glc * k_atp + k_glc * atp + k_atp * glc + glc * atp)
# Compute concentration changes from rate and timestep
delta_glc = -rate * timestep
delta_atp = -rate * timestep
delta_g6p = rate * timestep
delta_adp = rate * timestep
# Compile changes into an update
update = {
'cytoplasm': {
'GLC': delta_glc,
'G6P': delta_g6p,
# We exclude HK because it doesn't change
},
'nucleoside_phosphates': {
'ATP': delta_atp,
'ADP': delta_adp,
},
}
return update
Now let’s test this update function by seeing how it changes a state we provide. Replace the testing code we added to the bottom of the file with this:
if __name__ == '__main__':
parameters = {
'k_cat': 1.5,
}
my_process = GlucosePhosphorylation(parameters)
state = {
'cytoplasm': {
'GLC': 1.0,
'G6P': 0.0,
'HK': 0.1,
},
'nucleoside_phosphates': {
'ATP': 2.0,
'ADP': 0.0,
},
}
update = my_process.next_update(3.0, state)
print(update['cytoplasm']['G6P'])
With these parameters, we can calculate the reaction rate:
Therefore, we expect the change in concentration of G6P to be:
Let’s see if our process models this reaction as we expect:
$ python glucose_phosphorylation.py
0.4221388367729832
Hooray! This is what we expected.
Ports Schema¶
Our process works, but we had to manually set the state. We also haven’t
shown yet how to apply the update we generate to the simulation state.
Luckily for us, these steps will be handled automatically by Vivarium.
We just need to create a ports_schema
method that provides a
schema. A schema is a nested dictionary that describes each
variable the process will interact with. Each variable is defined by a
dictionary of schema keys that specify its default value, how it
should be updated, and other properties.
For this example, our updates are expressed as deltas that should be added to the old value of the variable. This is the default, so the schema can leave out the updater specification. Still, we’ll specify one of the updaters for demonstration.
def ports_schema(self):
return {
'cytoplasm': {
'GLC': {
# accumulate means to add the updates
'_updater': 'accumulate',
'_default': 1.0,
'_properties': {
'mw': 1.0 * units.g / units.mol,
},
'_emit': True,
},
# accumulate is the default, so we don't need to specify
# updaters for the rest of the variables
'G6P': {
'_default': 0.0,
'_properties': {
'mw': 1.0 * units.g / units.mol,
},
'_emit': True,
},
'HK': {
'_default': 0.1,
'_properties': {
'mw': 1.0 * units.g / units.mol,
},
},
},
'nucleoside_phosphates': {
'ATP': {
'_default': 2.0,
'_emit': True,
},
'ADP': {
'_default': 0.0,
'_emit': True,
}
},
}
By convention, we usually put the ports_schema
method before the
next_update
method.
Now, we can run a simulation using Vivarium’s Engine
vivarium.core.engine.Engine()
like this:
from vivarium.core.engine import Engine
from vivarium.plots.simulation_output import plot_simulation_output
...
if __name__ == '__main__':
parameters = {
'k_cat': 1.5,
}
my_process = GlucosePhosphorylation(parameters)
# make a composite
my_composite = my_process.generate()
# run a simulation
sim = Engine(composite=my_composite)
total_time = 10
sim.update(total_time)
# get the timeseries
timeseries = sim.emitter.get_timeseries()
plot_simulation_output(timeseries, {}, './')
We use
vivarium.plots.simulation_output.plot_simulation_output
to
plot the output from our simulation. In simulation.png
you should
see an output plot like this:
Tip
If a process is erroneously reporting negative values, try decreasing the timestep.
Oops, it looks like the cytoplasmic GLC concentration dropped below zero
around time 8! This happens when the timestep is too long and so our
approximation doesn’t adjust fast enough to dropping concentrations. To
fix this, let’s change the timestep to 0.1
.
Note
You may be wondering, “What unit is the timestep in?” The answer is that it doesn’t matter! We just need the parameters and timestep to use the same unit of time.
Here’s the parameters
dictionary with the updated timestep:
parameters = {
'k_cat': 1.5,
'time_step': 0.1,
}
Now if we run the file again, we should get a simulation.png
like
this:
You can download the completed process file here
.
Great job; you’ve written a new process! Now consider writing one to model a mechanism you are familiar with.