Numerical simulatin of β-catenin destruction cycle¶
(c) 2019 Justin Bois. With the exception of pasted graphics, where the source is noted, this work is licensed under a Creative Commons Attribution License CC-BY 4.0. All code contained herein is licensed under an MIT license.
This document was prepared at Caltech with financial support from the Donna and Benjamin M. Rosen Bioengineering Center.
This document was generated from a Jupyter notebook, which may be downloaded here.
import numpy as np
import pandas as pd
import scipy.integrate
import bokeh.io
import bokeh.models
import bokeh.plotting
import bokeh.application
import bokeh.application.handlers
bokeh.io.output_notebook()
In this document, I present a numerical solution for the dynamical equations describing the β-catenin destruction cycle, as described in the Lee, et al. paper. As a reminder, the destruction cycle is part of the broader Wnt signaling pathway, as shown in this figure from the Lee, et al. paper.
To aid in interpreting the dynamical equations, the chemical species and their corresponding numerical label that are involved in the cycle are:
| number | species |
|---|---|
| 3 | APC*/Axin*/GSK-3 |
| 8 | β-catenin/APC*/Axin*/GSK-3 |
| 9 | β-catenin*/APC*/Axin*/GSK-3 |
| 10 | β-catenin* |
| 11 | β-catenin |
The mass action equations describing the dynamics are
\begin{align} &\frac{\mathrm{d}c_3}{\mathrm{d}t} = -k_8 c_3 c_{11} + k_{\text{-}8}c_8 + k_{10}c_9, \label{eq:41} \\ &\frac{\mathrm{d}c_8}{\mathrm{d}t} = k_8 c_3 c_{11} - k_{\text{-}8}c_8 - k_9c_8, \label{eq:42}\\ &\frac{\mathrm{d}c_9}{\mathrm{d}t} = k_9c_8 - k_{10} c_9,\label{eq:43} \\ &\frac{\mathrm{d}c_{10}}{\mathrm{d}t} = k_{10} c_9 - k_{11}c_{10}, \label{eq:44}\\ &\frac{\mathrm{d}c_{11}}{\mathrm{d}t} = k_{12} - k_8 c_3 c_{11} + k_{\text{-}8}c_8. \end{align}
Most of the parameters are measured and reported, at least approximately, in the Lee, et al. paper. The parameters $k_8$, $k_{\text{-}8}$ and $k_{12}$ are unknown. However, $K_\mathrm{d} = k_{\text{-}8}/k_8$ is known, so we only have two parameters to set. We will develop an interactive plotting app to allow us to investigate the dynamics of the respective species as the two parameters $k_8$ and $k_{-8}$ are varied. In doing so, we will see how the scipy.integrate.solve_ivp() function works to numerically solve systems of ordinary differential equations.
Numerical integration using SciPy¶
There are three main APIs for solving real-valued initial value problems in the SciPy library. They are scipy.integrate.solve_ivp(), ode(), and odeint(). According to the SciPy developers, scipy.integrate.solve_ivp() is the preferred method, with the others labeled as having an "old" API. The solve_ivp() function has the flexibility of allowing choice of multiple numerical algorithms for solving ODEs. However, for the kinds of problems we encounter in this class, I find that the generic LSODA algorithm developed by Linda Petzold and Alan Hindmarsh that handles both stiff and non-stiff problems with variable time stepping is the best option. This is the only solver offered in the odeint() function. If we compare the two solvers, solve_ivp() and odeint(), the former has a large overhead, which can lead to performance issues for small problems (for large problems, this is not a big deal). Since most of our problems are small, we will use odeint(). It has much better performance, and though its API is different, it is still intuitive.
The basic call signature for scipy.integrate.odeint() is
scipy.integrate.odeint(func, y0, t, args=())
There are many other keyword arguments to set algorithmic parameters, but we will generally not need them (and you can read about them in the docs.
Importantly, fun is a vector-valued function with call signature f(y, t, *args) that specifies the right hand side of the system of ODEs to be solved. t is a scalar time point and y is a one-dimensional array (though multidimensional arrays are possible; see the SciPy docs). y0 is an array with the initial conditions.
Let's put together some code to solve the ODEs describing the β-catenin destruction cycle!
Numerical solution of the β-catenin destruction cycle dynamics¶
First, the preliminaries. We'll sed up parameters and names of the respective species.
# Key for names
names = ['Axin complex', 'Axin-βcat', 'Axin-βcat*', 'βcat*', 'βcat']
# Define known parameters
c_A = 50 # nM (given by fixed GSK-3 concentration)
k9 = 206 # 1/min
k10 = 206 # 1/min
k11 = 0.417 # 1/min
Kd8 = 120 # nM
# Unknown parameters
km8 = 1 # 1/min
k12 = 100 # nM/min
# k8 determined form Kd8 and km8
k8 = km8 / Kd8 # 1/nM-min
# Set up params for ODE
params = (k8, km8, k9, k10, k11, k12)
For our initial conditions, we will assume that we have no $\beta$-catenin and that we have a fixed amount of APC/Axin/GKC-3.
# Initial conditions
c0 = np.array([c_A, 0, 0, 0, 0])
Now we define the right hand side of the ODEs. We do need to pass in parameters. We list the parameters after the required arguments, t and y (for the present case, y is c). The function needs to return an array of time derivatives, in the present case five of them.
def dcdt(c, t, k8, km8, k9, k10, k11, k12):
"""
Time derivative of concentrations.
c = (c3, c8, c9, c10, c11)
"""
# Unpack concentrations
c3, c8, c9, c10, c11 = c
# Build derivatives
deriv = np.empty(5)
deriv[0] = -k8*c3*c11 + km8*c8 + k10*c9
deriv[1] = k8*c3*c11 - (km8 + k9)*c8
deriv[2] = k9*c8 - k10*c9
deriv[3] = k10*c9 - k11*c10
deriv[4] = -k8*c3*c11 + km8*c8 + k12
return deriv
Now, we just have to specify what time points we want.
# Set up time points
t = np.linspace(0, 20, 300)
We are finally ready to run the solver.
c = scipy.integrate.odeint(dcdt, c0, t, args=params)
The result is a Numpy array where each column corresponds to one of the variables we are solving for in the ODE, in this case the concentrations. We can use these to plot the solution.
# Set up canvas on which to paint plot
p = bokeh.plotting.figure(plot_width=650,
plot_height=400,
x_axis_label='time (min)',
y_axis_label='conc (nM)')
# Specify colors
colors = bokeh.palettes.d3['Category10'][5]
# Add glyphs
for c_res, name, color in zip(c.transpose(), names, colors):
p.line(t, c_res, line_width=3, color=color, legend=name)
# Place legend
p.legend.location = 'center_right'
bokeh.io.show(p)
The Axin-β-catenin complexes, both activated and not activated, remain at low levels and closely track each other (if you zoom, you can see that the orange and green curves do not quite perfectly overlap). The unbound Axin complex does not vary much from its initial concentration of 50 nM.
Interactive plotting¶
When investigating how parameter values affect dynamics, it is convenient to be able to quickly adjust parameter values and see how the response changes. Interactive plotting is a very useful tool in this regard. We can make sliders to vary parameter values and see how the dynamics change.
We will use Bokeh to generate an interactive plot. We write a function to make the plot, with the data being sourced with a Bokeh ColumnDataSource. We then create sliders and write a callback function that dictates how the data source changes with the slide values.
def plot_app(doc):
# Set up canvas on which to paint plot
p = bokeh.plotting.figure(plot_width=650,
plot_height=400,
x_axis_label='time (min)',
y_axis_label='conc (nM)')
# Specify colors
colors = bokeh.palettes.d3['Category10'][5]
# Set up a data source
df = pd.DataFrame(data=c, columns=names)
df['t'] = t
source = bokeh.models.ColumnDataSource(data=df)
# Add glyphs
for name, color in zip(names, colors):
p.line(x='t',
y=name,
source=source,
line_width=3,
color=color,
legend=bokeh.core.properties.value(name))
# Place legend
p.legend.location = 'top_left'
# Callback function to update curves with new solve on slider change
def callback(attr, old, new):
# Extract parameter values from sliders
km8 = 10**log10_km8_slider.value
k8 = km8 / Kd8
k12 = k12_slider.value
params = (k8, km8, k9, k10, k11, k12)
# Compute curves
c = scipy.integrate.odeint(dcdt, c0, t, args=params)
# Update data source
df = pd.DataFrame(c, columns=names)
df['t'] = t
source.data = dict(df)
# Set up sliders
log10_km8_slider = bokeh.models.Slider(
title='log₁₀ k₋₈ [log₁₀(1/min)]',
value=np.log10(km8),
start=-2.0,
end=4.0,
step=0.1)
log10_km8_slider.on_change('value', callback)
k12_slider = bokeh.models.Slider(
title='k₁₂ [1/nm-min]',
value=k12,
start=20.0,
end=100.0,
step=1.0)
k12_slider.on_change('value', callback)
# Build app
widgets = bokeh.layouts.widgetbox([log10_km8_slider, k12_slider])
doc.add_root(bokeh.layouts.column(p, widgets))
# Deploy (change notebook handle appropriately)
handler = bokeh.application.handlers.FunctionHandler(plot_app)
bokeh.io.show(bokeh.application.Application(handler),
notebook_handle='localhost:8888')
The interactive plotting shows that $k_{12}$ simply sets the scale of the concentration of β-catenin, with the axic complex concentrations being essentially unchanged. The parameter $k_{-8}$ adjusts the relative amounts of active and inactive β-catenin.