Technical Appendix 7a: Signaling cascades

Open in Google Colab | Download notebook

# Colab setup ------------------
import os, sys, subprocess
if "google.colab" in sys.modules:
    cmd = "pip install --upgrade colorcet 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 colorcet

bokeh.io.output_notebook()

In this technical appendix, we work out the details of modeling a simple signaling cascade where a signal catalytically double-phosphorylates a signaling molecule. The doubly-phosphorylated signaling molecule catalyzes double-phosphoryltation of a second signaling molecule, whose doubly-phorphorylated form ultimately serves a regulatory role. We will build the mathematical expressions step by step.

Phosphorylation catalyzed by a signaling molecule

Let us consider the phosphorylation of a signaling molecule X catalyzed by a signaling molecule S.

\[ \begin{align} \require{mhchem} \ce{X_0<=>[\mathrm{S}] X_1}. \end{align} \tag{1}\]

Here, the subscript on the chemical species represents the number of phosphorylated sites on the molecule. To ascribe a rate of phosphorylation, we use a Michaelis-Menten expression for the enzyme-catalyzed reaction.

\[ \begin{align} \text{rate of phosphorylation} = k_\mathrm{cat} s \,\frac{x_0/K_M}{1 + x_0/K_M}, \end{align} \tag{2}\]

where \(K_M\) is the Michaelis constant. We will consider the limit where the concentration of X is small, \(x_0 \ll K_M\), such that

\[ \begin{align} \text{rate of phosphorylation} \approx k_+ s x_0, \end{align} \tag{3}\]

where \(k_+ = k_\mathrm{cat}/K_M\).

The total amount of X is conserved; letting \(x_\mathrm{tot}\) be the total X concentration, we have \(x_\mathrm{tot} = x_0 + x_1\). We can write the dynamics of \(x_1\) as

\[ \begin{align} \frac{\mathrm{d}x_1}{\mathrm{d}t} = k_+ s (x_\mathrm{tot} - x_1) - k_- x_1. \end{align} \tag{4}\]

which is solved by integrating factor to give

\[ \begin{align} x_1(t) = x_1^\mathrm{init}\,\mathrm{e}^{-(k_+s + k_-)t} + x_\mathrm{tot}\,\frac{k_+ s}{k_+ s + k_-}\left(1 - \mathrm{e}^{-(k_+s + k_-)t}\right), \end{align} \tag{5}\]

where \(x_1^\mathrm{init}\) is the initial concentration of \(x_1\). Evidently, the time scale of the response to a step is signal is \(1/(k_+s + k_-)\).

The steady state is achieved as \(t\) goes toward infinity, giving

\[ \begin{align} x_1 = x_\mathrm{tot}\,\frac{k_+ s}{k_+ s + k_-} = x_\mathrm{tot}\,\frac{s/K}{1 + s/K}, \end{align} \tag{6}\]

where \(K \equiv k_- / k_+\), which is akin to a dissociation constant, but is not a dissociation constant because it is related to a ratio of out-of-equilibrium rates.

Catalyzed double-phosphorlyation

We now apply a similar analysis to the case where the molecule X may be phosphorylated twice, stating with the dynamics.

Dynamics of catalyzed two-state phosphorlyation

In this case, we have a reaction scheme

\[ \begin{align} \require{mhchem} \ce{X_0 <=>[\mathrm{S}] X_1 <=>[\mathrm{S}] X_2}. \end{align} \tag{7}\]

For simplicity, we again assume that the catalyzed reactions operate far from saturayion and that the \(k_+\) and \(k_-\) are the same for both phosphorylation and dephosphorylation reactions, respectively. We can write down a system of ODEs for the concentrations 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_1}{\mathrm{d}t} = - k_+ s x_1 + k_- x_2 + k_+sx_0 - k_- x_1,\\[1em] &\frac{\mathrm{d}x_n}{\mathrm{d}t} = k_+ s x_{1} - k_- x_2. \end{align} \tag{8}\]

