# Colab setup ------------------
import os, sys, subprocess
if "google.colab" in sys.modules:
cmd = "pip install --upgrade biocircuits watermark"
process = subprocess.Popen(cmd.split(), stdout=subprocess.PIPE, stderr=subprocess.PIPE)
stdout, stderr = process.communicate()
# ------------------------------
import numpy as np
import biocircuits
import bokeh.io
import bokeh.layouts
import bokeh.models
import bokeh.plotting
bokeh.io.output_notebook()In the previous technical appendix, we analyzed signaling cascades with doubly-phosphorylated molecules and two levels of a cascade. Here, we model the general case where each step in the cascade may have arbitrarily many phophorylation sites and the cascade may have arbitrarily many levels. Again, the maximally phorphorylated signaling molecule at the end of the cascade leads to a change in gene expression in response to an inputted signal. We will again restrict ourselves to the case where all enzyme-catalyzed reactions are far from saturation and phosphorylation rate constants are the same, as are dephosphorylation rate constants.
As in the previous technical appendix, we will approximate the rate of phosphorylation of species X catalyzed by species S as
\[ \begin{align} \text{rate of phosphorylation} \approx k_+ s x, \end{align} \tag{1}\]
and the dephosphorylation of species X as
\[ \begin{align} \text{rate of dephosphorylation} \approx k_- x. \end{align} \tag{2}\]
If a signaling molecule X has \(n\) sites that may be phosophorylated, each catalyzed by S, we have a reaction scheme
\[ \begin{align} \require{mhchem} \ce{X_0 <=>[\mathrm{S}] X_1 <=>[\mathrm{S}] X_2 $\cdots$ <=>[\mathrm{S}] X$_\mathrm{n}$}, \end{align} \tag{3}\]
where the subscripts represent the degree of phosphorylation. We can write down a system of ODEs describing the dynamics of each species.
\[ \begin{align} &\frac{\mathrm{d}x_0}{\mathrm{d}t} = - k_+ s x_0 + k_- x_1,\\[1em] &\frac{\mathrm{d}x_i}{\mathrm{d}t} = - k_+ s x_i + k_- x_{i+1} + k_+sx_{i-1} - k_- x_i\;\;\text{for } 1 \le i \le n-1,\\[1em] &\frac{\mathrm{d}x_n}{\mathrm{d}t} = k_+ s x_{n-1} - k_- x_n. \end{align} \tag{4}\]
Similarly as in the previous technical appendix, we can nondimensionalize by scaling all concentrations by \(K = k_-/k_+\) and time by \(1/k_-\), to give
\[ \begin{align} &\frac{\mathrm{d}x_0}{\mathrm{d}t} = - s x_0 + x_1,\\[1em] &\frac{\mathrm{d}x_i}{\mathrm{d}t} = - s x_i + x_{i+1} + sx_{i-1} - x_i\;\;\text{for } 1 \le i \le n-1,\\[1em] &\frac{\mathrm{d}x_n}{\mathrm{d}t} = s x_{n-1} - x_n. \end{align} \tag{5}\]
We will henceforth always work with dimensionless quantities. We will not solve these equations here, but rather will do numerical solutions of these equations when coupled to other reactions in a cascade.
We wish to compute the dimensionless steady state concentrations of each species \(x_0, \ldots x_n\). We accomplish this by noting at steady state the forward and reverse rate of each reaction are equal. Thus, we have
\[ \begin{align} sx_{i-1} = x_i\;\;\text{for } 1 \le i \le n, \end{align} \tag{6}\]
Therefore, for \(i = n\),
\[ \begin{align} x_{n-1} = \frac{x_n}{s}. \end{align} \tag{7}\]
For \(i = n-1\), we have
\[ \begin{align} x_{n-2} = \frac{x_{n-1}}{s} = \frac{x_n}{s^2}, \end{align} \tag{8}\]
where we have substituted the above expression for \(x_{n-1}\). We continue this process all the way until \(i = 0\) to get the result for all \(i \in [0, n]\).
\[ \begin{align} x_i = s^{-(n-i)}\,x_n. \end{align} \tag{9}\]
To determine \(x_n\) (and therefore all of the \(x_i\)’s via the above equation for \(x_i\)), we note that the total amount of X molecules of all degrees of phosphorylation is conserved. Letting \(x_\mathrm{tot}\) be the total X concentration, we have
\[ \begin{align} x_\mathrm{tot} = \sum_i x_i = \sum_{i=0}^n s^{-(n-i)}\,x_n = \frac{x_n}{s^n}\sum_{i=0}^n s^i. \end{align} \tag{10}\]
We recognize the sum as a partial sum of a geometric series such that
\[ \begin{align} \sum_{i=0}^n s^i = \frac{1 - s^{n+1}}{1 - s}. \end{align} \tag{11}\]
Therefore,
\[ \begin{align} x_\mathrm{tot} = \frac{1 - s^{n+1}}{s^n(1-s)}\,x_n. \end{align} \tag{12}\]
yielding
\[ \begin{align} x_n = x_\mathrm{tot}\,\frac{s^n(1-s)}{1 - s^{n+1}}. \end{align} \tag{13}\]
with the special case of \(x_n = x_\mathrm{tot}/(1+n)\) for \(s = 1\).
For convenience in coding up a plot and for analysis we will do momentarily, we define
\[ \begin{align} f_n(s) = \left\{\begin{array}{cl} \displaystyle{\frac{s^n(1-s)}{1-s^{n+1}}} & \text{for }s \ne 1 \\[1em] (1+n)^{-1} & \text{for }s = 1. \end{array} \right. \end{align} \tag{14}\]
Note that in the limit of large \(n\),
\[ \begin{align} f_{\infty}(s) = \lim_{n\to\infty} \frac{s^n(1-s)}{1-s^{n+1}} = \left\{ \begin{array}{ll} 0 & \text{for } s \le 1,\\[1em] 1-1/s & \text{for } s > 1.\\[1em] \end{array} \right. \end{align} \tag{15}\]
With \(f_n(s)\) so defined, are have,
\[ \begin{align} x_n = x_\mathrm{tot}\,f_n(s). \end{align} \tag{16}\]
Since \(x_n\) is our quantity of interest, let us look at how it varies with signal intensity.
def f_infinite_n(s):
# The result
res = np.empty_like(s)
# Indices where s greater than one
inds = s > 1
# Compute result
res[~inds] = 0.0
res[inds] = (s[inds] - 1) / s[inds]
return res
def f(s, n):
if np.isinf(n):
res = f_infinite_n(s)
else:
res = np.empty_like(s)
# Indices where s is not one
inds = ~np.isclose(s, 1)
# Compute result for close to one and otherwise
res[inds] = s[inds] ** n * (1 - s[inds]) / (1 - s[inds] ** (n + 1))
res[~inds] = s[~inds] ** n / (1 + n)
return res
# Values of s and n for which we want a plot
s = np.logspace(-2, 2, 200)
s_near_1 = np.logspace(0.0001, np.log10(1.011), 25)
s = np.sort(np.concatenate((s, s_near_1)))
n = [1, 2, 3, 4, np.inf]
# Styling for figures
fig_kwargs = dict(
frame_width=400,
frame_height=200,
x_axis_type="log",
x_axis_label=r"$$s/K$$",
y_axis_label=r"$$x_n/x_\mathrm{tot}$$",
x_range=[s.min(), s.max()],
)
# Linear and log scale figures
p = bokeh.plotting.figure(**fig_kwargs)
p_log = bokeh.plotting.figure(
**fig_kwargs, y_axis_type="log", y_range=[1e-8, 3], visible=False
)
# Select log or linear y-scale
radio_button_group = bokeh.models.RadioButtonGroup(
labels=["log", "linear"], active=1, width=60, orientation="vertical"
)
radio_button_group.js_on_change(
"active",
bokeh.models.CustomJS(
args=dict(p_log=p_log, p=p),
code="""
if (p_log.visible == true) {
p_log.visible = false;
p.visible = true;
}
else {
p_log.visible = true;
p.visible = false;
}
""",
),
)
# Populate glyphs
palette = list(bokeh.palettes.Blues7[:4][::-1]) + ["black"]
for n_val, color in zip(n, palette):
label = "n → ∞" if np.isinf(n_val) else f"n = {n_val}"
yn = f(s, n_val)
p.line(s, yn, line_width=2, color=color, legend_label=label)
p_log.line(s, yn, line_width=2, color=color, legend_label=label)
# Handle zero values on log plot
if np.isinf(n_val):
s1_ind = np.where(s >= 1)[0][0]
p_log.ray(1.0, yn[s1_ind], angle=-np.pi / 2, line_width=2, color=color)
# Position legend
p.legend.location = "bottom_right"
p_log.legend.location = "bottom_right"
# Build and display layout
layout = bokeh.layouts.row(
bokeh.models.Spacer(width=10),
bokeh.layouts.column(bokeh.models.Spacer(height=75), radio_button_group),
bokeh.models.Spacer(width=10),
bokeh.layouts.column(p, p_log),
)
bokeh.io.show(layout)We have thus far considered a single multi-phosphorylation reaction. Now, imagine that the final product of this reaction catalyzes a separate multi-phosphorylation reaction. We label the first signaling molecule X⁽¹⁾ with a superscript (1) and the second, X⁽²⁾ with a superscript (2). Signaling molecule (1) has \(n_1\) phosphorylation sites and signaling molecule (2) has \(n_2\).
\[\begin{align} \require{mhchem} &\ce{X_0^{(1)} <=>[\mathrm{S}] X_1^{(1)} <=>[\mathrm{S}] X_2^{(1)} $\cdots$ <=>[\mathrm{S}] X$_\mathrm{n_1}$^{(1)}},\\[1em] &\ce{X_0^{(2)} <=>[\mathrm{X_{n_1}^{(1)}}] X_1^{(2)} <=>[\mathrm{X_{n_1}^{(1)}}] X_2^{(2)} $\cdots$ <=>[\mathrm{X_{n_1}^{(1)}}] X$_\mathrm{n_2}$^{(2)}}. \end{align}\]
We can extend this to a cascade of \(m\) signaling molecules.
\[\begin{align} \require{mhchem} &\ce{X_0^{(1)} <=>[\mathrm{S}] X_1^{(1)} <=>[\mathrm{S}] X_2^{(1)} $\cdots$ <=>[\mathrm{S}] X$_\mathrm{n_1}$^{(1)}}\\[1em] &\quad\quad\quad\quad\quad\vdots\\[1em] &\ce{X_0^{(j)} <=>[\mathrm{X_{n_{j-1}}^{(j-1)}}] X_1^{(j)} <=>[\mathrm{X_{n_{j-1}}^{(j-1)}}] X_2^{(j)} $\cdots$ <=>[\mathrm{X_{n_{j-1}}^{(j-1)}}] X$_\mathrm{n_i}$^{(j)}}\\[1em] &\quad\quad\quad\quad\quad\vdots\\[1em] &\ce{X_0^{(m)} <=>[\mathrm{X_{n_{m-1}}^{(m-1)}}] X_1^{(m)} <=>[\mathrm{X_{n_{m-1}}^{(m-1)}}] X_2^{(m)} $\cdots$ <=>[\mathrm{X_{n_{m-1}}^{(m-1)}}] X$_\mathrm{n_2}$^{(m)}}. \end{align}\]
We can use the same arguments as above to write down the dimensionless dynamical equations (again with all concentrations in units of \(K\) and time in units of \(1/k_-\)) for the concentration of species in this set of reactions.
\[\begin{align} &\frac{\mathrm{d}x^{(1)}_0}{\mathrm{d}t} = - s x_0^{(1)} + x_1^{(1)},\\[1em] &\frac{\mathrm{d}x^{(1)}_i}{\mathrm{d}t} = - s x^{(1)}_i + x^{(1)}_{i+1} + sx^{(1)}_{i-1} - x^{(1)}_i\;\;\text{for } 1 \le i \le n_1-1,\\[1em] &\frac{\mathrm{d}x^{(1)}_{n_1}}{\mathrm{d}t} = s x^{(1)}_{n_1-1} - x^{(1)}_{n_1},\\[1em] &\frac{\mathrm{d}x^{(j)}_0}{\mathrm{d}t} = - x^{(j-1)}_{n_{j-1}} x_0 + x^{(j)}_1\;\;\text{for }2 \le j \le m\\[1em] &\frac{\mathrm{d}x^{(j)}_i}{\mathrm{d}t} = - x^{(j-1)}_{n_{j-1}} x^{(j)}_i + x^{(j)}_{i+1} + x^{(j-1)}_{n_{j-1}} x^{(j)}_{i-1} - x^{(j)}_i\;\;\text{for } 1 \le i \le n-1,\;2 \le j \le m\\[1em] &\frac{\mathrm{d}x^{(j)}_{n_{j}}}{\mathrm{d}t} = x^{(j-1)}_{n_{j-1}} x^{(j)}_{n_j-1} - x^{(j)}_{n_{j}}\;\;\text{for }2 \le j \le m. \end{align}\]
Our goal is to study the dynamics and steady state of \(x^{(m)}_{n_m}\) (the steady state being the transfer function of the signaling circuit), since this is the end of the signaling cascade leading ultimately to regulation of gene expression. In looking at the above equations, we see that the equations for \(x^{(m)}_{n_m}\) dynamics are analogous to those of \(x_n\) in the case we have already worked out for a single multiphosphorlyation reaction, except with \(x^{(m-1)}_{n_{m-1}}\) replacing \(s\). So, we have, at steady state,
\[\begin{align} x^{(m)}_{n_m} = x^{(m)}_\mathrm{tot}f_{n_m}\left(x^{(m-1)}_{n_{m-1}}\right). \end{align}\]
But by the same analogy,
\[\begin{align} x^{(m-1)}_{n_{m-1}} = x^{(m-1)}_\mathrm{tot}f_{n_{m-1}}\left(x^{(m-2)}_{n_{m-2}}\right), \end{align}\]
giving
\[\begin{align} x^{(m)}_{n_m} = x^{(m)}_\mathrm{tot}f_{n_m}\left(x^{(m-1)}_\mathrm{tot}f_{n_{m-1}}\left(x^{(m-2)}_{n_{m-2}}\right)\right). \end{align}\]
We continue up to the top of the cascade, giving
\[\begin{align} x^{(m)}_{n_m} = x^{(m)}_\mathrm{tot}f_{n_m}\left(x^{(m-1)}_\mathrm{tot}f_{n_{m-1}}\left(x^{(m-2)}_\mathrm{tot}f_{n_{m-2}}\left(\cdots x^{(1)}_\mathrm{tot} f_{n_1}\left(s\right)\right) \right)\right). \end{align}\]
It is easy to get confused reading this not-quite-composite (because of the presence of the \(x^{(\cdot)}_\mathrm{tot}\)’s) function, so here is the function for \(m = 1, 2\), and \(3\) as examples.
\[\begin{align} &x^{(1)}_{n_1} = x^{(1)}_\mathrm{tot}\,f_{n_1}\left(s\right),\\[1em] &x^{(2)}_{n_2} = x^{(2)}_\mathrm{tot}\,f_{n_2}\left(x^{(1)}_\mathrm{tot}\,f_{n_1}(s)\right),\\[1em] &x^{(3)}_{n_3} = x^{(3)}_\mathrm{tot}\,f_{n_3}\left(x^{(2)}_\mathrm{tot}\,f_{n_2}\left(x^{(1)}_\mathrm{tot}\,f_{n_1}(s)\right)\right). \end{align}\]
We have written down the dynamical equations and worked out the steady state transfer function for a multi-level signaling cascade. To gain insight to the nature of this signaling circuit, we wish to quantify
In order assess the speed of the response of the circuit, we need to solve the dynamical equations numerically, which we will do momentarily, and also compute the derivative of the transfer function. For the latter task, we can take advantage of the nested form of the transfer function using the chain rule. We define
\[\begin{align} f_n'(s) = \frac{\mathrm{d}f_n}{\mathrm{d}s} = \left\{\begin{array}{ll} \displaystyle{\frac{n}{2(1+n)}} & \text{for }s = 1,\\[1em] \displaystyle{\frac{s^{n-1}}{\left(1-s^{n+1}\right)^2}\left(n(1-s)-s\left(1-s^n\right)\right)} & \text{otherwise}. \end{array} \right. \end{align}\]
For \(n\to\infty\), this is
\[\begin{align} f_\infty'(s) = \frac{\mathrm{d}f_\infty}{\mathrm{d}s} = \left\{ \begin{array}{ll} 0 & \text{for } s < 1,\\[1em] 1/2 & \text{for } s = 1,\\[1em] 1/s^2 & \text{for } s > 1.\\[1em] \end{array} \right. \end{align}\]
With this definition, we can differentiate the transfer function using the chain rule as
\[\begin{align} \frac{\mathrm{d} x^{(m)}_{n_m}}{\mathrm{d}s} = \left(\prod_{j=1}^m x^{(j)}_\mathrm{tot}\right)f_{n_1}'(s) \cdot f_{n_2}'\left(x^{(1)}_\mathrm{tot}f_{n_1}(s)\right) \cdots f_{n_m}'\left(x^{(m-1)}_\mathrm{tot}f_{n_{m-1}}\left(x^{(m-2)}_\mathrm{tot}f_{n_{m-2}}\left(\cdots x^{(1)}_\mathrm{tot} f_{n_1}(s) \right)\right)\right). \end{align}\]
This expression is, again, eye-crossing. Here are the derivatives for \(m = 1, 2\), and \(3\).
\[\begin{align} &\frac{\mathrm{d} x^{(1)}_{n_1}}{\mathrm{d}s} = x^{(1)}_\mathrm{tot}\,f_{n_1}'(s),\\[1em] &\frac{\mathrm{d} x^{(2)}_{n_2}}{\mathrm{d}s} = x^{(1)}_\mathrm{tot}\,x^{(2)}_\mathrm{tot}\,f_{n_1}'(s) \cdot f_{n_2}'\left(x^{(1)}_\mathrm{tot}f_{n_1}(s)\right), \\[1em] &\frac{\mathrm{d} x^{(3)}_{n_3}}{\mathrm{d}s} = x^{(1)}_\mathrm{tot}\,x^{(2)}_\mathrm{tot}\,x^{(3)}_\mathrm{tot}\,f_{n_1}'(s) \cdot f_{n_2}'\left(x^{(1)}_\mathrm{tot}f_{n_1}(s)\right) \cdot f_{n_3}'\left(x^{(2)}_\mathrm{tot}\,f_{n_2}\left(x^{(1)}_\mathrm{tot}f_{n_1}(s)\right)\right). \end{align}\]
As before, with these transfer functions and their derivatives in hand, we can compute the sensitivity as
\[\begin{align} \text{sensitivity} = \frac{\mathrm{d}\ln x^{(m)}_{n_m}}{\mathrm{d} \ln s} = \frac{s}{x^{(m)}_{n_m}}\,\frac{\mathrm{d}x^{(m)}_{n_m}}{\mathrm{d}s}, \end{align}\]
and the gain as \(x^{(m)}_{n_m} / s\).
As is often the case, it is valuable to investigate limiting behavior of the signaling circuit. First, we consider \(s \ll 1\) (that is that the input signal concentration is much less than \(K = k_- / k_+\)). In this regime, for finite \(n\),
\[\begin{align} f_n(s) = \frac{s^n(1-s)}{1-s^{n+1}} \approx s^n, \end{align}\]
so that
\[\begin{align} x^{(m)}_{n_m} \approx s^{\prod_{j=1}^m n_j} \prod_{j=1}^m \left(x_\mathrm{tot}^{(j)}\right)^{\prod_{k=j+1}^m n_k}. \end{align}\]
So, the transfer function can decay very rapidly with decreasing input signal, e.g. like \(s^{-27}\) for a three-level cascade of triply phosphorylated signaling molecules. This leads to a threshold-like behavior. The input signal evokes almost no response from the signaling circuit until the input signal is of magnitude \(K\). However, if the signaling molecules only have a single phosphorylation site, such that \(n_j = 1\) for all levels of the cascade, then
\[\begin{align} x^{(m)}_{n_m} \approx s \prod_{j=1}^m x_\mathrm{tot}^{(j)}, \end{align}\]
meaning that at low input signal, the gain is given by the product of the concentrations of the signaling molecules in the cascade.
Still considering the small input signal limit, the derivative of the transfer function is
\[\begin{align} \frac{\mathrm{d} x^{(m)}_{n_m}}{\mathrm{d}s} \approx s^{\left(\prod_{j=1}^m n_j\right)-1} \prod_{j=1}^m n_j\left(x_\mathrm{tot}^{(j)}\right)^{\prod_{k=j+1}^m n_k} = \frac{x^{(m)}_{n_m}}{s}\,\prod_{j=1}^m n_j, \end{align}\]
from which we can compute the sensitivity for \(s \ll 1\) as
\[\begin{align} \text{sensitivity} = \prod_{j=1}^m n_j. \end{align}\]
Consider now the case where \(s\) is large (dimensional input signal must larger than \(K\)). In this case,
\[\begin{align} f_n(s) = \frac{s^n(1-s)}{1-s^{n+1}} \approx \frac{s^n(s-1)}{s^{n+1}} = 1 - \frac{1}{s}. \end{align}\]
Thus, for large \(s\), the dependence on the number of phosphorylation sites is eliminated and the transfer function is the same as for all \(n_j\to\infty\). We can investigate limiting behavior for large \(s\) be defining \(\xi = 1/s\) and writing the transfer function as a function of \(\xi\). The result is
\[\begin{align} x^{(m)}_{n_m} = x^{(m)}_\mathrm{tot} \left(1 - x^{(m-1)}_\mathrm{tot}\left(1 - x^{(m-2)}_\mathrm{tot} \cdots \left(1 - x^{(1)}_\mathrm{tot}(1-\xi)\right)\right)\right). \end{align}\]
Again for clarity, the approximate transfer functions for \(m = 1, 2\), and \(3\) are
\[\begin{align} &x^{(1)}_{n_1} = x^{(1)}_\mathrm{tot}(1 - \xi), \\[1em] &x^{(2)}_{n_2} = x^{(2)}_\mathrm{tot}\left(1 - x^{(1)}_\mathrm{tot}(1 - \xi)\right), \\[1em] &x^{(3)}_{n_3} = x^{(3)}_\mathrm{tot}\left(1 - x^{(2)}_\mathrm{tot}\left(1 - x^{(1)}_\mathrm{tot}(1 - \xi)\right)\right). \end{align}\]
It is clear that the transfer function is linear in \(\xi\) such that the following scaling holds.
\[\begin{align} \frac{x^{(m)}_{n_m}}{x^{(m)}_\mathrm{tot}} - 1 \sim \frac{1}{s}. \end{align}\]
This means that for very large input signal, nearly all of the signaling molecules are completely phosphorylated and the signal is given by how many molecules for the last level of the cascade are available.
According to the scaling of the transfer function, its derivative scales like
\[\begin{align} \frac{\mathrm{d} x^{(m)}_{n_m}}{\mathrm{d}s} \sim -\frac{1}{s^2}, \end{align}\]
and the sensitivity therefore scales like \(1/s\) for large \(s\). To consider the gain, we divide both sides of the scaling relation of the transfer function by \(s\) to give
\[\begin{align} \frac{x^{(m)}_{n_m}/s}{x^{(m)}_\mathrm{tot}} - \frac{1}{s} \sim \frac{1}{s^2}. \end{align}\]
Rearranging,
\[\begin{align} \frac{x^{(m)}_{n_m}/s}{x^{(m)}_\mathrm{tot}} \sim \frac{1}{s}, \end{align}\]
where we have noted that the \(1/s\) term dominates the \(1/s^2\) term for large \(s\). Thus,
\[\begin{align} \text{gain} = x^{(m)}_{n_m}/s \sim \frac{1}{s}. \end{align}\]
In summary, in the large \(s\) limit, the transfer function approaches unity, the gain decays like \(1/s\), as does the sensitivity, and the derivative of the transfer function scales like \(-1/s^2\).
The above expressions for the transfer function, its derivative, the gain, and the sensitivity are challenging to code, and we will not show the details here. Rather, the interactive plot below shows how these functions vary with dimensionless input signal \(s\). The curves are colored for \(n = 1, 2, 3,\) and \(4\) going from light to dark blue. The gray curve represents the baseline response corresponding to the case where an input signal catalyzes a single phosphorylation of a single signaling molecule. (Note that some discontinuities will appear in plots of sensitivity for large cascade depths for small \(s\). This is an artifact due to machine precision and underflow issues.)
bokeh.io.show(biocircuits.jsplots.unsat_phosphorylation_signal_cascade())%load_ext watermark
%watermark -v -p numpy,bokeh,biocircuits,jupyterlab