18. Excitability enables probabilistic, transient differentiation
Design principles
Fast positive feedback combined with slow negative feedback can combine to give excitable circuits.
Excitability enables probabilistic, transient differentiation
Concepts
Excitable circuits have a unique, globally attracting rest state and a large stimulus can send the system on a long excursion before it returns to the resting state.
[30]:
# Colab setup 
import os, sys, subprocess
if "google.colab" in sys.modules:
cmd = "pip install upgrade biocircuits iqplot numba multiprocess watermark"
process = subprocess.Popen(cmd.split(), stdout=subprocess.PIPE, stderr=subprocess.PIPE)
stdout, stderr = process.communicate()
# 
import collections
import tqdm
import pandas as pd
import numpy as np
import scipy.integrate
import numba
import biocircuits
import bokeh.io
import bokeh.plotting
import iqplot
# Toggle for fully interactive plots
fully_interactive_plots = True
notebook_url = "localhost:8888"
bokeh.io.output_notebook()
Altered states
Cells can exist in a broad repertoire of physiologically distinct states with different properties. In the previous chapter we saw how spontaneous transitions to slow growing states can give rise to antibiotic persistence in bacteria, effectively sacrificing the growth of some cells to increase the survival probability of their clone. In fact, cells are constantly switching in and out of a myriad of different physiological states.
We have now encountered several different features of dynamical systems that each enable a distinct cellular capability:
Dynamical system property 
Cellular capability 