This linear system has an analytical solution, but we will not explore that here, but rather will do numerical solutions when coupled to other reactions in a cascade.

Steady state of catalyzed two-step phosphorylation

If the time scale of phosphorylation in a signaling response is fast compared to cellular response (e.g., regulation of gene expression), then the steady state of the phosphorylation process dictates the cell response and is therefore of particular interest. We can compute the steady state concentrations of each species \(x_0\), \(x_1\), and \(x_2\). We accomplish this by noting at steady state the forward and reverse rate of each reaction are equal. Thus, we have

\[ \begin{align} &k_+ s x_0 = k_- x_1,\\[1em] &k_+sx_1 = k_- x_2. \end{align} \tag{9}\]

The second expression gives

\[ \begin{align} x_{1} = \frac{k_-}{k_+ s}\,x_2 = \frac{x_2}{s/K}. \end{align} \tag{10}\]

We can then write

\[ \begin{align} x_{1} = \frac{K}{s}\,x_2. \end{align} \tag{11}\]

Similarly, we have

\[ \begin{align} x_{0} = \frac{K}{s}\,x_{1} = \left(\frac{K}{s}\right)^2\,x_2, \end{align} \tag{12}\]

where we have substituted the above expression for \(x_{1}\). To determine \(x_2\) (and therefore also \(x_1\) and \(x_0\) via the above equations), 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} = x_0 + x_1 + x_2 = x_2\left(1 + \frac{K}{s} + \left(\frac{K}{s}\right)^2\right). \end{align} \tag{13}\]

We can solve for \(x_2\) to get

\[ \begin{align} x_2 = x_\mathrm{tot}\,\frac{(s/K)^2}{1 + s/K + (s/K)^2}. \end{align} \tag{14}\]

Two simple cascades

We have worked out the dynamical equations and steady states for catalyzed single- and double-phoshorylation. We now consider linking the product of these reactions to another phorphorylation. In the simplest example, we can link a single phosphorylation to another single phosphorylation.

\[ \begin{align} \require{mhchem} &\ce{X_0<=>[\mathrm{S}] X_1},\\[1em] &\ce{Y_0<=>[\mathrm{X}_1] Y_1}. \end{align} \tag{15}\]

In this case, the output of the circuit is the concentration of Y₁. We will again, and henceforth, assume that all \(k_+\) values are the same, as are the \(k_-\) values. At this point, it is worth using dimensionless variables, taking \(t \leftarrow k_-t\), \(x_i \leftarrow x_i/K\), \(y_i \leftarrow y_i/K\), and \(s \leftarrow s/K\). That is to say, we nondimensionalize time using \(k_-\) and all concentrations with \(K = k_-/k_+\), which we assumed to be the same for all reactions for simplicity. The dimensionless dynamical equations are

\[ \begin{align} &\frac{\mathrm{d}x_1}{\mathrm{d}t} = s (x_\mathrm{tot} - x_1) - x_1,\\[1em] &\frac{\mathrm{d}y_1}{\mathrm{d}t} = x_1 (y_\mathrm{tot} - y_1) - y_1. \end{align} \tag{16}\]

Note here that \(x_\mathrm{tot}\) and \(y_\mathrm{tot}\) are dimensionless, given in units of \(K\). We can work out that the steady state dimensionless concentration of \(y_1\) is

\[ \begin{align} y_1 = y_\mathrm{tot}\,\frac{x_\mathrm{tot}\,s}{1 + s + x_\mathrm{tot}\,s}. \end{align} \tag{17}\]

Now imagine that both X and Y can undergo double-phosphorylation such that doubly phosphorylated X catalyzes phosphorylation of Y, and doubly phosphorylated Y is the output of the circuit. The reaction system is

\[ \begin{align} \require{mhchem} &\ce{X_0 <=>[\mathrm{S}] X_1 <=>[\mathrm{S}] X_2},\\[1em] &\ce{Y_0 <=>[\mathrm{X}_2] Y_1 <=>[\mathrm{X}_2] Y_2}. \end{align} \tag{18}\]

