Numerical simulatin of β-catenin destruction cycle¶
(c) 2020 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.plotting
import panel as pn
pn.extension()
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=450,
plot_height=350,
x_axis_label="time (min)",
y_axis_label="conc (nM)",
x_range=[0, 20],
)
# 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_label=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 plots with Panel¶
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.
Panel really facilitates this. It allows specification of interactive sliders that can then be used to change appearance of plots.
First, we will set up sliders for our six parameters, $k_8$, $k_{-8}$, $k_9$, $k_{10}$, $k_{11}$, and $k_{12}$. These allow us to change the values of the parameters.
k8_slider = pn.widgets.FloatSlider(
name="k8", start=0.001, end=0.1, step=0.0001, value=k8
)
km8_slider = pn.widgets.FloatSlider(
name="km8", start=0.001, end=100, step=0.001, value=km8
)
k9_slider = pn.widgets.FloatSlider(
name="k9", start=10.0, end=1000, step=0.1, value=k9
)
k10_slider = pn.widgets.FloatSlider(
name="k10", start=10.0, end=1000, step=0.1, value=k10
)
k11_slider = pn.widgets.FloatSlider(
name="k11", start=0.001, end=1, step=0.001, value=k11
)
k12_slider = pn.widgets.FloatSlider(
name="k12", start=10, end=1000, step=0.1, value=k12
)
Next, all we have to do is write a function to make a plot of the solution for a given set of parameters. This function essentially just contains code we already wrote. To link the parameters of the function to the values of the slider, we use the @pn.depends decorator from Panel.
@pn.depends(
k8_slider.param.value,
km8_slider.param.value,
k9_slider.param.value,
k10_slider.param.value,
k11_slider.param.value,
k12_slider.param.value,
)
def plot_solution(k8, km8, k9, k10, k11, k12):
# Unpack parameters
params = (k8, km8, k9, k10, k11, k12)
# Run the solver
c = scipy.integrate.odeint(dcdt, c0, t, args=params)
# Set up canvas on which to paint plot
p = bokeh.plotting.figure(
plot_width=450,
plot_height=350,
x_axis_label="time (min)",
y_axis_label="conc (nM)",
x_range=[0, 20],
)
# 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_label=name)
# Place legend
p.legend.location = "center_right"
return p
Finally, we can set up a layout. We stack the slider widgets as a column and place them next to the plot. Note that this interactivity only works in a running notebook, not in an HTML rendering of the notebook. Note also that we have not included units for the parameters on the slider, but the units are indicated above.
widgets = pn.Column(
k8_slider, km8_slider, k9_slider, k10_slider, k11_slider, k12_slider, width=300
)
pn.Row(plot_solution, pn.Spacer(width=15), widgets)
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.
Computing environment¶
%load_ext watermark
%watermark -v -p numpy,scipy,bokeh,panel,jupyterlab