14. MultiFate enables expandable and controllable multistability
Design principles
Combinatorial dimerization and positive feedback generate multistability
Techniques
Phase portraits
Calculating separatrices
Key Papers
This notebook is still very much in draft form.
This draft notebook requires additional dependencies to run. In particular, you will need to install the package intersection in order to run the notebook fully. You can do this by opening a terminal window and typing
pip install intersect
and then reloading this notebook.
[1]:
# Colab setup ------------------
import os, sys, subprocess
if "google.colab" in sys.modules:
cmd = "pip install --upgrade biocircuits intersect watermark"
process = subprocess.Popen(cmd.split(), stdout=subprocess.PIPE, stderr=subprocess.PIPE)
stdout, stderr = process.communicate()
data_path = "https://biocircuits.github.io/chapters/data/"
else:
data_path = "data/"
# ------------------------------
import numpy as np
import pandas as pd
import scipy.integrate
import scipy.optimize
import biocircuits
import biocircuits.jsplots
import bokeh.io
import bokeh.layouts
import bokeh.models
import bokeh.plotting
import colorcet
from skimage import measure
# Third-party package for finding fixed points.
# https://github.com/sukhbinder/intersection
from intersect import intersection
# Set to True to have fully interactive plots with Python;
# Set to False to use pre-built JavaScript-based plots
interactive_python_plots = False
notebook_url = "localhost:8888"
from bokeh.util.warnings import BokehUserWarning
import warnings
warnings.simplefilter(action='ignore', category=BokehUserWarning)
bokeh.io.output_notebook()
Multicellular organisms rely on a repertoire of cell states
The cells in our bodies exist in thousands of different states, differing in morphology, functional capabilities, gene expression, and other features. This multitude of cell states is generated through the astonishing process of multicellular development. As the embryo develops, cells generally commit themselves and their descendants to ever more specific cell fates: initially they decide among the germ layers—endoderm, ectoderm, and mesoderm. Later on, they will further commit to more various tissue and organ-specific cell types.
The late 2010s and early 2020s have witnessed a revolution in our understanding of cell states. Single cell “omics” techniques like single-cell RNA sequencing (scRNAseq) have revealed the molecular states of individual cells, allowing the creation of “cell atlases” of entire adult organisms. Spatial transcriptomics approaches have further shown where cells in those states are located within complex spatially organized tissues. Furthermore, emerging molecular recording techniques such as GESTALT, MEMOIR, and others are beginning to reveal the lineage relationships of individual cells, revealing how cells reached their fates. As a result, we have begun to learn what states exist, where they are located, and how they arose.
In parallel with these developments, in order to treat disease, the fields of stem cell biology and regenerative medicine have identified procedures that direct stem cells to differentiate into desired cell fates. One strategy is to recapitulate in a cell culture dish the signals that cells are likely to receive in the embryo as they differentiate into the fate of interest. This involves stimulating or inhibiting fate-controlling signaling pathways such as Notch, Wnt, Sonic hedgehog, FGF, EGF, BMP (which we covered in a previous chapter) and others, in specific sequences. A sense of the power and breadth of these approaches can be gleaned from this wall chart, published by Cell in 2012. Since then, these approaches have continued to become ever more powerful.
As a result, we now know, more than ever before, what cell types are out there and how to obtain them. Nevertheless, it has remained challenging to understand the cell fate control circuits that makes all of these incredible behaviors possible. More specifically, what architectures do those circuits use, how do they determine the repertoire of cell state that can exist in an organism? How do they control which states can differentiate to which others? And so on. In short, how do circuits specify what states can exist and how they can be controlled?
Much of our thinking has been guided by an evocative and iconic illustration from Conrad Waddington’s 1957 book, The Strategy of the Genes.
This drawing provides a wonderful intuitive metaphor for the way in which a cell, figuratively represented by a marble, might progressively commit to ever more specific cell states, by rolling deeper and deeper into alternative valleys from which is cannot easily emerge, until it eventually comes to rest in a terminally differentiated state. In the right-hand illustration, Waddington further imagined how underlying genes might control the shape of this landscape, influencing which valleys the marble is more or less likely to fall into. He represented this idea abstractly with genes (black rectangles) attached to the landscape controlling the shape of the landscape via wires (black lines). However, while these drawings are intuitively appealing and conceptually useful, they do not by themselves explain how circuits of interacting molecuels could establish and control cell fate decision-making within an organism.
One way to address this questions is by identifying the genes that specify natural cell fates and working out the ways in which they interact with one another to control cell fates. Such approaches have provided critical insights into many cell fate control systems, such as the pluripotency system that maintains embryonic stem cells in a unique ‘pluripotent’ state with the potential to differentiate into almost any cell type in the organism. A complementary approach is to try to design and build a synthetic cell fate control system “from scratch” and to see what lessons that reveals about what types of circuits are sufficient to enable such complexity. In this chapter, we will explore this latter approach.
To see what we are up against in engineering a synthetic cell fate control system, consider some of the many capabilities it should, ideally, possess:
Multistability. The ideal system should be able to put cells into multiple distinct, and stable, cell states.
Progressive commitment. the system should also recapitulate the ubiquitous property of increasingly specific cell fate specification observed in so many contexts.
Reprogramming. In 2006 Takahashi and Yamanaka reported the ability to reprogram a fully differentiated cell back into an embryonic stem cell. This “rebooting” of the cell is also something a synthetic system should be able to achieve.
Expandability. During evolution, increases in organismal complexity coincided with concomitant increases in the number of different cell types in the organism. Natural cell fate control systems thus have the capacity to expand to accomodate more and more states. In the synthetic case, a similar property would enable one to expand the capacity of an engineered circuit to acheive ever more complex multicellular behaviors.
These features are illustrated here:
Can you ever imagine engineering a synthetic differentiation circuit that could rival the complexity of the complex multicellular organisms we see in nature? How would you design it?
Circuit architectures for cell differentiation
To begin to explore how these capabilities emerge from underlying circuits, we will shift from Waddington’s metaphorical diagram to a more directly applicable dynamical systems perspective.
We imagine each cell state as a stable attractor in a phase space of concentrations of transcription factors or other regulators. This representation is more general and does not require an overall “slope” down which the marble will roll.
Within this general framework, the ideal circuit should:
Allow external inputs to switch cells between states.
Enable progressive commitment and reprogramming.
Expand to generate more states with additional genes.
This last feature—expandability—means that one could scale up the circuit by adding to it, e.g. incorporating additional genes, without re-engineering existing circuitry. It is important if we want the system to be able to scale to produce complex tissues or organisms with many states.
These features are diagrammed schematically here:
The first multistable circuit that may come to mind is the classic two component toggle switch, which we analyzed in Chapter 3. In principle, this circuit might be extended by incorporating many transcription factors that all repress one another. However, this would be difficult to achieve. With \(N\) transcription factors, the promoter controlling each repressor gene would have to be repressed by \(N-1\) repressors. This could become impractical given the limited available space in a promoter for effective repression. More fundamentally, however, it does not expand easily: incrememting the system from \(N\) to \(N+1\) repressors would require adding binding sites for the new repressor to all of the \(N\) existing repressors.
Before reading further, take a moment to think about what kind of architecture you might use to create a cell fate control system providing the features above.
Inspiration: natural circuits use combinatorial dimerization and positive autoregulation
Inspiration for the design of a synthetic cell fate control system can come from looking at what is known about natural cell fate control systems. Zhu et al. (2022, Science) noticed that many natural fate control systems share two critical features:
First, they are based on sets of transcription factors that combinatorially dimerize to form many possible complexes. (This is reminescent of the combinatorial complex formation explored in the previous chapter on BMP.) Certain dimers activate downstream genes, whie others may be inactive. Formation of inactive complexes sequesters certain transcription factors, preventing them from forming alternative, active dimers with other proteins. This is analogous to the molecular titration motif we analyzed in Chapter 12.
Second, these circuits contain positive autoregulatory feedback loops, in which transcription factor dimers activate expression of one or more of their own constituents.
For example, in the circuit controlling muscle cell differentiation, the master transcription factor MyoD can dimerize with Id or E47. The MyoD-Id dimer is inactive. But te MyoD-E47 dimer transcriptionall activates many genes involved in myogenesis, including MyoD itself. Similarly, embryonic stem cells can remain in a “pluripotent” state in which they have the potential to differentiate into virtually all cell types in the embryo. The circuit that controls this state includes the key transcription factor Oct4, which can dimerize with either Sox2 or Sox17. Sox2-Oct4 and Sox17-Oct4 respectively activate expression of Sox2 and Sox17, generating autoregulation.
MultiFate: a naturally inspired synthetic cell fate control circuit
Intrigued and inspired by these architectures, Zhu et al. desigend a circuit called MultiFate, which incorporates the two themes of combinatorial dimerizaton among a set of transcription factors and positive autoregulation. We will see that the MultiFate architecture can generate multistability and the other features enough, but is simple enough to be constructed and analyzed synthetically within cells.
We will start by analzying a minimal version of the MultiFate architecture, termed MultiFate-2, comprising just two transcription factors, which we label A and B. A and B can each form homodimers. They can also heterodimerize with one another. The homodimers positively activate their own gene, while heterodimers are inactive. The reason that homodimers can activate their own genes, while the corresponding monomers cannot is that the DNA binding site for the dimer is effectively double the length of that of the monomer, providing approximately twice the binding energy.
The schematic drawing below explains the structure of MultiFate-2.
A two-dimensional model of MultiFate dynamics
How many states can this architecture generate? Can it be expanded? And what other dynamic properties does it possess? To find out, we will apply the approaches developed previously with additional phase plane analysis.
To begin, we will define a few variables. The A and B transcription factors can each exist either in monomeric, homodimeric, or heterodimeric form. We will represent their concentrations by the variables \([A]\) (or \([B]\)), \([A_2]\) (or \([B_2]\)), and \([AB]\), respectively.
Transcription factors can interconvert between these forms via the binding reactions:
\begin{align} & A + A \rightleftharpoons A_2 \\ & A + B \rightleftharpoons AB \\ & B + B \rightleftharpoons B_2. \end{align}
All of the binding reactions, whether they lead to homodimers or heterodimers, occur through the same type of protein-protein interaction mediated by a single engineered homodimerization motif that is incorporated in both of the synthetic transcription factors. The homodimerization reactions will therefore share a single equilibrium binding constant \(K_d\). The factor of 2 that appears in the heterodimerization reaction is a correction factor to account for the fact that the same heterodimer can form in two different ways (A+B and B+A). It prevents double-counting.
When we set these binding reactions to equilibrium, we have
\begin{align} & [A_2] = [A]^2/K_d \\ & [B_2] = [B]^2/K_d \\ & [AB] = 2[A][B]/K_d. \end{align}
We will now apply separation of timescales, making the quasi-steady-state assumption that the binding reactions are approximately at equilibrium over the timescales over which the total concentrations, \([A_{\mathrm{tot}}] = [A] + [AB] + 2[A_2]\), change due to their production and removal. We can then write down a set of ODEs that describe the slower timescale dynamics, incorporating positive autoregulation by the homodimers:
\begin{align} & \frac{\mathrm{d}[A_{\mathrm{tot}}]}{\mathrm{d}t} = \alpha + \frac{\beta [A_2]^n}{K_M^n + [A_2]^n} - \delta [A_{\mathrm{tot}}] \\ & \frac{\mathrm{d}[B_{\mathrm{tot}}]}{\mathrm{d}t} = \alpha + \frac{\beta [B_2]^n}{K_M^n + [B_2]^n} - \delta [B_{\mathrm{tot}}] \end{align}
Finally, we note that binding and unbinding reactions conserve mass:
\begin{align} & [A_{\mathrm{tot}}] = [A] + 2 [A_2] + [AB] \\ & [B_{\mathrm{tot}}] = [B] + 2 [B_2] + [AB] \end{align}
Conservation of mass reduces the model to two dimensions
These equations make a few choices and assumptions. First, we chose to use \([A_{\mathrm{tot}}]\) and \([B_{\mathrm{tot}}]\) as the dynamical variables appearing in the \(d/dt\) terms on the left hand side. This decision will simplify the analysis below, and is reasonable because it is these “total” concentrations that vary over the slower timescsales we will be interested in. Second, we represented transcriptional regulation using a leaky Hill function regulated by the homodimeric form of the transcription factor, with leak term \(\alpha\) and half-activation constant \(K_M\). The inclusion of this basal transcription is important for allowing the system to self-activate and avoiding a potential stable state in which both factors are off. (Such a state would be difficult to “read out” but could be physiologically useful in some contexts.) Finally, we assumed that removal of the two proteins occurs at the same rate regardless of whether they are in a monomeric or dimeric form, leading to the removal terms \(\delta [A_{\mathrm{tot}}]\) and \(\delta [B_{\mathrm{tot}}]\).
How many steady states can this circuit generate? And how does the answer depend on its parameters? To find out, we will solve for the steady states by setting the time derivatives to zero. Our two ODEs contain four variables: \([A_{\mathrm{tot}}], [B_{\mathrm{tot}}], [A_2], [B_2]\). To make this solvable (2 equations with 2 unknowns) we need to eliminate two of these variables. Using the equilibrium binding relations we obtained above, alongside the conservation of mass expressions, we can obtain an expression that represents \([A_2]\) and \([B_2]\) as a function of \([A_{\mathrm{tot}}]\) and \([B_{\mathrm{tot}}]\).
We accomplish this by the following procedure:
Rewrite the equilibrium binding relations so that everything is a function of the homodimer concentrations, i.e. \begin{align}&[A] = \sqrt{K_d[A_2]} \\ &[AB] = 2[A][B]/K_d = 2\sqrt{[A_2][B_2]}\end{align}
Substitute these expressions into the conservation law to obtain an equation relating the total concentrations to the homodimer concentrations only.
Rearrange the terms and solve the resulting system of equations to obtain expressions for the homodimer concentrations as a function of just the total concentrations.
Following the above process, we obtain \begin{align} & [A_2] = \frac{2 [A_{\mathrm{tot}}]^2}{K_d + 4\left([A_{\mathrm{tot}}] + [B_{\mathrm{tot}}]\right) + \sqrt{K_d^2 + 8\left([A_{\mathrm{tot}}] + [B_{\mathrm{tot}}]\right)K_d}} \\ & [B_2] = \frac{2 [B_{\mathrm{tot}}]^2}{K_d + 4\left([A_{\mathrm{tot}}] + [B_{\mathrm{tot}}]\right) + \sqrt{K_d^2 + 8\left([A_{\mathrm{tot}}] + [B_{\mathrm{tot}}]\right)K_d}}. \end{align}
If we substitute these expression back into our ODEs, then we now have a closed-form system where the rate of change of \([A_{\mathrm{tot}}]\) and \([B_{\mathrm{tot}}]\) are written as functions of only themselves.
Nondimensionalizing further simplifies the MultiFate model
Since the above expressions contain many parameters, it will be useful to nondimensionalize the system. We can make the substitutions
\begin{align} \tilde{t} &\leftarrow t\delta \\ \tilde{[A_{tot}]} &\leftarrow [A_{tot}]/K_M \\ \tilde{[B_{tot}]} &\leftarrow [B_{tot}]/K_M \\ \tilde{[A_2]} &\leftarrow [A_2]/K_M \\ \tilde{[B_2]} &\leftarrow [B_2]/K_M \\ \tilde{\alpha} &\leftarrow \alpha/(K_M \delta) \\ \tilde{\beta} &\leftarrow \beta/(K_M \delta) \\ \tilde{K_d} &\leftarrow K_d/K_M \end{align}
and then, dropping the tildes for notational convenience, arrive at a set of nondimensionalized equations:
\begin{align} \frac{\mathrm{d}[A_{tot}]}{\mathrm{d}t} &= \alpha + \frac{\beta [A_2]^n}{1 + [A_2]^n} - [A_{tot}] \\ \frac{\mathrm{d}[B_{tot}]}{\mathrm{d}t} &= \alpha + \frac{\beta [B_2]^n}{1 + [B_2]^n} - [B_{tot}] \\ [A_2] &= \frac{2[A_{tot}]^2}{K_d + 4\left( [A_{tot}] + [B_{tot}]\right) + \sqrt{K_d^2 + 8\left([A_{tot}] + [B_{tot}]\right)K_d}} \\ [B_2] &= \frac{2[B_{tot}]^2}{K_d + 4\left( [A_{tot}] + [B_{tot}]\right) + \sqrt{K_d^2 + 8\left([A_{tot}] + [B_{tot}]\right)K_d}} \end{align}
The MultiFate Phase portrait
Calculating nullclines
Plotting nullclines can be used to identify steady states and provide insight into the structure of a dynamical system. Even better, we can make the nullcline plot interactive, allowing us to see directly how different features depend on biochemical parameters. Recall from our analysis of the Toggle Switch in Chapter 3 that for a two-dimensional ODE system there will be two nullclines, one corresponding to all of the points in the plane where \(\mathrm{d}[A_{tot}]/\mathrm{d}t = 0\), and one corresponding to all of the points in the plane where \(\mathrm{d}[B_{tot}]/\mathrm{d}t = 0\). When the nullclines intersect, both ODEs are equal to zero and the system is at a fixed point.
The functional forms of the ODEs are complex, making it difficult to obtain a closed-form analytical expression for the nullclines. Instead, we will use numerical methods. Specifically, in order to find the \(\mathrm{d}[A_{tot}]/\mathrm{d}t\) nullcline we will calculate the value of \(\mathrm{d}[A_{tot}]/\mathrm{d}t\) across a grid of points in the plane, and then use the find_countours function from the scikit-measure package to find the points where
\(\mathrm{d}[A_{tot}]/\mathrm{d}t\) takes a particular value (in this case, zero). The same procedure can be used to find theh \(\mathrm{d}[B_{tot}]/\mathrm{d}t\) nullcline. This is way of doing it represents a kind of hack—its accuracy is limited—but it is convenient and can work in many situations.
The result can be seen in the interactive plot below. (Note that because this contour-finding method is not fast, the sliders on the plot are throttled so that the plot only updates when you let go of the mouse button as you move the slider.)
[2]:
class BSParams(object):
"""
Container for parameters for homodimeric bistable circuits.
Includes asymmetry parameters r, m, kappa, gamma.
"""
def __init__(self, **kwargs):
# Dimensionless parameters
self.alpha = 0.1*2 # leakiness
self.beta = 10*2 # Maximal production
self.Kd = 1 # Dissociation constant of dimers
self.n = 1.5 # Additional cooperativity of transcriptional activation
self.r = 1 # Leakiness ratio
self.m = 1 # A and B copy number ratio
self.kappa = 1 # Ratio of A and B's binding affinity with DNA
self.gamma = 1 # Ratio of stability between A and B
# Put in params that were specified in input
for entry in kwargs:
setattr(self, entry, kwargs[entry])
def df_dt(y, t, p):
"""
ODE of A and B based on parameters p.
Includes asymmetry parameters r, m, kappa, gamma
"""
Atot, Btot = y
A2 = 2*Atot**2 / (p.Kd + 4*(Atot+Btot) + np.sqrt(p.Kd**2 + 8 * (Atot+Btot) * p.Kd))
B2 = 2*Btot**2 / (p.Kd + 4*(Atot+Btot) + np.sqrt(p.Kd**2 + 8 * (Atot+Btot) * p.Kd))
return [p.alpha + p.beta * A2**p.n / (1**p.n + A2**p.n) - Atot,
p.alpha + p.beta * B2**p.n / (1**p.n + B2**p.n) - Btot]
def get_nullclines(params):
"""
Calculate nullclines for Atot and Btot based on parameters.
Includes asymmetry parameters r, m, kappa, gamma
"""
density = 1000
rangescale = 1.5
Amin = params.alpha*0.01
Amax = rangescale*params.beta
Bmin = params.r*params.alpha*0.01
Bmax = rangescale*params.m*params.beta
A_vals = np.linspace(Amin, Amax, density)
B_vals = np.linspace(Bmin, Bmax, density)
dAdt_vals = np.zeros((len(A_vals), len(B_vals)))
dBdt_vals = np.zeros((len(A_vals), len(B_vals)))
AA, BB = np.meshgrid(A_vals, B_vals)
# dummy time variable
t = np.linspace(0,100,10)
dAdt_vals, dBdt_vals = df_dt([AA, BB], t, params)
Anull_inds = measure.find_contours(dAdt_vals, 0)
Bnull_inds = measure.find_contours(dBdt_vals, 0)
assert len(Anull_inds) == 1
assert len(Bnull_inds) == 1
Anull_inds = Anull_inds[0]
Bnull_inds = Bnull_inds[0]
A_nullcline = np.empty((len(Anull_inds), 2))
A_nullcline[:,0] = A_vals[[int(Anull_inds[i,0]) for i in range(len(Anull_inds))]]
A_nullcline[:,1] = Amin + (Amax-Amin)*Anull_inds[:,1]/density
B_nullcline = np.empty((len(Bnull_inds), 2))
B_nullcline[:,0] = Bmin + (Bmax-Bmin)*Bnull_inds[:,0]/density
B_nullcline[:,1] = B_vals[[int(Bnull_inds[i,1]) for i in range(len(Bnull_inds))]]
return A_nullcline, B_nullcline
[3]:
# initial plotting
params = BSParams()
params.alpha = 0.2
params.beta = 15
params.Kd = 1
params.n = 1
params.r = 1
params.m = 1
params.kappa = 1
params.gamma = 1
A_nullcline, B_nullcline = get_nullclines(params)
cdsA = bokeh.models.ColumnDataSource(dict(Anull_Avals=A_nullcline[:,0], Anull_Bvals=A_nullcline[:,1]))
cdsB = bokeh.models.ColumnDataSource(dict(Bnull_Avals=B_nullcline[:,0], Bnull_Bvals=B_nullcline[:,1]))
# Build the plot
p = bokeh.plotting.figure(
frame_height=300,
frame_width=300,
x_axis_label="dimensionless A_tot",
y_axis_label="dimensionless B_tot",
)
p.line(source=cdsA, x="Anull_Avals", y="Anull_Bvals", line_width=2, legend_label='A_tot nullcline')
p.line(source=cdsB, x="Bnull_Avals", y="Bnull_Bvals", line_width=2, color='orange', legend_label='B_tot nullcline')
alpha_slider = bokeh.models.Slider(
title="alpha", start=0.1, end=0.5, step=0.1, value=0.2, width=150
)
beta_slider = bokeh.models.Slider(
title="beta", start=0.5, end=17, step=0.1, value=15, width=150
)
Kd_slider = bokeh.models.Slider(
title="Kd", start=0.5, end=10, step=0.1, value=1, width=150
)
n_slider = bokeh.models.Slider(
title="n", start=0.5, end=3, step=0.1, value=1, width=150
)
r_slider = bokeh.models.Slider(
title="r", start=0.1, end=10, step=0.1, value=1, width=150
)
m_slider = bokeh.models.Slider(
title="m", start=0.1, end=10, step=0.1, value=1, width=150
)
kappa_slider = bokeh.models.Slider(
title="kappa", start=0.1, end=10, step=0.1, value=1, width=150
)
gamma_slider = bokeh.models.Slider(
title="gamma", start=0.1, end=10, step=0.1, value=1, width=150
)
def callback(attr, old, new):
params = BSParams()
params.alpha = alpha_slider.value
params.beta = beta_slider.value
params.Kd = Kd_slider.value
params.n = n_slider.value
params.r = r_slider.value
params.m = m_slider.value
params.kappa = kappa_slider.value
params.gamma = gamma_slider.value
Anull, Bnull = get_nullclines(params)
cdsA.data["Anull_Avals"] = Anull[:,0]
cdsA.data["Anull_Bvals"] = Anull[:,1]
cdsB.data["Bnull_Avals"] = Bnull[:,0]
cdsB.data["Bnull_Bvals"] = Bnull[:,1]
alpha_slider.on_change("value_throttled", callback)
beta_slider.on_change("value_throttled", callback)
Kd_slider.on_change("value_throttled", callback)
n_slider.on_change("value_throttled", callback)
r_slider.on_change("value_throttled", callback)
m_slider.on_change("value_throttled", callback)
kappa_slider.on_change("value_throttled", callback)
gamma_slider.on_change("value_throttled", callback)
layout = bokeh.layouts.row(
p, bokeh.models.Spacer(width=30), bokeh.layouts.column(alpha_slider, beta_slider, Kd_slider, n_slider,
# r_slider, m_slider, kappa_slider, gamma_slider
)
)
def app(doc):
doc.add_root(layout)
if interactive_python_plots:
bokeh.io.show(app, notebook_url=notebook_url)
else:
bokeh.io.show(p)
Determining the stability of the fixed points
With the default parameters in the plot, you can count five nullcline intersections, or fixed points. For a cell fate control system we are interested in those that are stable. To determine their stability we return to Linear Stability Analysis (see Chapter 9).
Recall that to determine stability we evaluate the Jacobian matrix of partial derivatives, \(J\), at the fixed point of interest \((x_0, y_0)\):
\begin{align} J(x_0, y_0) = \left.\begin{bmatrix}J_{0,0} & J_{0,1} \\ J_{1,0} & J_{1,1} \end{bmatrix} \right|_{(x_0, y_0)} = \left.\begin{bmatrix}\frac{\partial f(x,y)}{\partial x} & \frac{\partial f(x,y)}{\partial y} \\ \frac{\partial g(x,y)}{\partial x} & \frac{\partial g(x,y)}{\partial y} \end{bmatrix} \right|_{(x_0, y_0)} \end{align}
for \(f(x,y) = \mathrm{d}[A_{tot}]/\mathrm{d}t\) and \(g(x,y) = \mathrm{d}[B_{tot}]/\mathrm{d}t\), where we have substituted \([A_{tot}]\) and \([B_{tot}]\) with \(x\) and \(y\), respectively, for notational convenience.
In order to find the locations of these fixed points, we can once again turn to numerical approaches— in this case, we will use the intersect package developed by Sukhbinder Singh, which identifies the intersection points of two curves.
Once we find a given intersection \((x_0, y_0)\), we know it is a fixed point so we can proceed to evaluating the Jacobian matrix at that point. However, the functions \(f(x,y)\) and \(g(x,y)\) have complicated functional forms. Since we are using these derivatives in the context of an otherwise fully numerical analysis of the system, we can take advantage of a symbolic mathematical calculator to determine closed-form expressions for these derivatives:
\begin{align} &J_{0,0} = -1 - \beta \left(1 + \left(\frac{K_d + 4(x_0 + y_0) + Q}{2x_0^2}\right)^n\right)^{-2} n \left(\frac{K_d + 4(x_0 + y_0) + Q}{2x_0^2}\right)^{n-1} \left(\frac{x_0 \left(4 + \frac{4K_d}{Q}\right) - 2(K_d + 4(x_0 + y_0) + Q)}{2x_0^3} \right) \\ &J_{0,1} = - \beta \left(1 + \left(\frac{K_d + 4(x_0 + y_0) + Q}{2x_0^2}\right)^n\right)^{-2} n \left(\frac{K_d + 4(x_0 + y_0) + Q}{2x_0^2}\right)^{n-1}\left(\frac{4 + \frac{4K_d}{Q} }{2x_0^2} \right) \\ &J_{1,0} = - \beta \left(1 + \left( \frac{K_d + 4(x_0 + y_0) + Q}{2y_0^2}\right)^n\right)^{-2}n \left(\frac{K_d + 4(x_0 + y_0) + Q}{2y_0^2}\right)^{n-1}\left(\frac{4 + \frac{4K_d}{Q} }{2y_0^2} \right) \\ &J_{1,1} = -1 - \beta \left(1 + \left( \frac{K_d + 4(x_0 + y_0) + Q}{2y_0^2}\right)^n\right)^{-2}n \left(\frac{K_d + 4(x_0 + y_0) + Q}{2y_0^2}\right)^{n-1} \left( \frac{y_0 \left(4 + \frac{4K_d}{Q}\right) -2 (K_d + 4(x_0 + y_0) + Q)}{2y_0^3} \right), \end{align}
where
\begin{align} Q \equiv \sqrt{K_d^2 + 8(x_0 + y_0) + K_d}. \end{align}
Let’s now overlay an indicator of fixed point stability over the nullcline plots that we made.
[4]:
# Code to calcaulte the jacobian symbolically each time to confirm
# the above equations are correct
# def symbolic_jacobian(v_str, f_list):
# varss = sym.symbols(v_str)
# f = sym.sympify(f_list)
# J = sym.zeros(len(f),len(varss))
# for i, fi in enumerate(f):
# for j, s in enumerate(varss):
# J[i,j] = sym.diff(fi, s)
# return J
# def get_jacobian_func(p):
# v_str = 'Atot Btot'
# A2 = f'(2*Atot**2 / ({p.Kd} + 4*(Atot+Btot) + sqrt({p.Kd}**2 + 8 * (Atot+Btot) * {p.Kd})))'
# B2 = f'(2*Btot**2 / ({p.Kd} + 4*(Atot+Btot) + sqrt({p.Kd}**2 + 8 * (Atot+Btot) * {p.Kd})))'
# f_list = [f'{p.alpha} + {p.beta}*{A2}**{p.n} / (1**{p.n} + {A2}**{p.n}) - Atot',
# f'{p.alpha} + {p.beta}*{B2}**{p.n} / (1**{p.n} + {B2}**{p.n}) - Btot']
# J = symbolic_jacobian(v_str, f_list)
# return lambdify(sym.symbols(v_str), J, modules='numpy')
# def stability(x0, y0, jacobian_func):
# jmat = jacobian_func(x0, y0)
# eig = np.linalg.eigvals(jmat)
# print(eig)
# return (np.real(eig)<0).all()
# def get_fixed_points(Anull, Bnull, params):
# # get jacobian
# jacobian_func = get_jacobian_func(params)
# # find fixed points
# x, y = intersection(Anull[:,0], Anull[:,1], Bnull[:,0], Bnull[:,1])
# stable_fp = np.empty((0,2))
# unstable_fp = np.empty((0,2))
# # for each fixed point, use the surrounding 4 points to determine the stability
# for i in range(len(x)):
# print(x[i], y[i])
# if stability(x[i], y[i], jacobian_func)==1:
# stable_fp = np.vstack([stable_fp, [x[i],y[i]]])
# print('stable exissted')
# else:
# print('unstable existed')
# unstable_fp = np.vstack([unstable_fp, [x[i],y[i]]])
# return stable_fp, unstable_fp
####################
# Using the Jacobian function Ron calculated
def jacobian(x0, y0, p):
# Calculate the jacobian matrix
# Fixed point = (x0, y0)
# p = parameters
jmat = np.zeros([2,2])
sqtemp = np.sqrt(p.Kd**2 + 8 * (x0 + y0) * p.Kd)
jmat[0,0] = - 1 - p.beta * (1 + ((p.Kd + 4 * (x0 + y0) + sqtemp) / (2 * x0**2))**p.n)**(-2) * p.n * ((p.Kd + 4 * (x0 + y0) + sqtemp) / (2 * x0**2))**(p.n-1) * (x0 * (4 + 4 * p.Kd / sqtemp) - 2 * (p.Kd + 4 * (x0 + y0) + sqtemp)) / (2 * x0**3)
jmat[0,1] = - p.beta * (1 + ((p.Kd + 4 * (x0 + y0) + sqtemp) / (2 * x0**2))**p.n)**(-2) * p.n * ((p.Kd + 4 * (x0 + y0) + sqtemp) / (2 * x0**2))**(p.n-1) * (4 + 4 * p.Kd / sqtemp) / (2 * x0**2)
jmat[1,0] = - p.m * p.beta * (1 + (p.kappa * (p.Kd + 4 * (x0 + y0) + sqtemp) / (2 * y0**2))**p.n)**(-2) * p.n * (p.kappa * (p.Kd + 4 * (x0 + y0) + sqtemp) / (2 * y0**2))**(p.n-1) * p.kappa * (4 + 4 * p.Kd / sqtemp) / (2 * y0**2)
jmat[1,1] = - p.gamma - p.m * p.beta * (1 + (p.kappa * (p.Kd + 4 * (x0 + y0) + sqtemp) / (2 * y0**2))**p.n)**(-2) * p.n * (p.kappa * (p.Kd + 4 * (x0 + y0) + sqtemp) / (2 * y0**2))**(p.n-1) * (y0 * p.kappa * (4 + 4 * p.Kd / sqtemp) - 2 * p.kappa * (p.Kd + 4 * (x0 + y0) + sqtemp)) / (2 * y0**3)
return jmat
def stability(x0, y0, p):
"""
Returns True if stable
"""
# x0, y0 = fixed point coord
# Calculate the stability of fixed point based on eigenvalue of jacobian matrix
jmat = jacobian(x0, y0, p)
eig = np.linalg.eigvals(jmat)
return (np.real(eig)<0).all()
def get_fixed_points(Anull, Bnull, params):
# find fixed points
x, y = intersection(Anull[:,0], Anull[:,1], Bnull[:,0], Bnull[:,1])
stable_fp = np.empty((0,2))
unstable_fp = np.empty((0,2))
# for each fixed point, use the surrounding 4 points to determine the stability
for i in range(len(x)):
if stability(x[i], y[i], params)==1:
stable_fp = np.vstack([stable_fp, [x[i],y[i]]])
else:
unstable_fp = np.vstack([unstable_fp, [x[i],y[i]]])
return stable_fp, unstable_fp
[5]:
# initial plotting
params = BSParams()
params.alpha = 0.2
params.beta = 15
params.Kd = 1
params.n = 1
params.r = 1
params.m = 1
params.kappa = 1
params.gamma = 1
A_nullcline, B_nullcline = get_nullclines(params)
stable_fps, unstable_fps = get_fixed_points(A_nullcline, B_nullcline, params)
cdsA = bokeh.models.ColumnDataSource(dict(Anull_Avals=A_nullcline[:,0], Anull_Bvals=A_nullcline[:,1]))
cdsB = bokeh.models.ColumnDataSource(dict(Bnull_Avals=B_nullcline[:,0], Bnull_Bvals=B_nullcline[:,1]))
cds_stable_fps = bokeh.models.ColumnDataSource(dict(x=stable_fps[:,0], y=stable_fps[:,1]))
cds_unstable_fps = bokeh.models.ColumnDataSource(dict(x=unstable_fps[:,0], y=unstable_fps[:,1]))
# Build the plot
p = bokeh.plotting.figure(
frame_height=300,
frame_width=300,
x_axis_label="dimensionless A_tot",
y_axis_label="dimensionless B_tot",
)
p.line(source=cdsA, x="Anull_Avals", y="Anull_Bvals", line_width=2, legend_label='A_tot nullcline')
p.line(source=cdsB, x="Bnull_Avals", y="Bnull_Bvals", line_width=2, color='orange', legend_label='B_tot nullcline')
p.circle(
source=cds_stable_fps,
x="x",
y="y",
color="black",
size=10,
line_width=1
)
p.circle(
source=cds_unstable_fps,
x="x",
y="y",
color="white",
line_color="black",
size=10,
line_width=1
)
alpha_slider = bokeh.models.Slider(
title="alpha", start=0.1, end=0.6, step=0.1, value=0.2, width=150
)
beta_slider = bokeh.models.Slider(
title="beta", start=0.5, end=17, step=0.1, value=15, width=150
)
Kd_slider = bokeh.models.Slider(
title="Kd", start=0.5, end=10, step=0.1, value=1, width=150
)
n_slider = bokeh.models.Slider(
title="n", start=0.5, end=3, step=0.1, value=1, width=150
)
r_slider = bokeh.models.Slider(
title="r", start=0.1, end=10, step=0.1, value=1, width=150
)
m_slider = bokeh.models.Slider(
title="m", start=0.1, end=10, step=0.1, value=1, width=150
)
kappa_slider = bokeh.models.Slider(
title="kappa", start=0.1, end=10, step=0.1, value=1, width=150
)
gamma_slider = bokeh.models.Slider(
title="gamma", start=0.1, end=10, step=0.1, value=1, width=150
)
def callback(attr, old, new):
params = BSParams()
params.alpha = alpha_slider.value
params.beta = beta_slider.value
params.Kd = Kd_slider.value
params.n = n_slider.value
params.r = r_slider.value
params.m = m_slider.value
params.kappa = kappa_slider.value
params.gamma = gamma_slider.value
Anull, Bnull = get_nullclines(params)
cdsA.data["Anull_Avals"] = Anull[:,0]
cdsA.data["Anull_Bvals"] = Anull[:,1]
cdsB.data["Bnull_Avals"] = Bnull[:,0]
cdsB.data["Bnull_Bvals"] = Bnull[:,1]
stable_fps, unstable_fps = get_fixed_points(Anull, Bnull, params)
cds_stable_fps.data["x"] = stable_fps[:,0]
cds_stable_fps.data["y"] = stable_fps[:,1]
cds_unstable_fps.data["x"] = unstable_fps[:,0]
cds_unstable_fps.data["y"] = unstable_fps[:,1]
alpha_slider.on_change("value_throttled", callback)
beta_slider.on_change("value_throttled", callback)
Kd_slider.on_change("value_throttled", callback)
n_slider.on_change("value_throttled", callback)
r_slider.on_change("value_throttled", callback)
m_slider.on_change("value_throttled", callback)
kappa_slider.on_change("value_throttled", callback)
gamma_slider.on_change("value_throttled", callback)
layout = bokeh.layouts.row(
p, bokeh.models.Spacer(width=30), bokeh.layouts.column(alpha_slider, beta_slider, Kd_slider, n_slider,
# r_slider, m_slider, kappa_slider, gamma_slider
)
)
def app(doc):
doc.add_root(layout)
if interactive_python_plots:
bokeh.io.show(app, notebook_url=notebook_url)
else:
bokeh.io.show(p)
Plotting a flow field
To investigate the behavior near the fixed points, and indeed away from them as well, we can plot the phase portrait. A phase portrait is a plot displaying possible trajectories of a dynamical system in phase space, in this case on the \(A_\mathrm{tot}\)-\(B_\mathrm{tot}\) plane. Let’s generate a phase portrait using the biocircuits package (with the nullclines overlaid).
The biocircuits.phase_portrait() function requires explicit functions for \(\mathrm{d}[A_{tot}]/\mathrm{d}t\) and \(\mathrm{d}[B_{tot}]/\mathrm{d}t\).
[6]:
def dAtot_dt(Atot, Btot, alpha, beta, n, Kd):
A2 = (
2
* Atot ** 2
/ (Kd + 4 * (Atot + Btot) + np.sqrt(Kd ** 2 + 8 * (Atot + Btot) * Kd))
)
return alpha + beta * A2 ** n / (1 ** n + A2 ** n) - Atot
def dBtot_dt(Atot, Btot, alpha, beta, n, Kd):
B2 = (
2
* Btot ** 2
/ (Kd + 4 * (Atot + Btot) + np.sqrt(Kd ** 2 + 8 * (Atot + Btot) * Kd))
)
return alpha + beta * B2 ** n / (1 ** n + B2 ** n) - Btot
Next, we will build a plot with the fixed points and nullclines, as before.
[7]:
# initial plotting
alpha = 0.2
beta = 15
Kd = 1
n = 1
r = 1
m = 1
kappa = 1
gamma = 1
density = 1000
rangescale = 1.5
Amin = params.alpha*0.01
Amax = rangescale*params.beta
Bmin = params.r*params.alpha*0.01
Bmax = rangescale*params.m*params.beta
params = BSParams()
params.alpha = alpha
params.beta = beta
params.Kd = Kd
params.n = n
params.r = r
params.m = m
params.kappa = kappa
params.gamma = gamma
A_nullcline, B_nullcline = get_nullclines(params)
stable_fps, unstable_fps = get_fixed_points(A_nullcline, B_nullcline, params)
cdsA = bokeh.models.ColumnDataSource(dict(Anull_Avals=A_nullcline[:,0], Anull_Bvals=A_nullcline[:,1]))
cdsB = bokeh.models.ColumnDataSource(dict(Bnull_Avals=B_nullcline[:,0], Bnull_Bvals=B_nullcline[:,1]))
cds_stable_fps = bokeh.models.ColumnDataSource(dict(x=stable_fps[:,0], y=stable_fps[:,1]))
cds_unstable_fps = bokeh.models.ColumnDataSource(dict(x=unstable_fps[:,0], y=unstable_fps[:,1]))
p = bokeh.plotting.figure(
frame_height=300,
frame_width=300,
x_axis_label="dimensionless A_tot",
y_axis_label="dimensionless B_tot",
)
p.line(source=cdsA, x="Anull_Avals", y="Anull_Bvals", line_width=2, legend_label='A_tot nullcline')
p.line(source=cdsB, x="Bnull_Avals", y="Bnull_Bvals", line_width=2, color='orange', legend_label='B_tot nullcline')
p.circle(
source=cds_stable_fps,
x="x",
y="y",
color="black",
size=10,
line_width=1
)
p.circle(
source=cds_unstable_fps,
x="x",
y="y",
color="white",
line_color="black",
size=10,
line_width=1
);
Finally, we overlay the phase portrait.
[8]:
p = biocircuits.phase_portrait(
dAtot_dt,
dBtot_dt,
[Amin, Amax],
[Bmin, Bmax],
(alpha, beta, n, Kd),
(alpha, beta, n, Kd),
x_axis_label="dimensionless A_tot",
y_axis_label="dimensionless B_tot",
color="gray",
p=p,
)
bokeh.io.show(p)
Plotting separatrices
Notice how two trajectories can start very close to each other yet diverge and end up at different stable fixed points. The set of all starting points that eventually land on a given attractor (stable fixed point) is known as its basin of attraction. The boundary between two basins of attraction is called a separatrix.
We can calculaute the separatrices by making use of the fact that a separatrix is the one-dimensional stable manifold of a particular saddle fixed point— while trajectories on either side of a separatrix will converge to a stable fixed point, trajectories that start exactly on a separatrix will converge to a saddle point.
Recall from our exploration of Linear Stability Analysis in Chapter 9 that a saddle fixed point has one stable eigendirection and one unstable eigendirection. Near the fixed point, one eigenvector, \(\boldsymbol{v}\), and its inverse \(-\boldsymbol{v}\), point outwards from the fixed point. This eigenvector is parallel to the separatrix close to the fixed point. However, moving away from the fixed point, the eigenvector and the separatrix will begin to diverge. In order to follow the separatrix, therefore, we must start at the fixed point, take a very small stepp along the \(\boldsymbol{v}\) direction, and then follow the ODEs themselves by integrating backwards in time to find all of the trajectories that would eventually converge to the fixed point. We can then repeat this procedure in the \(-\boldsymbol{v}\) direction to find the other arm of the separatrix.
We demonstrate this procedure below.
[9]:
def df_dt_tuple(concs, t, alpha, beta, Kd, n, r, m, kappa, gamma):
Atot, Btot = concs
A2 = 2*Atot**2 / (Kd + 4*(Atot+Btot) + np.sqrt(Kd**2 + 8 * (Atot+Btot) * Kd))
B2 = 2*Btot**2 / (Kd + 4*(Atot+Btot) + np.sqrt(Kd**2 + 8 * (Atot+Btot) * Kd))
return [-1 * (alpha + beta * A2**n / (1**n + A2**n) - Atot),
-1 * (alpha + beta * B2**n / (1**n + B2**n) - Btot)]
def plot_separatrix(
x0,
y0,
A_range,
B_range,
params,
p,
t_max=100,
eps=1e-6,
color="#d95f02",
line_width=2,
):
"""
Plot separatrix on phase portrait.
Takes bokeh plot p as an input and returns p with
separatrix plot overlaid.
"""
# Negative time function to integrate to compute separatrix
def rhs(concs, t):
# Unpack variables
Atot, Btot = concs
# return df_dt_tuple(concs, t, alpha, beta, Kd, n, r, m, kappa, gamma)
# Stop integrating if we get the edge of where we want to integrate
if A_range[0] < Atot < A_range[1] and B_range[0] < Btot < B_range[1]:
return df_dt_tuple(concs, t, alpha, beta, Kd, n, r, m, kappa, gamma)
else:
return np.array([0, 0])
# Compute eigenvalues and eigenvectors
evals, evecs = np.linalg.eig(jacobian(x0, y0, params))
# Get eigenvector corresponding to attraction
evec = evecs[:, evals < 0].flatten()
# Parameters for building separatrix
t = np.linspace(0, t_max, 400)
# Put parameters in a tuple
args = (params.alpha, params.beta, params.Kd, params.n, params.r, params.m, params.kappa, params.gamma)
alpha = params.alpha
beta = params.beta
Kd = params.Kd
n = params.n
r = params.r
m = params.m
kappa = params.kappa
gamma = params.gamma
# Build upper right branch of separatrix
c0 = np.array([x0, y0]) + eps * evec
x_upper = scipy.integrate.odeint(rhs, c0, t)
# Build lower left branch of separatrix
c0 = np.array([x0, y0]) - eps * evec
x_lower = scipy.integrate.odeint(rhs, c0, t)
# Concatenate, reversing lower so points are sequential
sep_x = np.concatenate((x_lower[::-1, 0], x_upper[:, 0]))
sep_y = np.concatenate((x_lower[::-1, 1], x_upper[:, 1]))
# Plot
p.line(sep_x, sep_y, color=color, line_width=line_width,legend_label="seperatrix")
return p
[10]:
params = BSParams()
params.alpha = 0.2
params.beta = 15
params.Kd = 1
params.n = 1
params.r = 1
params.m = 1
params.kappa = 1
params.gamma = 1
Amin = params.alpha*0.01
Amax = rangescale*params.beta
Bmin = params.r*params.alpha*0.01
Bmax = rangescale*params.m*params.beta
A_nullcline, B_nullcline = get_nullclines(params)
stable_fps, unstable_fps = get_fixed_points(A_nullcline, B_nullcline, params)
cdsA = bokeh.models.ColumnDataSource(dict(Anull_Avals=A_nullcline[:,0], Anull_Bvals=A_nullcline[:,1]))
cdsB = bokeh.models.ColumnDataSource(dict(Bnull_Avals=B_nullcline[:,0], Bnull_Bvals=B_nullcline[:,1]))
cds_stable_fps = bokeh.models.ColumnDataSource(dict(x=stable_fps[:,0], y=stable_fps[:,1]))
cds_unstable_fps = bokeh.models.ColumnDataSource(dict(x=unstable_fps[:,0], y=unstable_fps[:,1]))
# Build the plot
p = biocircuits.phase_portrait(
dAtot_dt,
dBtot_dt,
[Amin, Amax],
[Bmin, Bmax],
(alpha, beta, n, Kd),
(alpha, beta, n, Kd),
x_axis_label="dimensionless A_tot",
y_axis_label="dimensionless B_tot",
color="gray",
plot_width=450,
)
p.line(source=cdsA, x="Anull_Avals", y="Anull_Bvals", line_width=2, legend_label='A_tot nullcline')
p.line(source=cdsB, x="Bnull_Avals", y="Bnull_Bvals", line_width=2, color='orange', legend_label='B_tot nullcline')
p.circle(
source=cds_stable_fps,
x="x",
y="y",
color="black",
size=10,
line_width=1
)
p.circle(
source=cds_unstable_fps,
x="x",
y="y",
color="white",
line_color="black",
size=10,
line_width=1
)
# Calculate separatrices
for ufp in unstable_fps:
p = plot_separatrix(ufp[0], ufp[1], [Amin, Amax],[Bmin, Bmax], params, p, eps=1e-1)
bokeh.io.show(p)