The dimensionless dynamical equations are

\[ \begin{align} &\frac{\mathrm{d}x_0}{\mathrm{d}t} = - s x_0 + x_1,\\[1em] &\frac{\mathrm{d}x_1}{\mathrm{d}t} = - s x_1 + x_{2} + sx_{0} - x_1,\\[1em] &\frac{\mathrm{d}x_2}{\mathrm{d}t} = s x_{1} - x_2,\\[1em] &\frac{\mathrm{d}y_0}{\mathrm{d}t} = - x_2 y_0 + y_1,\\[1em] &\frac{\mathrm{d}y_1}{\mathrm{d}t} = - x_2 y_1 + y_{2} + x_2 y_{0} - y_1,\\[1em] &\frac{\mathrm{d}y_2}{\mathrm{d}t} = x_2 y_{1} - y_2. \end{align} \tag{19}\]

The steady state concentration of Y₂ is (you can verify on your own if you like)

\[ \begin{align} y_2 = y_\mathrm{tot}\,\frac{x_\mathrm{tot}^2 s^4}{(1+s+s^2)^2 + x_\mathrm{tot} s^2 (1 + s + s^2) + x_\mathrm{tot}^2 s^4}. \end{align} \tag{20}\]

Analysis of the cascading circuits

We have considered four circuits, which we will label with number codes.

    1. The input signal catalyzes phosphorylation of one site of a signaling molecule.
    1. The input signal catalyzes phosphorylation of two sites of a signaling molecule.
  • (1, 1) The input signal catalyzes phosphorylation of one site of a signaling molecule, which in turn phorsphorylates one site of a subsequent signaling molecule.
  • (2, 2) The input signal catalyzes phosphorylation of two sites of a signaling molecule, which in turn phorsphorylates two sites of a subsequent signaling molecule.

For each circuit, we have written down the dynamical equations and worked out the steady state transfer function. As a convenient summary, the dimensionless transfer functions are:

\[ \begin{align} &(1): x_1 = x_\mathrm{tot}\,\frac{s}{1 + s},\\[1em] &(2): x_2 = x_\mathrm{tot}\,\frac{s^2}{1 + s + s^2},\\[1em] &(1, 1): y_1 = y_\mathrm{tot}\,\frac{x_\mathrm{tot}\,s}{1 + s + x_\mathrm{tot}\,s}, \\[1em] &(2, 2): y_2 = y_\mathrm{tot}\,\frac{x_\mathrm{tot}^2 s^4}{(1+s+s^2)^2 + x_\mathrm{tot} s^2 (1 + s + s^2) + x_\mathrm{tot}^2 s^4}. \end{align} \tag{21}\]

To gain insight to the nature of these signaling circuits, we wish to compute the quantities we defined in Chapter 7.

  1. The transfer function,
  2. The speed of the circuit response to a jump in input signal \(s\),
  3. The gain of the circuit,
  4. The derivative of the transfer function,
  5. The sensitivity.

In order assess the speed of the response of the circuit, we need to solve the dynamical equations numerically, which we will do momentarily. For the other quantities, we need the derivative of the transfer functions. They are listed below.

\[ \begin{align} &(1): \frac{\mathrm{d}x_1}{\mathrm{d}s} = x_\mathrm{tot}\,\frac{1}{(1 + s)^2},\\[1em] &(2): \frac{\mathrm{d}x_2}{\mathrm{d}s} = x_\mathrm{tot}\,\frac{s(2+s)}{\left(1 + s + s^2\right)^2},\\[1em] &(1, 1): \frac{\mathrm{d}y_1}{\mathrm{d}s} = y_\mathrm{tot}\,\frac{x_\mathrm{tot}}{\left(1 + s + x_\mathrm{tot}\,s\right)^2}, \\[1em] &(2, 2): \frac{\mathrm{d}y_2}{\mathrm{d}s} = y_\mathrm{tot}\,\frac{x_\mathrm{tot}^2 s^3(2+s)\left(2 + 2 s + (2+x_\mathrm{tot})s^2\right)}{\left((1+s+s^2)^2 + x_\mathrm{tot} s^2 (1 + s + s^2) + x_\mathrm{tot}^2 s^4\right)^2}. \end{align} \tag{22}\]