Monostable fixed point 
Remain in a constant state (homeostasis) 
Bistability / multistability 
Exist in two or more epigenetically stable states (e.g., toggle switch, MultiFate) 
Limit cycles 
Periodic oscillations (e.g., repressilator) 
In this chapter, we will add a new line to this table, showing how excitability takes advantage of noise to probabilistically initiate transient episodes of differentiation. Furthermore, we will see how the circuit architecture makes it possible for bacteria to independently modulate two key properties of these episodes: their frequency and their duration. We will discuss excitability in the context of cellular competence.
Excitability is tricky to define. Borrowing from Strogatz (in Nonlinear Dynamics and Chaos) an excitable system has two key properties:
It has a unique, globally attracting rest state;
A large enough stimulus can send the system on a long excursion before it returns to the resting state.
Excitability is prevalent
A forest is a classic example of an excitable system. Its resting state consists of many trees peacefully photosynthesizing and providing shade for hikers. Occasionally, though, a relatively small perturbation, perhaps an errant lightning strike or a careless camper, can dramatically alter the landscape, initiating a transient excursion wherein the forest burns, recovers, and then grows back again.
The most familiar example of excitability may be a toilet. Normally the toilet sits at rest in a more or less constant state. Pressing the handle strongly enough initiates a stereotyped sequence of events: the bowl drains, refills with water, and returns to its initial resting state. This is a deterministic system. But imagine adding noise by jiggling the handle. Most of the jiggles have no effect. But occasionally one can be strong enough to cross the threshold and initiate a flush. When that happens, the resulting flush cycle is the same as one initiated by a deliberate push. It does not matter how the threshold is crossed. In an excitable system like this, any thresholdcrossing perturbation—whether external or from noise—generates a similar dynamic trajectory.
Closer to molecular biology, the membrane potentials of neurons also represent excitable systems. When resting, they sit at their postsynaptic potential (PSP) and a stimulus results in them relaxing immediately back to their PSP. However, when they are excitable, a stimulus can lead to a much larger excursion, called an action potential, or spike, in which the membrane potential grows far beyond the PSP, before returning back to the PSP.
(Image adapted from Izhikevich, Dynamical Systems in Neuroscience, 2007.)
Cellular competence
Now let’s see how bacteria take advantage of excitability. When stressed, some bacteria, such as Bacillus subtilis, can activate a variety of stress response pathways. This activation can be heterogeneous, with different cells activating different responses, even in the same conditions.
One of these responses is to enter a competent state in which the cell can take up environmental DNA and recombine it into its own genome. This process represents a type of genetic exchange that could provide evolutionary advantages similar to sexual reproduction. But competence also carries risks. The cell might take up viral DNA, for instance. Hence, from a bet hedging point of view, it may be advantageous to activate it in some, but not all, cells within the population.
Some of the heterogeneity of stress response can be observed in the following timelapse movie. In the movie below, cells have been engineered to light up as pink/purple when they become competent. But in the same conditions, cells can also transform into dormant spores, which appear as white ovals in the microscope. So at least two stress responses are clearly evident in the same conditions:
From this movie, it is evident that competence is both probabilistic (only a fraction of cells become competent) and transient (competence is a temporary condition). These are both hallmarks of excitability.
We will now examine the genetic circuit responsible for this excitability and seek to address the following questions:
How do cells enter and exit competence?
How do cells regulate the frequency and duration of competence?
Can these properties be tuned?
Is noise necessary to induce competence?
The B. subtilis competence circuit
The genetic circuit responsible for competence is connected to receptors responsible for sensing the environment and other circuits within the cell. It can appear quite complicated, as shown in the figure below.
Image taken fromHamoan, et al., Microbiology, 2003
To follow the dynamics of the circuit, Süel, et al. inserted combinations of yellow (YFP) and cyan (CFP) fluorescent reporter genes for various competence genes such as comK, comG and comS in B. subtilis, enabling them to monitor their dynamics in single cells. ComK (blue channel) and ComG (red channel) activate together in competent cells.
Expression of ComK, which has many inputs, largely follows expression of ComG, whose primary input is ComK. (This is why the competent cells appeared purple (blue + red) in the movie above.) This suggested that ComK is the primary regulator of its own expression, at least under these conditions. A second key discovery came from observing ComS and ComG act the same time in two colors. In the movie below, you can see how ComS, in green, declines just as the cell becomes competent (turns on ComK, as shown by the red reporter). This provided compelling evidence for the negative regulatory feedback loop above.
Through these and other experiments, Süel, et al were able to identify a simpler 3gene core circuit, comprising ComK, ComS, and MecA, that could account for the probabilistic and transient nature of competence differentiation.
Here, ComK activates its own transcription in a positive feedback loop (orange) and indirectly represses ComS (dashed line). Both ComK and ComS are degraded by the MecA–ClpP–ClpC complex (which we will call MecA for short). Competition for degradation by MecA means that high levels of ComS can effectively inhibit degradation of ComK by MecA. In this way, ComS effectively stabilizes ComK. ComG is simply a target gene directly activated by ComK, whose promoter can be used to read out ComK activity.
How does positive autoregulation by ComK and negative feedback through ComS lead to excitable dynamics? Imagine a stochastic fluctuation in ComK expression that produces a sudden, transient increase in ComK levels. The positive feedback loop (orange arrows) amplifies the fluctuation through increased production of ComK. Initially, high basal levels of ComS prevent MecA from efficiently degrading ComK. On a longer time scale, however, the increased ComK inhibits ComS production, causing ComS levels to decline, which in turn allows MecA to begin degrading ComK instead of ComS. This brings ComK back down to its preexcitedment level. In this sense, ComS serves to help initiate competence by reducing degradation of ComK, but also serves to bring the exit of competence when it is repressed by the excessive ComK.
In what follows, we will explore this competence circuit, and discover a new circuit design principle: Fast positive feedback combined with slow negative feedback can combine to give noiseexcitable dynamics. Those excitable dynamics, in turn, enable transient, probabilistic differentiation.
Modeling the ComK circuit
To model the ComK circuit, we will follow Süel, et al. (Science, 2007). We start with the MecA complex, which degrades both ComK and ComS. ComS can inhibit degradation of ComK through competition, tying up the MecA complex so much that little is left to degrade ComK. More specifically, there are two competing reactions: degradation of ComK and degradation of ComS. We represent both using MichaelisMenten kinetics, giving the reaction schemes
\begin{align} \require{mhchem} &\ce{MecA + ComK <=>[\gamma_{+a}][\gamma_{a}] MecAComK >[\gamma_1] MecA} \\[1em] &\ce{MecA + ComS <=>[\gamma_{+b}][\gamma_{b}] MecAComS >[\gamma_2] MecA}. \end{align}
ComK also indirectly inhibits expression of ComS, as denoted by the dashed arrow. Although it is indirect, and could be more complex, for simplicity we represent this regulatory interaction with a Hill function. Finally, ComK has a basal expression level, which we will denote \(\alpha_k\). With these assumptions, we can write the dynamical equations for ComK, \(k\), and the MecAComK complex, \(m_K\), as
\begin{align} \frac{\mathrm{d}k}{\mathrm{d}t}&= \alpha_k + \beta_k\,\frac{(k/k_k)^n}{1 + (k/k_k)^n}  \gamma_a m_f k + \gamma_{a}m_K  \gamma_k k,\\[1em] \frac{\mathrm{d}m_K}{\mathrm{d}t} &= (\gamma_{a} + \gamma_1)m_K + \gamma_a m_f K, \end{align}
where we have denoted the concentration of MecAComK as \(m_K\) and of free MecA as \(m_f\). We have also assumed that ComK undergoes noncatalyzed degradation, with rate constant \(\gamma_k\).
With our assumption of MichaelisMenten kinetics on ComS degradation by MecA, we can now also write the dynamical equations for ComS, \(s\), and the ComSMecA complex, \(m_S\).
\begin{align} \frac{\mathrm{d}s}{\mathrm{d}t}&= \frac{\beta_s}{1 + (k/k_s)^p}  \gamma_b m_f s + \gamma_{b}m_S  \gamma_s s,\\[1em] \frac{\mathrm{d}m_S}{\mathrm{d}t} &= (\gamma_{b} + \gamma_2)m_S + \gamma_b m_f s. \end{align}
If the conversion of the MecACom complexes are fast compared to the transcriptional dynamics (\(\gamma_1\) and \(\gamma_2\) are large), we can make the quasisteady state approximation that \(\mathrm{d}m_K/\mathrm{d}t \approx \mathrm{d}m_S/\mathrm{d}t \approx 0\). Additionally, assuming the total amount of MecA complex is conserved, i.e., \(m_f + m_K + m_S = m_\mathrm{total} \approx \text{constant}\), we can simply the equations:
\begin{align} \frac{\mathrm{d}k}{\mathrm{d}t}&= \alpha_k + \beta_k\,\frac{(k/k_k)^n}{1 + (k/k_k)^n} \gamma_1 m_K  \gamma_k k,\\[1em] \frac{\mathrm{d}s}{\mathrm{d}t}&= \alpha_s + \frac{\beta_s}{1 + (k/k_s)^p}  \gamma_2 m_S  \gamma_s s. \end{align}
We can get expressions for \(m_K\) and \(m_S\) in terms of \(m_\mathrm{total}\), \(k\), and \(s\) by using the relations following our two assumptions,
\begin{align} (\gamma_{a} + \gamma_1)m_K + \gamma_a(m_\mathrm{total}  m_K  m_S) k &= 0,\\[1em] (\gamma_{b} + \gamma_2)m_S + \gamma_b(m_\mathrm{total}  m_K  m_S) s &= 0. \end{align}
We can write these in a more convenient form to eliminate \(m_K\) and \(m_S\).
\begin{align} (\gamma_{a} + \gamma_1 + \gamma_a k)m_K + \gamma_a k m_S &= \gamma_a m_\mathrm{total} k \\[1em] \gamma_b s m_K + (\gamma_{b} + \gamma_2 + \gamma_b s)m_S &= \gamma_b m_\mathrm{total} s. \end{align}
Solving this system of linear equations for \(m_K\) and \(m_S\) yields
\begin{align} m_K &= m_\mathrm{total}\,\frac{k/\Gamma_k}{1 + k/\Gamma_k + s/\Gamma_s}, \\[1em] m_S &= m_\mathrm{total}\,\frac{s/\Gamma_s}{1 + k/\Gamma_k + s/\Gamma_s}, \end{align}
where \(\Gamma_k \equiv (\gamma_{a} + \gamma_1)/\gamma_a\) and \(\Gamma_k \equiv (\gamma_{b} + \gamma_2)/\gamma_b\). Inserting these expressions back into our ODEs for the \(k\) and \(s\) dynamics and defining \(\delta_k\equiv \gamma_1 m_\mathrm{total}/\Gamma_k\) and \(\delta_s \equiv\gamma_2 m_\mathrm{total}/\Gamma_s\) for notational convenience yields
\begin{align} \frac{\mathrm{d}k}{\mathrm{d}t}&= \alpha_k + \beta_k\,\frac{(k/k_k)^n}{1 + (k/k_k)^n}  \frac{\delta_k k}{1 + k/\Gamma_k + k/\Gamma_s}  \gamma_k k\\[1em] \frac{\mathrm{d}S}{\mathrm{d}t}&= \alpha_s + \frac{\beta_s}{1 + (k/k_s)^p}  \frac{\delta_s s}{1 + k/\Gamma_k + s/\Gamma_s}  \gamma_s s. \end{align}
We can further nondimensionalize time, using the rate of ComK degradation by the MecA complex, \(\delta_k\). The parameters \(\Gamma_k\) and \(\Gamma_s\) are natural choices for nondimensionalizing \(K\) and \(S\).
After all this, we arrive at a set of just two dimensionless equations to approximately describe the system:
\begin{align} \frac{\mathrm{d}k}{\mathrm{d}t}&= a_k + \frac{b_k\,(k/\kappa_k)^n}{1 + (k/\kappa_k)^n}  \frac{k}{1 + k + s}  \gamma_k k, \\[1em] \frac{\mathrm{d}s}{\mathrm{d}t}&= a_s + \frac{b_s}{1 + (k/\kappa_s)^p}  \frac{s}{1 + k + s}  \gamma_s s. \end{align}
It is left to the reader to work out the relationships between the dimensionless and dimensional parameters.
Computational analysis of the competence dynamical system
With the differential equations finally in hand, we can now take this dynamical system out for a spin and see what it can do. We start by defining functions that compute the time derivatives in the dimensionless ComK and ComS equations immediately above.
[61]:
def dk_dt(k, s, p):
"""p is a CompetenceParams instance"""
return (
p.a_k
+ p.b_k * k ** p.n / (p.kappa_k ** p.n + k ** p.n)
 k / (1 + k + s)
 p.gamma_k * k
)
def ds_dt(k, s, p):
return (
p.a_s
+ p.b_s / (1 + (k / p.kappa_s) ** p.p)
 s / (1 + k + s)
 p.gamma_s * s
)
Identifying nullclines and fixed points
As we saw in previous chapters, plotting the nullclines, defined by setting each of the time derivatives to zero, can identify fixed points and provide qualitative insights into the dynamics possible in the system. Here, the nullclines are the solutions to \(\mathrm{d}k/\mathrm{d}t=0\) and \(\mathrm{d}s/\mathrm{d}t=0\). They are most easily represented in the form \(s(k)\), that is, \(s\) as a function of \(k\). First, consider the \(\mathrm{d}k/\mathrm{d}t\) nullcline:
\begin{align} &\frac{\mathrm{d}k}{\mathrm{d}t} = 0 \;\Rightarrow s\; = \frac{k}{a_k  \gamma_k k + b_k k^n / (\kappa_k^n + k^n)}  k  1. \end{align}
The nullcline given by \(\mathrm{d}s/\mathrm{d}t = 0\) is the positive root of the quadratic polynomial equation \(as^2 + b s + c = 0\) with
\begin{align} a &= \gamma_s,\\[1em] b &= 1 + \gamma_s(1+k)  A,\\[1em] c &= A(1+k), \end{align}
with
\begin{align} A = a_s + \frac{b_s}{(1 + (k/\kappa_s)^p}. \end{align}
We can code these up as functions:
[57]:
def k_nullcline(k, p):
abk = (
p.a_k
+ p.b_k * k ** p.n / (p.kappa_k ** p.n + k ** p.n)
 p.gamma_k * k
)
return k / abk  k  1
def s_nullcline(k, p):
A = p.a_s + p.b_s / (1 + (k / p.kappa_s) ** p.p)
a = p.gamma_s
b = 1  A + p.gamma_s * (1 + k)
c = A * (1 + k)
return 2 * c / (b + np.sqrt(b ** 2  4 * a * c))
We will begin our graphical analysis by plotting the nullclines, as we have done in past chapters. For flexibility, we will write a function that can plot the nullclines on either a logarithmic scale or a linear scale over a specified range of \(s\) and \(k\) values:
[59]:
def plot_nullclines(k_range, s_range, params, log=True, p=None, **kwargs):
"""
Compute and plot the nullclines.
"""
if log:
x_range = np.log10(k_range)
y_range = np.log10(s_range)
else:
x_range = k_range
y_range = s_range
if p is None:
p = bokeh.plotting.figure(
frame_width=260,
frame_height=260,
x_axis_label="log₁₀ k",
y_axis_label="log₁₀ s",
x_range=x_range,
y_range=y_range,
**kwargs
)
def filter_out_of_range(k, s_nck, s_ncs):
k_k = k.copy()
k_s = k.copy()
inds_k = np.logical_or(s_nck < s_range[0], s_nck > s_range[1])
inds_s = np.logical_or(s_ncs < s_range[0], s_ncs > s_range[1])
s_nck[inds_k] = np.nan
s_ncs[inds_s] = np.nan
k_k[inds_k] = np.nan
k_s[inds_s] = np.nan
return k_k, k_s, s_nck, s_ncs
if log:
k = np.logspace(np.log10(k_range[0]), np.log10(k_range[1]), 200)
s_nck = k_nullcline(k, params)
s_ncs = s_nullcline(k, params)
k_k, k_s, s_nck, s_ncs = filter_out_of_range(k, s_nck, s_ncs)
p.line(np.log10(k_k), np.log10(s_nck), color="#ff7f0e", line_width=3)
p.line(np.log10(k_s), np.log10(s_ncs), color="#9467bd", line_width=3)
else:
k = np.linspace(k_range[0], k_range[1], 200)
s_nck = k_nullcline(k, params)
s_ncs = s_nullcline(k, params)
k_k, k_s, s_nck, s_ncs = filter_out_of_range(k, s_nck, s_ncs)
p.line(k_k, s_nck, color="#ff7f0e")
p.line(k_s, s_ncs, color="#9467bd")
return p
We can now use this function to begin constructing the phase portrait using the parameter values from the Süel, et al. paper and make the plot with \(k\) and \(s\) on a logarithmic scale. Note that we are labeling our axes as “log₁₀ k” instead of choosing x_axis_type = 'log'
as a keyword argument when building the figure. The reason for this will become clear when we add a phase portrait to the plot.
What about the parameters?
We are just about ready to plot the nullclines. But, first, there is just one more thing we need: approximate values for the many biochemical parameters in the system. We will use the set identified in Süel, et al., Science, 2007, which we will store in a named tuple:
[31]:
CompetenceParams = collections.namedtuple(
"CompetenceParams",
[
"a_k",
"a_s",
"b_k",
"b_s",
"kappa_k",
"kappa_s",
"gamma_k",
"gamma_s",
"n",
"p",
],
)
params = CompetenceParams(
a_k=0.00035,
a_s=0.0,
b_k=0.3,
b_s=3.0,
kappa_k=0.2,
kappa_s=1 / 30,
gamma_k=0.1,
gamma_s=0.1,
n=2,
p=5,
)
The phase portrait of the competence system
At last, using these parameter values, we are able to plot the nullclines:
[62]:
k_range = np.array([1e3, 2])
s_range = np.array([1e3, 5000])
p = plot_nullclines(k_range, s_range, params, log=True)
bokeh.io.show(p)
The nullclines have three crossings, generating three fixed points. To investigate the behavior near the fixed points, and indeed away from them as well, we can plot the full phase portrait.
[66]:
# Turn on zoomability if applicable
if fully_interactive_plots:
zoomable = True
title = None
else:
zoomable = False
title = "Run a live notebook for zoomability"
# Make the plot
p = plot_nullclines(k_range, s_range, params, log=True, title=title)
p = biocircuits.phase_portrait(
dk_dt,
ds_dt,
k_range,
s_range,
(params,),
(params,),
log=True,
p=p,
zoomable=zoomable,
)
bokeh.io.show(p)
How beautiful! Let’s trace the dynamics. We start at the leftmost fixed point, at high concentratons of ComS and low concentrations of ComK. This is the “normal” or “vegetative” state, where the competence system is off and cells can divide and grow normally. One can see that this state is locally stable—nearby trajectories curve in towards it. But if \(k\) were to somehow increase far enough, past the middle fixed point, we can see how the dynamics will drive \(k\) to ever higher concentrations, and then descend to low \(s\), before eventually returning back to the leftmost fixed point. This feature—a large amplitude phase space trajectory following a smaller thresholdcrossing perturbation—is the hallmark of an excitable system.
We will refer to this excursion as the “competence loop.”
To better understand how cells get on the loop, and what happens to them during their voyage, we will now zoom in on the dynamics close to the fixed points. (Because zoomability is not available in the static rendering of this notebook, we will make plots around each of the fixed points individually.)
We start with the leftmost fixed point:
[38]:
k_range = [0.0025, 0.0030]
s_range = [15, 30]
p = plot_nullclines(k_range, s_range, params, log=True)
p = biocircuits.phase_portrait(
dk_dt, ds_dt, k_range, s_range, (params,), (params,), log=True, p=p,
)
bokeh.io.show(p)
This is a stable fixed point, with the system flowing into it from all directions (when \(k\) and \(s\) are close to the fixed point). Let’s now look at the fixed point for the next largest \(k\).
[67]:
k_range = [0.014, 0.02]
s_range = [10, 30]
p = plot_nullclines(k_range, s_range, params, log=True)
p = biocircuits.phase_portrait(
dk_dt, ds_dt, k_range, s_range, (params,), (params,), log=True, p=p,
)
bokeh.io.show(p)
We can see from the flow lines that this is an unstable fixed point, a saddle. We could perform linear stability analysis on this fixed point, and we would find a positive trace and a negative determinant of the linear stability matrix, which means one eigenvalue is positive and the other is negative. This means that the system comes toward the fixed point from one direction, but is pushed away along the other.
Finally, let’s look at the last fixed point, for largest \(k\).
[68]:
k_range = [0.02, 0.05]
s_range = [1, 15]
p = plot_nullclines(k_range, s_range, params, log=True)
p = biocircuits.phase_portrait(
dk_dt, ds_dt, k_range, s_range, (params,), (params,), log=True, p=p
)
bokeh.io.show(p)
This fixed point is a spiral source. It is unstable and the system spirals away from it. This makes competence a transient behavior rather than a stable state.
Numerical solution of dynamical equations
Now that we have a good picture of the dynamics, we can numerically solve the dynamical equations. For simplicity, as an initial condition, we will start with an unphysiological state in which the whole system is “off.” That is, where the concentrations of both ComK and ComS are zero.
[69]:
# Righthandside of ODEs
def rhs(ks, t, p):
"""
Righthand side of dynamical system.
"""
k, s = ks
return np.array([dk_dt(k, s, p), ds_dt(k, s, p)])
# Set up and solve
t = np.linspace(0, 200, 100)
ks_0 = np.array([0, 0])
ks = scipy.integrate.odeint(rhs, ks_0, t, args=(params,))
k, s = ks.transpose()
We will plot this in phase space, with the nullclines overlaid. The trajectory is in gray.
[70]:
k_range = np.array([1e3, 2])
s_range = np.array([1e3, 5000])
p = plot_nullclines(k_range, s_range, params, log=True)
p.line(np.log10(k[1:]), np.log10(s[1:]), line_width=2, color='gray')
bokeh.io.show(p)
We see a rise to the fixed point. Now, let’s say we had some fluctuation that kicked us up a bit off of the fixed point. We’ll start with \(s = 6\) and \(k = 0.05\) and we’ll color the trajectory in a blue color.
[43]:
# Set up and solve
t = np.linspace(0, 200, 1000)
ks_0 = np.array([1e2, 1e3])
ks = scipy.integrate.odeint(rhs, ks_0, t, args=(params,))
k2, s2 = ks.transpose()
p.line(np.log10(k2), np.log10(s2), line_width=2, color='steelblue')
bokeh.io.show(p)
The system takes a much longer trajectory to get to the fixed point, taking the excursion we suspected it would from the phase portrait.
Stochastic simulation
As we have just seen, the system can take a transient journey into competence provided it is perturbed far enough away from the stable fixed point. This raises the question, is intrinsic noise enough to give the “kick” needed to make the excursion to competence? To investigate this, we can perform a Gillespie simulation (SSA) to investigate the dynamics. To specify the master equation, we need the usual ingredients.
A definition of the species present with variables representing their copy number.
A definition of the unit moves that are allowed (the nonzero transition kernels).
Updates to copy numbers for each move.
The propensity for each move.
We can conveniently organize these in tables. First, the variables.
index 
description 
variable 

0 
comK mRNA copy number 

1 
ComK protein copy number 

2 
comS mRNA copy number 

3 
ComS protein copy number 

4 
Free MecA complex copy number 

5 
MecAComK complex copy number 

6 
MecAComS complex copy number 

With the species defined, we an write the move set with propensities.
index 
description 
update 
propensity 

0 
constitutive transcription of a comK mRNA 


1 
regulated transcription of a comK mRNA 


2 
translation of a ComK protein 


3 
constitutive transcription of a comS mRNA 


4 
regulated transcription of a comS mRNA 


5 
translation of a ComS protein 


6 
degradation of a comK mRNA 


7 
degradation of a comK protein 


8 
degradation of a comS mRNA 


9 
degradation of a comS protein 


10 
binding of MecA and ComK 


11 
unbinding of MecA and ComK 


12 
degradation ComK by MecA 


13 
binding of MecA and ComS 


14 
unbinding of MecA and ComS 


15 
degradation ComS by MecA 


In defining the propensities, we introduced parameters. We can define their values, again in accordance with the Süel, et al. paper.
parameter 
value 
units 


1 
– 

1 
– 

1 
– 

1 
– 

1 
molecule/nM 

0.00021875 
1/sec 

0.1875 
1/sec 

0.2 
1/sec 

0 
– 

0.0015 
1/sec 

0.2 
1/sec 

0.005 
1/sec 

0.0001 
1/sec 

0.005 
1/sec 

0.0001 
1/sec 

2.02×10\(^{6}\) 
1/nMsec 

0.0005 
1/sec 

0.05 
1/sec 

4.5×10\(^{6}\) 
1/nMsec 

5×10\(^{5}\) 
1/nMsec 

4×10\(^{5}\) 
1/sec 

5000 
nM 

833 
nM 

2 
– 

5 
– 
Additionally, we take the total MecA complex copy number to be 500. That is, af+ak+as_ = 500
.
If is useful to connect these parameters to those of the ODEs we have considered so far and have used to build the phase portraits. As a reminder, the dimensional version of those ODEs is
\begin{align} \frac{\mathrm{d}K}{\mathrm{d}t}&= \alpha_k + \beta_k\,\frac{K^n}{k_k^n + K^n}  \frac{\delta_k K}{1 + K/\Gamma_k + S/\Gamma_s}  \gamma_k K\\[1em] \frac{\mathrm{d}S}{\mathrm{d}t}&= \alpha_s + \frac{\beta_s}{1 + (K/k_s)^p}  \frac{\delta_s S}{1 + K/\Gamma_k + S/\Gamma_s}  \gamma_s S. \end{align}
By working through the propensities and performing some algebraic manipulations that we will not go through here (though you can on your own if you like), the parameters in the ODEs above are related to those of the discrete model as follows.
\begin{align} &\alpha_k = \frac{k_1\,k_3}{k_7}\,\frac{P_\mathrm{comK}^\mathrm{const}}{\Omega},\\[1em] &\alpha_s = \frac{k_4\,k_6}{k_9}\,\frac{P_\mathrm{comS}^\mathrm{const}}{\Omega},\\[1em] &\beta_k = \frac{k_2\,k_3}{k_7}\,\frac{P_\mathrm{comK}}{\Omega},\\[1em] &\beta_s = \frac{k_5\,k_6}{k_9}\,\frac{P_\mathrm{comS}}{\Omega},\\[1em] &\Gamma_k = \frac{k_{11}+k_{12}}{k_{11}} \\[1em] &\Gamma_s = \frac{k_{13}+k_{14}}{k_{13}} \\[1em] &\delta_k = \frac{k_{12}\,M}{\Gamma_k} \\[1em] &\delta_s = \frac{k_{14}\,M}{\Gamma_s} \\[1em] &\gamma_k = k_8,\\[1em] &\gamma_s = k_{10}, \end{align}
where \(M\) is the total MecA complex copy number (taken to be \(M = 500\)).
We now have all the ingredients, so we can proceed to code up the simulation. That is, we supply the propensity function and updates.
[44]:
@numba.njit
def propensity(
propensities,
population,
t,
P_comK_const,
P_comS_const,
P_comK,
P_comS,
omega,
k1,
k2,
k3,
k4,
k5,
k6,
k7,
k8,
k9,
k10,
k11,
km11,
k12,
k13,
km13,
k14,
kappa_k,
kappa_s,
n,
p,
):
mk, pk, ms, ps, af, ak, as_ = population
f_num = (pk / kappa_k * omega) ** n
f = k2 * f_num / (1 + f_num)
g = k5 / (1 + (pk / kappa_s * omega) ** p)
propensities[0] = k1 * P_comK_const
propensities[1] = P_comK * f
propensities[2] = k3 * mk
propensities[3] = k4 * P_comS_const
propensities[4] = P_comS * g
propensities[5] = k6 * ms
propensities[6] = k7 * mk
propensities[7] = k8 * pk
propensities[8] = k9 * ms
propensities[9] = k10 * ps
propensities[10] = k11 * af * pk / omega
propensities[11] = km11 * ak
propensities[12] = k12 * ak
propensities[13] = k13 * af * ps / omega
propensities[14] = km13 * as_
propensities[15] = k14 * as_
update = np.array(
[
# 0 1 2 3 4 5 6
[ 1, 0, 0, 0, 0, 0, 0], # 0
[ 1, 0, 0, 0, 0, 0, 0], # 1
[ 0, 1, 0, 0, 0, 0, 0], # 2
[ 0, 0, 1, 0, 0, 0, 0], # 3
[ 0, 0, 1, 0, 0, 0, 0], # 4
[ 0, 0, 0, 1, 0, 0, 0], # 5
[1, 0, 0, 0, 0, 0, 0], # 6
[ 0, 1, 0, 0, 0, 0, 0], # 7
[ 0, 0, 1, 0, 0, 0, 0], # 8
[ 0, 0, 0, 1, 0, 0, 0], # 9
[ 0, 1, 0, 0, 1, 1, 0], # 10
[ 0, 1, 0, 0, 1, 1, 0], # 11
[ 0, 0, 0, 0, 1, 1, 0], # 12
[ 0, 0, 0, 1, 1, 0, 1], # 13
[ 0, 0, 0, 1, 1, 0, 1], # 14
[ 0, 0, 0, 0, 1, 0, 1], # 15
],
dtype=np.int64,
)
Next, for convenience, we’ll write a function that will return the arguments for the propensity funtion, complete with defaults from the Süel, et al, 2007 paper. We will store the result in a named tuple for easy access of the parameter values when we need them.
[45]:
ParametersSSA = collections.namedtuple(
"ParametersSSA",
[
"P_comK_const",
"P_comS_const",
"P_comK",
"P_comS",
"omega",
"k1",
"k2",
"k3",
"k4",
"k5",
"k6",
"k7",
"k8",
"k9",
"k10",
"k11",
"km11",
"k12",
"k13",
"km13",
"k14",
"kappa_k",
"kappa_s",
"n",
"p",
],
)
params_ssa = ParametersSSA(
P_comK_const=1,
P_comS_const=1,
P_comK=1,
P_comS=1,
omega=1,
k1=0.00021875,
k2=0.1875,
k3=0.2,
k4=0.0,
k5=0.0015,
k6=0.2,
k7=0.005,
k8=0.0001,
k9=0.005,
k10=0.0001,
k11=2.02e6,
km11=0.0005,
k12=0.05,
k13=4.5e6,
km13=5e5,
k14=4e5,
kappa_k=5000.0,
kappa_s=833.0,
n=2,
p=5,
)
Next, we need to write a function to convert the counts of ComK and ComS protein to dimensionless concentration units.
[46]:
def convert_to_conc(pop, params):
k = pop[0, :, 1] / params.omega * params.k11 / (params.km11 + params.k12)
s = pop[0, :, 3] / params.omega * params.k13 / (params.km13 + params.k14)
return k, s
Finally, we’ll set up the initial populations and time points.
[47]:
population_0 = np.array([10, 100, 10, 400, 500, 0, 0], dtype=int)
time_points = np.linspace(0, 4000000, 4001)
All of the required arguments are set, so let’s now run a trajectory. When we plot the trajectory, we will show it along with the phase protrait (with counts adjusted to the same dimensionless concentration units).
[48]:
# Perform the Gillespie simulation
pop = biocircuits.gillespie_ssa(
propensity, update, population_0, time_points, args=tuple(params_ssa)
)
k, s = convert_to_conc(pop, params_ssa)
# Make plot
k_range = np.array([1e4, 2])
s_range = np.array([1e2, 5000])
# Plot trajectory, making sure to handle zeros in loglog plot
k_plot = k[(k > 0) & (s > 0)]
s_plot = s[(k > 0) & (s > 0)]
p = bokeh.plotting.figure(
frame_width=260,
frame_height=260,
x_axis_label="log₁₀ k",
y_axis_label="log₁₀ s",
x_range=k_range,
y_range=s_range,
)
p.line(np.log10(k_plot), np.log10(s_plot), color="gray")
# Nullclines and phase portrait
p = plot_nullclines(k_range, s_range, params, log=True, p=p)
p = biocircuits.phase_portrait(
dk_dt,
ds_dt,
k_range,
s_range,
(params,),
(params,),
log=True,
p=p,
)
bokeh.io.show(p)