Some tools to analyze dynamical systems¶
(c) 2018 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 an Jupyter notebook. You can download the notebook here. Not all of the interactive parts of this notebook will work in a web browser; you need to run the notebook with a Python kernel.
# Our engines
import numpy as np
import scipy.integrate
# Plotting modules
import bokeh.layouts
import bokeh.plotting
import bokeh.io
bokeh.io.output_notebook()
import ipywidgets
Reaction-diffusion models for morphogen patterning¶
A morphogen is a chemical species that determines cell fate in a concentration-dependent manner. Morphogens are spatially distributed throughout a developing organism, thereby imparting spatial variation in cell types. They present us with a great opportunity to connect our genetic network dynamics to spatial variation.
Today, we will consider reaction-diffusion (R-D) models for mophogen dynamics. Consider a set of biochemical species, $\mathbf{M} = \{M_1, M_2, \ldots\}$. The set of concentrations of these biochemical species is $\mathbf{c} = \{c_1, c_2, \ldots\}$. We can write down a system of partial differential equations describing the dynamics of morphogen concentrations (which is more or less a statement of conservation of mass) as
\begin{align} \frac{\partial c_i}{\partial t} &= \frac{\partial}{\partial x}\,D_i\,\frac{\partial c_i}{\partial x} + r_i(\mathbf{c}), \end{align}
valid for all $i$. Here, $r_i(\mathbf{c})$ describes the net rate of production of species $M_i$, dependent on the concentrations of all other species. The diffusion coefficient $D_i$ parametrizes the diffusion of species $i$. For our purposes today, we are restricting ourselves to one-dimensional systems. We also have implicitly assumed that the diffusivity tensor is diagonal. Finally, for all of our analyses in this lecture, we will assume $x \in [0,L]$ with constant flux (usually no flux) boundary conditions at the ends of the domain. Recall that the flux of morphogen is given by Fick's first law,
\begin{align} j_i &= -D_i\,\frac{\partial c_i}{\partial x}.\\ &\phantom{blah} \end{align}
A gradient from a point source¶
We will begin by analyzing the simplest case of morphogen patterning. We have a single morphogen, $M$, which has a source at $x = 0$. We model the source as a constant flux, $\eta$ at $x = 0$.
Simple decay of the morphogen¶
We begin with the case where our morphogen only undergoes a single reaction, degradation which we would represent as $M \to \emptyset$, and has a constant diffusion coefficient. Our PDE is
\begin{align} \frac{\partial c}{\partial t} &= D \,\frac{\partial^2 c}{\partial x^2} - \gamma c, \\[1em] \left.\frac{\partial c}{\partial x}\right|_{x=0} &= -\frac{\eta}{D}\\[1em] \left.\frac{\partial c}{\partial x}\right|_{x=L} &= 0. \end{align}
We can readily solve for the unique steady state in this case as
\begin{align} c(x) = \frac{\eta}{\sqrt{D\gamma}}\,\frac{\mathrm{e}^{-x/\lambda} + \mathrm{e}^{(x-2L)/\lambda}}{1 - \mathrm{e}^{-2L/\lambda}}, \end{align}
where $\lambda \equiv \sqrt{D/\gamma}$. So, we have an exponentially decaying morphogen profile with decay length $\lambda$. Typically, $L \gg \lambda$, so we have
\begin{align} c(x) &\approx \frac{\eta}{\sqrt{D\gamma}}\,\mathrm{e}^{-x/\lambda}. \end{align}
Two obvious problems with this expression¶
Considering that morphogens determine cell fate in a concentration-dependent manner, there are two obvious problems with this expression.
- The total level of the morphogen is set by the strength of the source, $\eta$, assuming $D$ and $\gamma$ are constant. This would mean that cell fate would not be robust to the strength of the source.
- The decay length of the concentration profile is independent of size. This means that a larger organism would have all of its morphogen highly concentrated near $x/L = 0$.
Issue 2 is most clearly seen if we plot to profiles together.
# x-values for plot
L_1 = 100
L_2 = 200
x_1 = np.linspace(0, L_1, 200)
x_2 = np.linspace(0, L_2, 200)
# Value of lambda
lam = 25
# Make plots
p1 = bokeh.plotting.figure(height=250, width=350, x_axis_label='x/L', y_axis_label='conc')
p2 = bokeh.plotting.figure(height=250, width=350, x_axis_label='x/L', y_axis_type='log')
p1.line(x_1/L_1, np.exp(-x_1/lam), color='dodgerblue', legend='L=100', line_width=2)
p1.line(x_2/L_2, np.exp(-x_2/lam), color='tomato', legend='L=200', line_width=2)
p2.line(x_1/L_1, np.exp(-x_1/lam), color='dodgerblue', line_width=2)
p2.line(x_2/L_2, np.exp(-x_2/lam), color='tomato', line_width=2)
bokeh.io.show(bokeh.layouts.gridplot([[p1, p2]]))
We will address issue 1 in the homework. Now, we will consider a solution to issue 2 proposed by Ben-Zvi and Barkai, PNAS, 107, 6924-6929, 2010.
The concept of expansion repression¶
Consider another molecule, called an expander, that facilitates expansion of the morphogen. Specifically, the expander enhances diffusion of the morphogen and also prevents its degradation. However, and this is important because we need (you guessed it) feedback, the morphogen represses production of the expander. Our adjusted differential equations are now
\begin{align} \frac{\partial c_m}{\partial t} &= \frac{\partial}{\partial x}\,D_m(c_m, c_e)\,\frac{\partial c_m}{\partial x} - \gamma_m(c_m, c_e)\, c_m \\[1em] \frac{\partial c_e}{\partial t} &= D_e \,\frac{\partial^2 c_e}{\partial x^2} + \frac{\beta_e}{1 + (c_m/k)^n} - \gamma_e\, c_e. \\ &\phantom{blah} \nonumber \end{align}
We will consider the case where $D_m$ is constant, where the morphogen enhances its own degradation, and the expander represses degradation, shown schematically below.
The morphogen and expander are thus in a negative feedback loop. How might this help us with the issue that the morphogen gradient should scale?
Numerical solution to R-D equations¶
To solve R-D systems, I wrote the two functions below, which are the main engines for computation. The first computes the time derivative of concentrations at each of the spatial grid points. The second is a wrapper around scipy.integrate.odeint() to do the time stepping. Their doc strings describe them in more detail.
def dc_dt(c, t, x, derivs_0, derivs_L, diff_coeff_fun, diff_coeff_params,
rxn_fun, rxn_params, n_species, h):
"""
Time derivative of concentrations in an R-D system
for constant flux BCs.
Parameters
----------
c : ndarray, shape (n_species * n_gridpoints)
The concentration of the chemical species interleaved in a
a NumPy array. The interleaving allows us to take advantage
of the banded structure of the Jacobian when using the
Hindmarsh algorithm for integrating in time.
t : float
Time.
derivs_0 : ndarray, shape (n_species)
derivs_0[i] is the value of the diffusive flux,
D dc_i/dx, at x = 0, the leftmost boundary of the domain of x.
derivs_L : ndarray, shape (n_species)
derivs_0[i] is the value of the diffusive flux,
D dc_i/dx, at x = L, the rightmost boundary of the domain of x.
diff_coeff_fun : function
Function of the form diff_coeff_fun(c_tuple, t, x, *diff_coeff_params).
Returns an tuple where entry i is a NumPy array containing
the diffusion coefficient of species i at the grid points.
c_tuple[i] is a NumPy array containing the concentrations of
species i at the grid poitns.
diff_coeff_params : arbitrary
Tuple of parameters to be passed into diff_coeff_fun.
rxn_fun : function
Function of the form rxn_fun(c_tuple, t, *rxn_params).
Returns an tuple where entry i is a NumPy array containing
the net rate of production of species i by chemical reaction
at the grid points. c_tuple[i] is a NumPy array containing
the concentrations of species i at the grid poitns.
rxn_params : arbitrary
Tuple of parameters to be passed into rxn_fun.
n_species : int
Number of chemical species.
h : float
Grid spacing (assumed to be constant)
Returns
-------
dc_dt : ndarray, shape (n_species * n_gridpoints)
The time derivatives of the concentrations of the chemical
species at the grid points interleaved in a NumPy array.
"""
# Tuple of concentrations
c_tuple = tuple([c[i::n_species] for i in range(n_species)])
# Compute diffusion coefficients
D_tuple = diff_coeff_fun(c_tuple, t, x, *diff_coeff_params)
# Compute reaction terms
rxn_tuple = rxn_fun(c_tuple, t, *rxn_params)
# Return array
conc_deriv = np.empty_like(c)
# Convenient array for storing concentrations
da_dt = np.empty(len(c_tuple[0]))
# Useful to have square of grid spacing around
h2 = h**2
# Compute diffusion terms (central differencing w/ Neumann BCs)
for i in range(n_species):
# View of concentrations and diffusion coeff. for convenience
a = np.copy(c_tuple[i])
D = np.copy(D_tuple[i])
# Time derivative at left boundary
da_dt[0] = D[0] / h2 * 2 * (a[1] - a[0] - h * derivs_0[i])
# First derivatives of D and a
dD_dx = (D[2:] - D[:-2]) / (2*h)
da_dx = (a[2:] - a[:-2]) / (2*h)
# Time derivative for middle grid points
da_dt[1:-1] = D[1:-1] * np.diff(a, 2) / h2 + dD_dx*da_dx
# Time derivative at left boundary
da_dt[-1] = D[-1] / h2 * 2 * (a[-2] - a[-1] + h * derivs_L[i])
# Store in output array with reaction terms
conc_deriv[i::n_species] = da_dt + rxn_tuple[i]
return conc_deriv
def rd_solve(c_0_tuple, t, L=1, derivs_0=0, derivs_L=0, diff_coeff_fun=None,
diff_coeff_params=(), rxn_fun=None, rxn_params=(),
rtol=1.49012e-8, atol=1.49012e-8):
"""
Parameters
----------
c_0_tuple : tuple
c_0_tuple[i] is a NumPy array of length n_gridpoints with the
initial concentrations of chemical species i at the grid points.
t : ndarray
An array of time points for which the solution is desired.
L : float
Total length of the x-domain.
derivs_0 : ndarray, shape (n_species)
derivs_0[i] is the value of dc_i/dx at x = 0.
derivs_L : ndarray, shape (n_species)
derivs_L[i] is the value of dc_i/dx at x = L, the rightmost
boundary of the domain of x.
diff_coeff_fun : function
Function of the form diff_coeff_fun(c_tuple, t, *diff_coeff_params).
Returns an tuple where entry i is a NumPy array containing
the diffusion coefficient of species i at the grid points.
c_tuple[i] is a NumPy array containing the concentrations of
species i at the grid poitns.
diff_coeff_params : arbitrary
Tuple of parameters to be passed into diff_coeff_fun.
rxn_fun : function
Function of the form rxn_fun(c_tuple, t, *rxn_params).
Returns an tuple where entry i is a NumPy array containing
the net rate of production of species i by chemical reaction
at the grid points. c_tuple[i] is a NumPy array containing
the concentrations of species i at the grid poitns.
rxn_params : arbitrary
Tuple of parameters to be passed into rxn_fun.
rtol : float
Relative tolerance for solver. Default os odeint's default.
atol : float
Absolute tolerance for solver. Default os odeint's default.
Returns
-------
c_tuple : tuple
c_tuple[i] is a NumPy array of shape (len(t), n_gridpoints)
with the initial concentrations of chemical species i at
the grid points over time.
Notes
-----
.. When intergrating for long times near a steady state, you
may need to lower the absolute tolerance (atol) because the
solution does not change much over time and it may be difficult
for the solver to maintain tight tolerances.
"""
# Number of grid points
n_gridpoints = len(c_0_tuple[0])
# Number of chemical species
n_species = len(c_0_tuple)
# Grid spacing
h = (L / (n_gridpoints - 1))
# Grid points
x = np.linspace(0, L, n_gridpoints)
# Set up boundary conditions
if np.isscalar(derivs_0):
derivs_0 = np.array(n_species * [derivs_0])
if np.isscalar(derivs_L):
derivs_L = np.array(n_species * [derivs_L])
# Set up parameters to be passed in to dc_dt
params = (x, derivs_0, derivs_L, diff_coeff_fun, diff_coeff_params,
rxn_fun, rxn_params, n_species, h)
# Set up initial condition
c0 = np.empty(n_species * n_gridpoints)
for i in range(n_species):
c0[i::n_species] = c_0_tuple[i]
# Solve using odeint, taking advantage of banded structure
c = scipy.integrate.odeint(
dc_dt, c0, t, args=params, ml=n_species, mu=n_species,
rtol=rtol, atol=atol)
return tuple([c[:,i::n_species] for i in range(n_species)])
A sample calculation: the activator-substrate depletion model (ASDM)¶
We will use the R-D solver to solve the ASDM model, which exhibits Turing patterns. For the ASDM, we call $a$ the "activator" and the second species $s$ the "substrate." The reaction terms are
\begin{align} r_a &= \rho_a a^2 s - \mu_a a \\[2mm] r_s &= -\rho_s a^2 s + \sigma_s. \end{align}
Upon nondimensionalizing, we get
\begin{align} \frac{\partial a}{\partial t} &= d\, \frac{\partial^2a}{\partial x^2} + a^2s - a \\[2mm] \frac{\partial h}{\partial t} &= \frac{\partial^2h}{\partial x^2} + \mu(1 - a^2 s), \end{align}
giving the two parameters $\mu$ and $d = D_a/D_s$.
The ASDM has a unique homogeneous steady state of $a_0 = s_0 = 1$.
We will now write our functions to compute the diffusion coefficients and reaction terms for the ASDM, using the required I/O formats for the functions.
def constant_diff_coeffs(c_tuple, t, x, diff_coeffs):
n = len(c_tuple[0])
return tuple([diff_coeffs[i] * np.ones(n) for i in range(len(c_tuple))])
def asdm_rxn(as_tuple, t, mu):
"""
Reaction expression for activator-substrate depletion model.
Returns the rate of production of activator and substrate, respectively.
r_a = a**2 * s - a
r_s = mu * (1 - a**2 * s)
"""
# Unpack concentrations
a, s = as_tuple
# Compute and return reaction rates
a2s = a**2 * s
return (a2s - a, mu * (1.0 - a2s))
# Set up steady state (using 500 grid points)
a_0 = np.ones(500)
s_0 = np.ones(500)
# Seed the random number generator for reproducibility
np.random.seed(42)
# Make a small perturbation to a_0
a_0 += 0.001 * np.random.rand(len(a_0))
# Time points
t = np.linspace(0.0, 100.0, 100)
# Physical length of system
L = 20.0
# Diffusion coefficient
diff_coeffs = (0.05, 1.0)
# Reaction parameter (must be a tuple of params, even though only 1 for ASDM)
rxn_params = (1.4,)
# Solve
a, s = rd_solve((a_0, s_0), t, L=L, derivs_0=0, derivs_L=0,
diff_coeff_fun=constant_diff_coeffs,
diff_coeff_params=(diff_coeffs,), rxn_fun=asdm_rxn,
rxn_params=rxn_params)
Now that we've done the solution, let's plot the result with a slider so that we can change the values of the time points.
# Didn't need x to compute, need it to plot (remember: it's dimensionless)
x = np.linspace(0.0, L, len(a_0))
y_max = max(a.max(), s.max())
y_range = [-0.02*y_max, 1.02*y_max]
p = bokeh.plotting.figure(width=600, height=300, x_axis_label='x', y_axis_label='conc',
y_range=y_range)
act = p.line(x, a[0,:], color='dodgerblue', line_width=2, legend='activator')
sub = p.line(x, s[0,:], color='tomato', line_width=2, legend='substrate')
def update(time=0.0):
i = np.searchsorted(t, time)
act.data_source.data['y'] = a[i,:]
sub.data_source.data['y'] = s[i,:]
bokeh.io.push_notebook()
bokeh.io.show(p, notebook_handle=True)
ipywidgets.interactive(update,
time=ipywidgets.FloatSlider(min=t[0], max=t[-1], value=t[0]))
Numerical solution for simple source + degradation¶
We can solve for the temporal dynamics of the simple source + degradation system.
def degradation_rxn(c_tuple, t, gamma):
"""
Reaction expression for simple degradation.
"""
return (-gamma * c_tuple[0],)
# Time points
t = np.linspace(0.0, 10.0, 100)
# Physical length of system
L = 5.0
# Diffusion coefficient
diff_coeffs = (1.0,)
# Gamma
rxn_params = (1.0,)
# Initial concentration (500 grid points)
c_0_tuple = (np.zeros(500), )
# Production rate
eta = 1.0
derivs_0 = -eta / np.sqrt(diff_coeffs[0] * rxn_params[0])
# Solve
c = rd_solve(c_0_tuple, t, L=L, derivs_0=derivs_0, derivs_L=0,
diff_coeff_fun=constant_diff_coeffs,
diff_coeff_params=(diff_coeffs,), rxn_fun=degradation_rxn,
rxn_params=rxn_params)[0]
And let's make a plot.
# Make a movie and plot
# Didn't need x to compute, need it to plot (remember: it's dimensionless)
x = np.linspace(0.0, L, len(c_0_tuple[0]))
y_range = [-0.02*c.max(), 1.02*c.max()]
p = bokeh.plotting.figure(width=600, height=300, x_axis_label='x', y_axis_label='conc',
y_range=y_range)
c_line = p.line(x, c[0,:], color='dodgerblue', line_width=2)
def update(time=0.0):
i = np.searchsorted(t, time)
c_line.data_source.data['y'] = c[i,:]
bokeh.io.push_notebook()
bokeh.io.show(p, notebook_handle=True)
ipywidgets.interactive(update,
time=ipywidgets.FloatSlider(min=t[0], max=t[-1], value=t[0]))
Numerical solution of degradation with expander¶
To study the dynamics of the expander system, we will consider constant diffusion coefficients, and we take
\begin{align} \gamma_m(c_m, c_e) = (\gamma_{m,1} + \gamma_{m,2}c_m)\,\frac{c_m}{1+c_e}. \end{align}
This means that in addition to natural degradation of the morphogen, we also have enhanced autodegradation, both of which are retarded by the expander. Thus, our PDEs are
\begin{align} \frac{\partial c_m}{\partial t} &= D_m \,\frac{\partial^2 c_m}{\partial x^2} - (\gamma_{m,1} + \gamma_{m,2}c_m)\,\frac{c_m}{1+c_e}, \\[1em] \frac{\partial c_e}{\partial t} &= D_e \,\frac{\partial^2 c_e}{\partial x^2} + \frac{\beta_e}{1 + (c_m/k)^n} - \gamma_e\, c_e. \end{align}
Note the squared term in the degradation of the morphogen. This serves to make the gradient more robust to the source strength, as we will see in the homework.
We will perform the calculation for two values of $L$ and see how the gradient scales.
def expander_gradient_rxn(c_tuple, t, gamma_m1, gamma_m2, gamma_e,
beta_e, k, n):
"""
Reaction expression for simple degradation.
"""
# Unpack
c_m, c_e = c_tuple
# Morphogen reaction rate
m_rate = -(gamma_m1 + gamma_m2 * c_m) * c_m / (1 + c_e)
# Expander production rate
e_rate = beta_e / (1 + (c_m / k)**n) - gamma_e * c_e
return (m_rate, e_rate)
# Time points
t = np.linspace(0.0, 5e5, 300)
# Physical length of system and grid points
L_1 = 100.0
L_2 = 200.0
n_gridpoints = 100
# Diffusion coefficients
diff_coeffs = (1, 10)
# Reaction params
gamma_m1 = 0
gamma_m2 = 100
gamma_e = 1e-6
beta_e = 1e-2
k = 1e-3
n = 4
rxn_params = (gamma_m1, gamma_m2, gamma_e, beta_e, k, n)
# Initial concentration
c_0_tuple = (np.zeros(n_gridpoints), np.zeros(n_gridpoints))
# Production rate
eta = 10.0
# Solve
c_tuple_L_1 = rd_solve(
c_0_tuple, t, L=L_1, derivs_0=(-eta, 0), derivs_L=0,
diff_coeff_fun=constant_diff_coeffs, diff_coeff_params=(diff_coeffs,),
rxn_fun=expander_gradient_rxn, rxn_params=rxn_params, atol=1e-7)
c_tuple_L_2 = rd_solve(
c_0_tuple, t, L=L_2, derivs_0=(-eta, 0), derivs_L=0,
diff_coeff_fun=constant_diff_coeffs, diff_coeff_params=(diff_coeffs,),
rxn_fun=expander_gradient_rxn, rxn_params=rxn_params)
With the results in hand, we can compute the morphogen levels and see if we have scaling.
# x-values for plot
x_1 = np.linspace(0, L_1, n_gridpoints)
x_2 = np.linspace(0, L_2, n_gridpoints)
# Make plots
p1 = bokeh.plotting.figure(height=250, width=350, x_axis_label='x/L', y_axis_type='log')
p2 = bokeh.plotting.figure(height=250, width=350, x_axis_label='x', y_axis_type='log')
p1.line(x_1/L_1, c_tuple_L_1[0][-1,:], color='dodgerblue', legend='L=100', line_width=2)
p1.line(x_2/L_2, c_tuple_L_2[0][-1,:], color='tomato', legend='L=200', line_width=2)
p2.line(x_1, c_tuple_L_1[0][-1,:], color='dodgerblue', line_width=2)
p2.line(x_2, c_tuple_L_2[0][-1,:], color='tomato', line_width=2)
bokeh.io.show(bokeh.layouts.gridplot([[p1, p2]]));
Let's make a movie.
# Didn't need x to compute, need it to plot (remember: it's dimensionless)
x = np.linspace(0.0, L_1, n_gridpoints)
m = c_tuple_L_1[0]
e = c_tuple_L_1[1]
m /= m.max()
e /= e.max()
y_max = max(m.max(), e.max())
y_range = [-0.02, 1.02]
p = bokeh.plotting.figure(width=600, height=300, x_axis_label='x', y_axis_label='conc',
y_range=y_range)
morph = p.line(x, m[0,:], color='dodgerblue', line_width=2, legend='morphogen')
expan = p.line(x, e[0,:], color='tomato', line_width=2, legend='expander')
def update(time=0.0):
i = np.searchsorted(t, time)
morph.data_source.data['y'] = m[i,:]
expan.data_source.data['y'] = e[i,:]
bokeh.io.push_notebook()
bokeh.io.show(p, notebook_handle=True)
ipywidgets.interactive(update,
time=ipywidgets.FloatSlider(min=t[0], max=t[-1], value=t[0]))
Numerical solution for expander with non-cooperative degradation¶
We now consider the case where we do not have cooperative degradation. I.e., we have
\begin{align} \gamma_m(c_m, c_e) = \gamma_{m,1}\,\frac{c_m}{1+c_e}. \end{align}
Thus, our PDEs are
\begin{align} \frac{\partial c_m}{\partial t} &= D_m \,\frac{\partial^2 c_m}{\partial x^2} - \gamma_{m,1}\frac{c_m}{1+c_e} \\[1em] \frac{\partial c_e}{\partial t} &= D_e \,\frac{\partial^2 c_e}{\partial x^2} + \frac{\beta_e}{1 + (c_m/k)^n} - \gamma_e\, c_e. \end{align}
We can reuse the functions we already wrote with $\gamma_{m,1} >0$ and $\gamma_{m,2} = 0$.
# Reaction params
gamma_m1 = 0.1
gamma_m2 = 0
gamma_e = 1e-6
beta_e = 1e-2
k = 1e-3
n = 4
rxn_params = (gamma_m1, gamma_m2, gamma_e, beta_e, k, n)
# Initial concentration
c_0_tuple = (np.zeros(n_gridpoints), np.zeros(n_gridpoints))
# Production rate
eta = 0.1
# Solve
c_tuple_L_1 = rd_solve(
c_0_tuple, t, L=L_1, derivs_0=(-eta, 0), derivs_L=0,
diff_coeff_fun=constant_diff_coeffs, diff_coeff_params=(diff_coeffs,),
rxn_fun=expander_gradient_rxn, rxn_params=rxn_params, atol=1e-6)
c_tuple_L_2 = rd_solve(
c_0_tuple, t, L=L_2, derivs_0=(-eta, 0), derivs_L=0,
diff_coeff_fun=constant_diff_coeffs, diff_coeff_params=(diff_coeffs,),
rxn_fun=expander_gradient_rxn, rxn_params=rxn_params, atol=1e-6)
# x-values for plot
x_1 = np.linspace(0, L_1, n_gridpoints)
x_2 = np.linspace(0, L_2, n_gridpoints)
# Make plots
p1 = bokeh.plotting.figure(height=250, width=350, x_axis_label='x/L', y_axis_type='log')
p2 = bokeh.plotting.figure(height=250, width=350, x_axis_label='x', y_axis_type='log')
p1.line(x_1/L_1, c_tuple_L_1[0][-1,:], color='dodgerblue', legend='L=100', line_width=2)
p1.line(x_2/L_2, c_tuple_L_2[0][-1,:], color='tomato', legend='L=200', line_width=2)
p2.line(x_1, c_tuple_L_1[0][-1,:], color='dodgerblue', line_width=2)
p2.line(x_2, c_tuple_L_2[0][-1,:], color='tomato', line_width=2)
bokeh.io.show(bokeh.layouts.gridplot([[p1, p2]]));
Again, we get scaling of the gradient slope, but we see that the absolute level varies.