As an example, for circuit (1), the sensitivity is computed from these derivatives as

\[ \begin{align} \text{sensitivity} = \frac{\mathrm{d}\ln x_1}{\mathrm{d} \ln s} = \frac{s}{x_1}\,\frac{\mathrm{d}x_1}{\mathrm{d}s}. \end{align} \tag{23}\]

The gain is easily calculated for this example as \(x_1 / s\).

def f(s, n):
    if n == np.inf:
        output = np.zeros_like(s)
        output[s > 1] = 1 - 1 / s[s > 1]
    else:
        denom = 1
        for i in range(1, n+1):
            denom += s**i
        output = s**n / denom

    return output

def f1(s):
    return s / (1 + s)

def f2(s):
    return s ** 2 / (1 + s * (1 + s))

def f11(s, xtot):
    return f1(xtot * f1(s))

def f22(s, xtot):
    return f2(xtot * f2(s))

def f1_prime(s):
    return 1 / (1 + s) ** 2

def f2_prime(s):
    return s * (2 + s) / (1 + s * (1 + s)) ** 2

def f11_prime(s, xtot):
    return xtot * f1_prime(s) * f1_prime(xtot * f1(s))

def f22_prime(s, xtot):
    return xtot * f2_prime(s) * f2_prime(xtot * f2(s))

def deriv(s, xtot, deriv_fun):
    args = (s,) if deriv_fun in (f1_prime, f2_prime) else (s, xtot)
    return deriv_fun(*args)

def gain(s, xtot, transfer_fun):
    args = (s,) if transfer_fun in (f1, f2) else (s, xtot)
    return transfer_fun(*args) / s

def sensitivity(s, xtot, transfer_fun, deriv_fun):
    args = (s,) if transfer_fun in (f1, f2) else (s, xtot)
    return s / transfer_fun(*args) * deriv_fun(*args)
s = np.logspace(-2, 2, 400)

# Styling for figures
fig_kwargs = dict(
    frame_width=400,
    frame_height=200,
    x_axis_label="s",
    x_axis_type="log",
    x_range=[s.min(), s.max()],
)
colors = bokeh.palettes.Paired6[:2] + bokeh.palettes.Paired6[4:6]

# xtot slider
xtot_slider = bokeh.models.Slider(
    title=r"$$x_\mathrm{tot}$$",
    start=0.1,
    end=100,
    value=1,
    step=0.1,
)

data = dict(
    s=s,
    tf_1=f1(s),
    tf_2=f2(s),
    tf_11=f11(s, xtot_slider.value),
    tf_22=f22(s, xtot_slider.value),
    deriv_1=f1_prime(s),
    deriv_2=f2_prime(s),
    deriv_11=f11_prime(s, xtot_slider.value),
    deriv_22=f22_prime(s, xtot_slider.value),
    gain_1=f1(s) / s,
    gain_2=f2(s) / s,
    gain_11=f11(s, xtot_slider.value) / s,
    gain_22=f22(s, xtot_slider.value) / s,
    sens_1=f1_prime(s) * s / f1(s),
    sens_2=f2_prime(s) * s / f2(s),
    sens_11=f11_prime(s, xtot_slider.value) * s / f11(s, xtot_slider.value),
    sens_22=f22_prime(s, xtot_slider.value) * s / f22(s, xtot_slider.value),
)
cds = bokeh.models.ColumnDataSource(data)
titles = dict(
    tf="transfer function", gain="gain", deriv="derivative", sens="sensitivity"
)
y_axis_labels = dict(
    tf="fraction phosphorylated", gain="gain", deriv="derivative", sens="sensitivity"
)

plots = {
    feature: bokeh.plotting.figure(
        **fig_kwargs, y_axis_label=y_axis_labels[feature], title=titles[feature]
    )
    for feature in ["tf", "gain", "deriv", "sens"]
}

