Software Tutorial: Traveling Down an Attractant Gradient
In the previous tutorial, we simulated the behavior of a bacterium moving up the concentration gradient. In this tutorial, we will simulate the opposite - when the bacterium is not in luck and moves down a concentration gradient.
To get started, open Visual Studio Code, and click File > Open Folder.... Open the EColiSimulations folder from the first tutorial. If you would rather not follow along below, you can download a completed BioNetGen file here: removal.bngl.
We also will build a Jupyter notebook in this tutorial for plotting the concentrations of different particles over time. To do so, you should save a copy of your plotter_up.ipynb file called plotter_down.ipynb; if you would rather not follow along, we provide a completed notebook here: plotter_down.ipynb
Before running this notebook, make sure the following dependencies are installed.
| Installation Link | Version | Check install/version |
|---|---|---|
| Python3 | 3.6+ | python --version |
| Jupyter Notebook | 4.4.0+ | jupyter --version |
| Numpy | 1.14.5+ | pip list \| grep numpy |
| Matplotlib | 3.0+ | pip list \| grep matplotlib |
| Colorspace or with pip | any | pip list \| grep colorspace |
Modeling a decreasing ligand gradient with a BioNetGen function
We have simulated how the concentration of phosphorylated CheY changes when the cell moves up the attractant gradient. The concentration dips, but over time, methylation states change so that they can compensate for the increased ligand-receptor binding and restore the equilibrium of phosphorylated CheY. What if instead ligands are removed, as we would see if the bacterium is traveling down an attractant gradient? We might imagine that we would see an increase in phosphorylated CheY to increase tumbling and change course, followed by a return to steady state. But is this what we will see?
To simulate the cell traveling down an attractant gradient, we will add a kill reaction removing unbound ligand at a constant rate. To do so, we will add the following rule within the reaction rules section.
#Simulate ligand removal
LigandGone: L(t) -> 0 k_gone
In the parameters section, we start by defining k_gone to be 0.3, so that d[L]/dt = -0.3[L]. The solution of this differential equation is [L] = 107e-0.3t. We will also change the initial ligand concentration (L0) to be 1e7. Thus, the concentration of ligand becomes so low that ligand-receptor binding reaches 0 within 50 seconds.
k_gone 0.3
L0 1e7
We will set the initial concentrations of all species to be the final steady state concentrations from the result for our adaptation.bngl model, and see if after reducing the concentration of unbound ligand gradually, the simulation can restore these concentrations to steady state.
First, visit the adaptation.bngl model and add the concentration for each combination of methylation state and ligand binding state of the receptor complex to the observables section. Then run this simulation with L0 equal to 1e7.
When the simulation is finished, visit the adaptation folder from the adaptation tutorial and find the simulation result at the final time point.
When the model finishes running, input the final concentrations of molecules to the species section of our removal.bngl model. Here is what we have.
begin species
@EC:L(t) L0
@PM:T(l!1,r,Meth~A,Phos~U).L(t!1) 1190
@PM:T(l!1,r,Meth~B,Phos~U).L(t!1) 2304
@PM:T(l!1,r,Meth~C,Phos~U).L(t!1) 2946
@PM:T(l!1,r,Meth~A,Phos~P).L(t!1) 2
@PM:T(l!1,r,Meth~B,Phos~P).L(t!1) 156
@PM:T(l!1,r,Meth~C,Phos~P).L(t!1) 402
@CP:CheY(Phos~U) CheY0*0.71
@CP:CheY(Phos~P) CheY0*0.29
@CP:CheZ() CheZ0
@CP:CheB(Phos~U) CheB0*0.62
@CP:CheB(Phos~P) CheB0*0.38
@CP:CheR(t) CheR0
end species
Running the BioNetGen model
We are now ready to run our BioNetGen model. To do so, first add the following after end model to run our simulation over 1800 seconds.
generate_network({overwrite=>1})
simulate({method=>"ssa", t_end=>1800, n_steps=>1800})
The following code contains our complete simulation, which can also be downloaded here: removal.bngl.
begin model
begin molecule types
L(t)
T(l,r,Meth~A~B~C,Phos~U~P)
CheY(Phos~U~P)
CheZ()
CheB(Phos~U~P)
CheR(t)
end molecule types
begin observables
Molecules bound_ligand L(t!1).T(l!1)
Molecules phosphorylated_CheY CheY(Phos~P)
Molecules low_methyl_receptor T(Meth~A)
Molecules medium_methyl_receptor T(Meth~B)
Molecules high_methyl_receptor T(Meth~C)
Molecules phosphorylated_CheB CheB(Phos~P)
end observables
begin parameters
NaV2 6.02e8 #Unit conversion to cellular concentration M/L -> #/um^3
miu 1e-6
L0 1e7
T0 7000
CheY0 20000
CheZ0 6000
CheR0 120
CheB0 250
k_lr_bind 8.8e6/NaV2 #ligand-receptor binding
k_lr_dis 35 #ligand-receptor dissociation
k_TaUnbound_phos 7.5 #receptor complex autophosphorylation
k_Y_phos 3.8e6/NaV2 #receptor complex phosphorylates Y
k_Y_dephos 8.6e5/NaV2 #Z dephosphorylates Y
k_TR_bind 2e7/NaV2 #Receptor-CheR binding
k_TR_dis 1 #Receptor-CheR dissociation
k_TaR_meth 0.08 #CheR methylates receptor complex
k_B_phos 1e5/NaV2 #CheB phosphorylation by receptor complex
k_B_dephos 0.17 #CheB autodephosphorylation
k_Tb_demeth 5e4/NaV2 #CheB demethylates receptor complex
k_Tc_demeth 2e4/NaV2 #CheB demethylates receptor complex
k_gone 0.3
end parameters
begin reaction rules
LigandReceptor: L(t) + T(l) <-> L(t!1).T(l!1) k_lr_bind, k_lr_dis
#Receptor complex (specifically CheA) autophosphorylation
#Rate dependent on methylation and binding states
#Also on free vs. bound with ligand
TaUnboundP: T(l,Meth~A,Phos~U) -> T(l,Meth~A,Phos~P) k_TaUnbound_phos
TbUnboundP: T(l,Meth~B,Phos~U) -> T(l,Meth~B,Phos~P) k_TaUnbound_phos*1.1
TcUnboundP: T(l,Meth~C,Phos~U) -> T(l,Meth~C,Phos~P) k_TaUnbound_phos*2.8
TaLigandP: L(t!1).T(l!1,Meth~A,Phos~U) -> L(t!1).T(l!1,Meth~A,Phos~P) 0
TbLigandP: L(t!1).T(l!1,Meth~B,Phos~U) -> L(t!1).T(l!1,Meth~B,Phos~P) k_TaUnbound_phos*0.8
TcLigandP: L(t!1).T(l!1,Meth~C,Phos~U) -> L(t!1).T(l!1,Meth~C,Phos~P) k_TaUnbound_phos*1.6
#CheY phosphorylation by T and dephosphorylation by CheZ
YPhos: T(Phos~P) + CheY(Phos~U) -> T(Phos~U) + CheY(Phos~P) k_Y_phos
YDephos: CheZ() + CheY(Phos~P) -> CheZ() + CheY(Phos~U) k_Y_dephos
#CheR binds to and methylates receptor complex
#Rate dependent on methylation states and ligand binding
TRBind: T(r) + CheR(t) <-> T(r!2).CheR(t!2) k_TR_bind, k_TR_dis
TaRUnboundMeth: T(r!2,l,Meth~A).CheR(t!2) -> T(r,l,Meth~B) + CheR(t) k_TaR_meth
TbRUnboundMeth: T(r!2,l,Meth~B).CheR(t!2) -> T(r,l,Meth~C) + CheR(t) k_TaR_meth*0.1
TaRLigandMeth: T(r!2,l!1,Meth~A).L(t!1).CheR(t!2) -> T(r,l!1,Meth~B).L(t!1) + CheR(t) k_TaR_meth*30
TbRLigandMeth: T(r!2,l!1,Meth~B).L(t!1).CheR(t!2) -> T(r,l!1,Meth~C).L(t!1) + CheR(t) k_TaR_meth*3
#CheB is phosphorylated by receptor complex, and autodephosphorylates
CheBphos: T(Phos~P) + CheB(Phos~U) -> T(Phos~U) + CheB(Phos~P) k_B_phos
CheBdephos: CheB(Phos~P) -> CheB(Phos~U) k_B_dephos
#CheB demethylates receptor complex
#Rate dependent on methyaltion states
TbDemeth: T(Meth~B) + CheB(Phos~P) -> T(Meth~A) + CheB(Phos~P) k_Tb_demeth
TcDemeth: T(Meth~C) + CheB(Phos~P) -> T(Meth~B) + CheB(Phos~P) k_Tc_demeth
#Simulate ligand removal
LigandGone: L(t) -> 0 k_gone
end reaction rules
begin compartments
EC 3 100 #um^3
PM 2 1 EC #um^2
CP 3 1 PM #um^3
end compartments
begin species
@EC:L(t) L0
@PM:T(l!1,r,Meth~A,Phos~U).L(t!1) 1190
@PM:T(l!1,r,Meth~B,Phos~U).L(t!1) 2304
@PM:T(l!1,r,Meth~C,Phos~U).L(t!1) 2946
@PM:T(l!1,r,Meth~A,Phos~P).L(t!1) 2
@PM:T(l!1,r,Meth~B,Phos~P).L(t!1) 156
@PM:T(l!1,r,Meth~C,Phos~P).L(t!1) 402
@CP:CheY(Phos~U) CheY0*0.71
@CP:CheY(Phos~P) CheY0*0.29
@CP:CheZ() CheZ0
@CP:CheB(Phos~U) CheB0*0.62
@CP:CheB(Phos~P) CheB0*0.38
@CP:CheR(t) CheR0
end species
end model
generate_network({overwrite=>1})
simulate({method=>"ssa", t_end=>1800, n_steps=>1800})
Now save your file and run the simulation by clicking Run BNG. The results will be saved in a new folder called removal/TIME contained in the current directory. Rename the folder from the timestamp to the value of k_gone, 0.3.
Open the newly created removal.gdat file and create a plot by clicking the Built-in plotting button.
What happens to the concentration of phosphorylated CheY? Are the concentrations of complexes at different methylation states restored to their levels before adding ligands to the adaptation.bngl model?
As we did in the tutorial simulating increasing ligand, we can try different values for k_gone. Change t_end in the simulate method to 1800 seconds, and run the simulation with k_gone equal to 0.01, 0.03, 0.05, 0.1, and 0.5.
Note: All simulation results are stored in the removal directory in your working directory. As you change the values of k_gone, rename the directory with the k_gone values instead of the timestamp for simplicity.
Visualizing the results of our simulation
We will use the jupyter notebook plotter_up.ipynb as a template for the plotter_down.ipynb file that we will use to visualize our results. First, we will specify the directories, model name, species of interest, and reaction rates. Put the Jupyter notebook in the same directory as removal.bngl or change the model_path accordingly.
model_path = "removal" #The folder containing the model
model_name = "removal" #Name of the model
target = "phosphorylated_CheY" #Target molecule
vals = [0.01, 0.03, 0.05, 0.1, 0.3, 0.5] #Gradients of interest
The second code block is the same as that provided in the previous tutorial. This code loads the simulation result at each time point from the .gdat file, which stores the concentration of all observables at all steps. It then plots the concentration of phosphorylated CheY over time.
import numpy as np
import sys
import os
import matplotlib.pyplot as plt
import colorspace
#Define the colors to use
colors = colorspace.qualitative_hcl(h=[0, 300.], c = 60, l = 70, pallete = "dynamic")(len(vals))
def load_data(val):
file_path = os.path.join(model_path, str(val), model_name + ".gdat")
with open(file_path) as f:
first_line = f.readline() #Read the first line
cols = first_line.split()[1:] #Get the col names (species names)
ind = 0
while cols[ind] != target:
ind += 1 #Get col number of target molecule
data = np.loadtxt(file_path) #Load the file
time = data[:, 0] #Time points
concentration = data[:, ind] #Concentrations
return time, concentration
def plot(val, time, concentration, ax, i):
legend = "k = - " + str(val)
ax.plot(time, concentration, label = legend, color = colors[i])
ax.legend()
return
fig, ax = plt.subplots(1, 1, figsize = (10, 8))
for i in range(len(vals)):
val = vals[i]
time, concentration = load_data(val)
plot(val, time, concentration, ax, i)
plt.xlabel("time")
plt.ylabel("concentration (#molecules)")
plt.title("Active CheY vs time")
ax.minorticks_on()
ax.grid(b = True, which = 'minor', axis = 'both', color = 'lightgrey', linewidth = 0.5, linestyle = ':')
ax.grid(b = True, which = 'major', axis = 'both', color = 'grey', linewidth = 0.8 , linestyle = ':')
plt.show()
Run the notebook. How does the value of k_gone impact the concentration of phosphorylated CheY? Why? Are the tumbling frequencies restored to the background frequency? As we return to the main text, we will show the resulting plots and discuss these questions.