log_plots = {
    feature: bokeh.plotting.figure(
        **fig_kwargs,
        y_axis_type="log",
        y_axis_label=y_axis_labels[feature],
        title=titles[feature],
        visible=False
    )
    for feature in ["tf", "gain", "deriv", "sens"]
}

for feature in ["tf", "gain", "deriv", "sens"]:
    for color, cascade_code in zip(colors, ["1", "2", "11", "22"]):
        kwargs = dict(
            source=cds, x="s", y=feature + "_" + cascade_code, color=color, line_width=2
        )
        if feature == "tf":
            kwargs["legend_label"] = cascade_code
        plots[feature].line(**kwargs)
        log_plots[feature].line(**kwargs)
    if feature == "tf":
        plots[feature].legend.location = "bottom_right"
        log_plots[feature].legend.location = "bottom_right"

# Select log or linear y-scale
radio_button_group = bokeh.models.RadioButtonGroup(
    labels=["log", "linear"], active=1, width=30, orientation="vertical"
)
radio_button_group.js_on_change(
    "active",
    bokeh.models.CustomJS(
        args=dict(**{'p_'+feature: plots[feature] for feature in ["tf", "gain", "deriv", "sens"]},
                  **{'p_'+feature+'_log': log_plots[feature] for feature in ["tf", "gain", "deriv", "sens"]}),
        code="""
  if (p_tf_log.visible == true) {
    p_tf_log.visible = false;
    p_tf.visible = true;
    p_gain_log.visible = false;
    p_gain.visible = true;
    p_deriv_log.visible = false;
    p_deriv.visible = true;
    p_sens_log.visible = false;
    p_sens.visible = true;
  }
  else {
    p_tf_log.visible = true;
    p_tf.visible = false;
    p_gain_log.visible = true;
    p_gain.visible = false;
    p_deriv_log.visible = true;
    p_deriv.visible = false;
    p_sens_log.visible = true;
    p_sens.visible = false;
  }
""",
    ),
)

def callback(attr, old, new):
    data = dict(
        s=s,
        tf_1=f1(s),
        tf_2=f2(s),
        tf_11=f11(s, xtot_slider.value),
        tf_22=f22(s, xtot_slider.value),
        deriv_1=f1_prime(s),
        deriv_2=f2_prime(s),
        deriv_11=f11_prime(s, xtot_slider.value),
        deriv_22=f22_prime(s, xtot_slider.value),
        gain_1=f1(s) / s,
        gain_2=f2(s) / s,
        gain_11=f11(s, xtot_slider.value) / s,
        gain_22=f22(s, xtot_slider.value) / s,
        sens_1=f1_prime(s) * s / f1(s),
        sens_2=f2_prime(s) * s / f2(s),
        sens_11=f11_prime(s, xtot_slider.value) * s / f11(s, xtot_slider.value),
        sens_22=f22_prime(s, xtot_slider.value) * s / f22(s, xtot_slider.value),
    )
    cds.data = data
    
xtot_slider.on_change('value', callback)

        
layout = bokeh.layouts.column(
    bokeh.layouts.row(bokeh.layouts.Spacer(width=20), radio_button_group, bokeh.layouts.Spacer(width=20), xtot_slider),
    bokeh.layouts.Spacer(height=15),
    bokeh.layouts.row(
        bokeh.layouts.column(plots["tf"], log_plots["tf"]),
        bokeh.layouts.column(plots["gain"], log_plots["gain"]),
    ),
    bokeh.layouts.row(
        bokeh.layouts.column(plots["deriv"], log_plots["deriv"]),
        bokeh.layouts.column(plots["sens"], log_plots["sens"]),
    ),
)

def app(doc):
    doc.add_root(layout)
    
bokeh.io.show(app)

Now we can make a plot for various values of \(n\).

# Values of s and n for which we want a plot
s = np.logspace(-2, 2, 400)
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=30, 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)

Cascade of two multi-phosphorylation reactions

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.

\[\begin{align} \require{mhchem} &\ce{X_0 <=>[\mathrm{S}] X_1 <=>[\mathrm{S}] X_2 $\cdots$ <=>[\mathrm{S}] X$_\mathrm{n_x}$},\\[1em] &\ce{Y_0 <=>[\mathrm{X_{n_x}}] Y_1 <=>[\mathrm{X_{n_x}}] Y_2 $\cdots$ <=>[\mathrm{X_{n_x}}] Y$_\mathrm{n_y}$}. \end{align}\]

For simplicity, we will take \(n_x = n_y\), meaning that both X and Y have the same number of phosphorylation sites. We can use the same arguments as above to write down the dynamical equations for the concentration of species in this set of reactions.

\[\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},\\[1em] &\frac{\mathrm{d}y_0}{\mathrm{d}t} = - k_+ x_{n} y_0 + k_- y_1,\\[1em] &\frac{\mathrm{d}y_i}{\mathrm{d}t} = - k_+ x_{n} y_i + k_- y_{i+1} + k_+x_{n} y_{i-1} - k_- y_i\;\;\text{for } 1 \le i \le n-1,\\[1em] &\frac{\mathrm{d}y_{n}}{\mathrm{d}t} = k_+ x_{n} y_{n-1} - k_- y_{n}. \end{align}\]

We have again for simplicity taken all of the forward rate constants to be \(k_+\) and all of the reverse rate constants to be \(k_-\). At this point, it is worth nondimensionalizing the equations, taking \(t \leftarrow k_-t\), \(x_i \leftarrow x_i/K\), \(y_i \leftarrow y_i/K\), and \(s \leftarrow y_i/K\), giving

\[\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,\\[1em] &\frac{\mathrm{d}y_0}{\mathrm{d}t} = - x_n y_0 + y_1,\\[1em] &\frac{\mathrm{d}y_i}{\mathrm{d}t} = - x_n y_i + y_{i+1} + x_n y_{i-1} - y_i\;\;\text{for } 1 \le i \le n-1,\\[1em] &\frac{\mathrm{d}y_n}{\mathrm{d}t} = x_n y_{n-1} - y_n. \end{align}\]

Going forward we will use dimensionless variables, and note that \(x_\mathrm{tot}\) and \(y_\mathrm{tot}\) are also dimensionless, given by \(x_\mathrm{tot} \leftarrow x_\mathrm{tot}/K\) and \(x_\mathrm{tot} \leftarrow x_\mathrm{tot}/K\).

In looking at the above equations, we see that the equations for Y dynamics are analogous to those of X, except with \(x_n\) replacing \(s\). We already worked out that the steady state for \(x_{n}\) is \(x_{n} = x_\mathrm{tot} f_{n}(s)\). Analogously, the steady state for \(y_{n}\) is \(y_{n} = y_\mathrm{tot} f_{n}(x_{n})\). We can thus compactly and conveniently write the dependence of \(y_{n}\) on the input signal \(s\) as

\[\begin{align} y_{n} = y_\mathrm{tot}f_{n}(x_\mathrm{tot} f_{n}(s)). \end{align}\]

Since \(y_n\) is the output of interest of the circuit, the above equation constitutes the transfer function for this signaling circuit.

Analysis of the two-level-cascade

We have written down the dynamical equations and worked out the steady state transfer function for the two-level signaling cascade. To gain insight to the nature of this signaling circuit, we wish to quantify

  1. The speed of the circuit response to a jump in input signal \(s\),
  2. The gain of the circuit,
  3. The derivative of the transfer function,
  4. The sensitivity.

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'(z) = \frac{\mathrm{d}f_n}{\mathrm{d}z} = \left\{\begin{array}{ll} \frac{n}{2(1+n)} & \text{for }z = 1,\\[1em] \frac{z^{n-1}}{\left(1+z^{n+1}\right)^2}\left(n(1-z)-z\left(1-z^n\right)\right) & \text{otherwise}. \end{array} \right. \end{align}\]

For \(n\to\infty\), this is

\[\begin{align} f_\infty'(z) = \frac{\mathrm{d}f_\infty}{\mathrm{d}z} = \left\{ \begin{array}{ll} 0 & \text{for } z < 1,\\[1em] 1/2 & \text{for } z = 1,\\[1em] 1/z^2 & \text{for } z > 1,\\[1em] \end{array} \right. \end{align}\]

with \(f'_\infty(0) \to \infty\). With this definition, we can differentiate the transfer function as

\[\begin{align} \frac{\mathrm{d} y_n}{\mathrm{d}s} = x_\mathrm{tot}\,y_\mathrm{tot}\,f_n'(s)\,f_n'(x_\mathrm{tot} f(s)). \end{align}\]

The sensitivity is then

\[\begin{align} \text{sensitivity} = \frac{\mathrm{d}\ln y_n}{\mathrm{d} \ln s} = \frac{s}{y_n}\,\frac{\mathrm{d}y_n}{\mathrm{d}s} = \frac{s\,x_\mathrm{tot}\,f_n'(s)\,f_n'(x_\mathrm{tot} f_n(s))}{f_{n}(x_\mathrm{tot} f_{n}(s))}. \end{align}\]

Finally, the gain is easily calculated as \(y_n / s\).

def f_prime(x, n):
    if np.isinf(n):
        res = f_infinite_n(x)
    else:
        res = np.empty_like(x)

        # Indices where x is not one
        inds = ~np.isclose(x, 1)

        # Compute result for close to one and otherwise
        res[inds] = x[inds] ** n * (1 - x[inds]) / (1 - x[inds] ** (n + 1))
        res[~inds] = x[~inds] ** n / (1 + n)

    return res

\[\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_x-1,\\[1em] &\frac{\mathrm{d}x_{n_x}}{\mathrm{d}t} = k_+ s x_{n_x-1} - k_- x_{n_x},\\[1em] &\frac{\mathrm{d}y_0}{\mathrm{d}t} = - k_+ x_{n_x} y_0 + k_- y_1,\\[1em] &\frac{\mathrm{d}y_i}{\mathrm{d}t} = - k_+ x_{n_x} y_i + k_- y_{i+1} + k_+x_{n_x} y_{i-1} - k_- y_i\;\;\text{for } 1 \le i \le n-1,\\[1em] &\frac{\mathrm{d}y_{n_y}}{\mathrm{d}t} = k_+ x_{n_x} y_{n_y-1} - k_- y_{n_y}. \end{align}\]

We have again for simplicity taken all of the forward rate constants to be \(k_+\) and all of the reverse rate constants to be \(k_-\). At this point, it is work nondimensionalizing the equations, taking \(t \leftarrow k_-t\), \(x_i \leftarrow x_i/K\), \(y_i \leftarrow y_i/K\), and \(s \leftarrow y_i/K\), giving

\[\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,\\[1em] &\frac{\mathrm{d}y_0}{\mathrm{d}t} = - x_n y_0 + y_1,\\[1em] &\frac{\mathrm{d}y_i}{\mathrm{d}t} = - x_n y_i + y_{i+1} + x_n y_{i-1} - y_i\;\;\text{for } 1 \le i \le n_y-1,\\[1em] &\frac{\mathrm{d}y_n}{\mathrm{d}t} = x_n y_{n-1} - y_n. \end{align}\]

Going forward we will use dimensionless variables, and note that \(x_\mathrm{tot}\) and \(y_\mathrm{tot}\) are also dimensionless, given by \(x_\mathrm{tot} \leftarrow x_\mathrm{tot}/K\) and \(x_\mathrm{tot} \leftarrow x_\mathrm{tot}/K\).

In looking at the above equations, we see that the equations for Y dynamics are analogous to those of X, except with \(x_n\) replacing \(s\). We already worked out that the steady state for \(x_{n_x}\) is \(x_{n_x} = x_\mathrm{tot} f_{n_x}(s)\). Analogously, the steady state for \(y_{n_y}\) is \(y_{n_y} = y_\mathrm{tot} f_{n_y}(x_{n_x})\). We can thus compactly and conveniently write the dependence of \(y_{n_y}\) on the input signal \(s\) as

\[\begin{align} y_{n_y} = y_\mathrm{tot}f_{n_y}(x_\mathrm{tot} f_{n_x}(s)). \end{align}\]

Though not quite a function composition (the presence of the \(x_\mathrm{tot}\) term means it is not a composition), this form is also very convenient for differentiation since the chain rule may still be directly applied. Defining

\[\begin{align} f_n'(z) = \frac{\mathrm{d}f_n}{\mathrm{d}z} = \left\{\begin{array}{ll} \frac{n}{2(1+n)} & \text{for }z = 1,\\[1em] \frac{z^{n-1}}{\left(1+z^{n+1}\right)^2}\left(n(1-z)-z\left(1-z^n\right)\right) & \text{otherwise}, \end{array} \right. \end{align}\]

we have

\[\begin{align} \frac{\mathrm{d} y_n}{\mathrm{d}s} = x_\mathrm{tot}\,y_\mathrm{tot}\,f'(s)\,f'(y_\mathrm{tot} f(s)). \end{align}\]

Toward a goal of plotting the response \(y_n\) as a function of input signal and also the sensitivity and gain, we can code up expressions for \(y_n\) and \(\mathrm{d}y_n/\mathrm{d}s\). We can do it for more general expressions.

def cascade_signal(s, x_tot, n):
    if np.isscalar(x_tot):
        return x_tot * f(s, n)

    if np.isscalar(n):
        n = [n] * len(x_tot)
    
    res = x_tot[0] * f(s, n[0])
    for x_tot_val, n_val in zip(x_tot[1:], n[1:]):
        res = x_tot_val * f(res, n_val)

    return res


def cascade_signal_deriv(s, x_tot, n):
    if np.isscalar(x_tot):
        return x_tot * f(s, n)

    if np.isscalar(n):
        n = [n] * len(x_tot)
    
    res = x_tot[0] * f(s, n[0])
    for x_tot_val, n_val in zip(x_tot[1:], n[1:]):
        res = x_tot_val * f(res, n_val)

    return res
x_tot = 2
y_tot = 2
s = np.logspace(-2, 2, 400)

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()],
)

p = bokeh.plotting.figure(**fig_kwargs)

p_log = bokeh.plotting.figure(**fig_kwargs, y_axis_type='log', y_range=[1e-8, 10], visible=False)

radio_button_group = bokeh.models.RadioButtonGroup(
        labels=["log", "linear"], active=1, width=30, orientation="vertical"
    )
col = bokeh.layouts.column(
    p_log, p, bokeh.layouts.row(bokeh.models.Spacer(width=100), radio_button_group)
)

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;
  }
""",
    )
)

x_tot_slider = bokeh.models.Slider(start=0.1, end=10, step=0.1, value=x_tot)
y_tot_slider = bokeh.models.Slider(start=0.1, end=10, step=0.1, value=y_tot)

data = dict(s=s)
for n in [1, 2, 3, 4, np.inf]:
    data[str(n)] = cascade_signal(s, [x_tot, y_tot], n)

cds = bokeh.models.ColumnDataSource(data=data)

for n, color in zip([1, 2, 3, 4, np.inf], list(bokeh.palettes.Blues7[:4][::-1]) + ['black']):
    label = "n → ∞" if np.isinf(n) else f"n = {n}"
    p.line(source=cds, x='s', y=str(n), line_width=2, color=color, legend_label=label)
    p_log.line(source=cds, x='s', y=str(n), line_width=2, color=color, legend_label=label)

p.legend.location = 'bottom_right'
p_log.legend.location = 'bottom_right'

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)
bokeh.io.show(biocircuits.jsplots.phosphorylation_signal_cascade())

Computing environment

%load_ext watermark
%watermark -v -p numpy,scipy,bokeh,biocircuits,jupyterlab