In [ ]:
%pylab inline
Populating the interactive namespace from numpy and matplotlib
Let's start with the double negative feedback loop, these are the equations for this system:
$$\begin{aligned} \frac{dx}{dt} = \frac{k_1*K^n}{y^n+K^n} - k_2*x\\ \frac{dy}{dt} = \frac{k_3*K^n}{x^n+K^n} - k_4*y \end{aligned}$$
Let's start by implementing euler's method to plot a couple of trayectories.
We are going to use $k_1= k_2= k_3= k_4 = 1$ and K = 0.5
Also we will be looking at the case of $n=1$
In [ ]:
dt = 1e-4
tmax = 10
t = np.arange(0,tmax+dt,dt)
numt = len(t)
n = 1
def dxdt(x,y,n,k1=1,k2=1,K=0.5):
return (k1*K**n)/(y**n+K**n) - (k2*x)
def dydt(x,y,n,k3=1,k4=1,K=0.5):
return (k3*K**n)/(x**n+K**n) - (k4*y)
In [ ]:
dt
Out[Â ]:
0.0001
In [ ]:
def dm(x_i, y_i,n):
x = np.zeros(numt)
y = np.zeros(numt)
x += x_i
y += y_i
for i in range(1,numt):
x[i] = x[i-1] + (dt * dxdt(x[i-1],y[i-1],n))
y[i] = y[i-1] + (dt * dydt(x[i-1],y[i-1],n))
return x,y
In [ ]:
init_conds = np.linspace(0,1,25)
fig,ax = subplots(ncols= 1,nrows =1)
fig.set_figwidth(10)
fig.set_figheight(10)
n=1
for i in range(25):
ini = np.random.choice(init_conds, 2)
ax.plot(dm(ini[0], ini[1],n)[0],dm(ini[0], ini[1],n)[1])
ax.scatter(ini[0], ini[1], facecolors = 'none',edgecolors = 'k', s=30)
ax.set_ylim(-0.1,1.1)
ax.set_title('25 different trayectories, double negative feedback, n:{}'.format(n));
ax.set_ylabel('concentration of y')
ax.set_xlabel('concentration of x')
Out[Â ]:
Text(0.5, 0, 'concentration of x')
The first thing that strikes to mind is that for the case of $n=1$ we have one stable steady state. all trayectories, depicted by the different color lines, start at the black circumferences and seem to end in the center. We could say that this has a stable steady state. However, to confirm these findings we can draw the nullclines and the vector field, which we'll do next:
$$0 = \frac{k_1 K^n}{y^n+K^n} - (k_2x) $$ $$(k_2x) = \frac{k_1 K^n}{y^n+K^n}$$ $$x = \frac{k_1 K^n}{k_2(y^n+K^n)}$$
Similarly, for $\frac{dy}{dt}$:
$$y = \frac{k_3 K^n}{k_4(x^n+K^n)}$$
In [ ]:
np.meshgrid(np.linspace(0,1,9),np.linspace(0,1,9))
Out[Â ]:
[array([[0. , 0.125, 0.25 , 0.375, 0.5 , 0.625, 0.75 , 0.875, 1. ],
[0. , 0.125, 0.25 , 0.375, 0.5 , 0.625, 0.75 , 0.875, 1. ],
[0. , 0.125, 0.25 , 0.375, 0.5 , 0.625, 0.75 , 0.875, 1. ],
[0. , 0.125, 0.25 , 0.375, 0.5 , 0.625, 0.75 , 0.875, 1. ],
[0. , 0.125, 0.25 , 0.375, 0.5 , 0.625, 0.75 , 0.875, 1. ],
[0. , 0.125, 0.25 , 0.375, 0.5 , 0.625, 0.75 , 0.875, 1. ],
[0. , 0.125, 0.25 , 0.375, 0.5 , 0.625, 0.75 , 0.875, 1. ],
[0. , 0.125, 0.25 , 0.375, 0.5 , 0.625, 0.75 , 0.875, 1. ],
[0. , 0.125, 0.25 , 0.375, 0.5 , 0.625, 0.75 , 0.875, 1. ]]),
array([[0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ],
[0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125],
[0.25 , 0.25 , 0.25 , 0.25 , 0.25 , 0.25 , 0.25 , 0.25 , 0.25 ],
[0.375, 0.375, 0.375, 0.375, 0.375, 0.375, 0.375, 0.375, 0.375],
[0.5 , 0.5 , 0.5 , 0.5 , 0.5 , 0.5 , 0.5 , 0.5 , 0.5 ],
[0.625, 0.625, 0.625, 0.625, 0.625, 0.625, 0.625, 0.625, 0.625],
[0.75 , 0.75 , 0.75 , 0.75 , 0.75 , 0.75 , 0.75 , 0.75 , 0.75 ],
[0.875, 0.875, 0.875, 0.875, 0.875, 0.875, 0.875, 0.875, 0.875],
[1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. ]])]
In [ ]:
dxn1
Out[Â ]:
array([[ 1. , 0.875 , 0.75 , 0.625 , 0.5 ,
0.375 , 0.25 , 0.125 , 0. ],
[ 0.8 , 0.675 , 0.55 , 0.425 , 0.3 ,
0.175 , 0.05 , -0.075 , -0.2 ],
[ 0.66666667, 0.54166667, 0.41666667, 0.29166667, 0.16666667,
0.04166667, -0.08333333, -0.20833333, -0.33333333],
[ 0.57142857, 0.44642857, 0.32142857, 0.19642857, 0.07142857,
-0.05357143, -0.17857143, -0.30357143, -0.42857143],
[ 0.5 , 0.375 , 0.25 , 0.125 , 0. ,
-0.125 , -0.25 , -0.375 , -0.5 ],
[ 0.44444444, 0.31944444, 0.19444444, 0.06944444, -0.05555556,
-0.18055556, -0.30555556, -0.43055556, -0.55555556],
[ 0.4 , 0.275 , 0.15 , 0.025 , -0.1 ,
-0.225 , -0.35 , -0.475 , -0.6 ],
[ 0.36363636, 0.23863636, 0.11363636, -0.01136364, -0.13636364,
-0.26136364, -0.38636364, -0.51136364, -0.63636364],
[ 0.33333333, 0.20833333, 0.08333333, -0.04166667, -0.16666667,
-0.29166667, -0.41666667, -0.54166667, -0.66666667]])
In [ ]:
nulldx(0,1)
Out[Â ]:
1.0
In [ ]:
y_ = np.linspace(0.5,1, 5)
x_ = nulldy(y_,1)
In [ ]:
def nulldx(y, n):
return (k1*K**n)/(k2*(y**n+K**n))
def nulldy(x, n):
return (k3*K**n)/(k4*(x**n+K**n))
x = np.linspace(0,1,200)
y = np.linspace(0,1,200)
X,Y =np.meshgrid(np.linspace(0,1,9),np.linspace(0,1,9))
k1,k2,k3,k4 =1,1,1,1
K= 0.5
n = 1
dxn1 = (k1*K**n)/(Y**n+K**n) - (k2*X)
dyn1 = (k3*K**n)/(X**n+K**n) - (k4*Y)
# nulldx(x, n)
# nulldy(x, n)
fig,ax = subplots(ncols= 1,nrows =1)
fig.set_figwidth(10)
fig.set_figheight(10)
for pair in zip(x_,y_):
ini = pair#np.random.choice(init_conds,2)
ax.plot(dm(ini[0], ini[1],n)[0],dm(ini[0], ini[1],n)[1])
ax.scatter(ini[0], ini[1], facecolors = 'none',edgecolors = 'k', s=30)
ax.set_ylim(-0.1,1.1)
ax.set_title('25 different trayectories, double negative feedback, n:{}'.format(n));
ax.plot(x,nulldy(x, n=1), c = 'y',ls= '--', alpha =0.5)
ax.plot(nulldx(x, n=1),y, c = 'g',ls= '--', alpha =0.5)
ax.set_title('nullclines in phase space')
ax.quiver(X,Y,dxn1,dyn1, alpha = 0.3, pivot='mid')
ax.set_ylabel('concentration of y')
ax.set_xlabel('concentration of x')
Out[Â ]:
Text(0.5, 0, 'concentration of x')
The PhoP/PhoQ system¶
The PhoP/PhoQ operon is in the bacterium Salmonella enterica and is a good example of a positive feedback loop.
- PhoP comes in two phosphorylation states: in our model, $P$ is phosphorylated PhoP, and $U$ is unphorphorylated
- There is a constitutive promoter (which is ON in all circumstances), $P_1$
- There is also an inducible promoter, which is induced by the presence of $P$ -- phosphorylated PhoP
Differential equations: $$ \begin{align} &\frac{dU}{dt} = k_1 + k_2\frac{P^n}{K^n+P^n}+k_a^-P-(k_d+k_a^+)U\\ &\frac{dP}{dt}=k_a^+U-(k_d+k_a^-)P \end{align} $$
Let's go through the terms slowly:
- $U$:
- $k_1$ describes the rate of the constitutive promoter $P_1$
- $k_2$ is the rate of the inducible promoter, which has positive feedback from phosphorylated $P$, which we model with a Hill function
- $U$ phosphorylates to $P$ with rate $k_a^+$
- $P$ dephosphorylates to $U$ with rate $k_a^-$
- $U$ degrades with rate $k_d$
- $P$:
- $U$ phosphorylates to $P$ with rate $k_a^+$
- $P$ dephosphorylates to $U$ with rate $k_a^-$
- $P$ degrades with rate $k_d$
In [ ]:
## function to solve deterministically using Euler's method
def solve_ode(sol0, params, n=1, tspan=[0,1], dt=0.001):
# retrieve our kinetic rates
[k1, k2, kaplus, kaminus, kd, K] = params
# our differential equations
def dudt(u, p):
return k1 + k2*p**n/(K**n+p**n) + kaminus*p - (kd + kaplus)*u
def dpdt(u, p):
return kaplus*u - (kd + kaminus)*p
# set up time array
T = np.arange(tspan[0],tspan[-1]+dt,dt)
# set our initial conditions
u = [sol0[0]]
p = [sol0[1]]
# go through timepoint by timepoint, and compute
# u_n+1 = u_n + du/dt (u_n, p_n),
# p_n+1 = p_n + dp/dt (p_n, p_n)
# (euler's method)
for t in T[0:-1]:
u.append(u[-1] + dt * dudt(u[-1],p[-1]))
p.append(p[-1] + dt * dpdt(u[-1],p[-1]))
return T,u,p
def plot_PhoP(params, n, T, u, p):
# plot time evolution and phase plane
fig, axs = plt.subplots(1, 2, figsize=(14,5), gridspec_kw={'width_ratios': [2, 1]})
plt.suptitle('$k_1 = {}, ~k_2 = {}, ~k_a^+ = {}, ~k_a^- = {}, ~k_d = {}, ~K = {}, ~n = {}$'.format(*params, n))
axs[0].set_title('time evolution')
axs[0].plot(T,u,label="U")
axs[0].plot(T,p,label="P")
axs[0].set_xlabel('time (min)')
axs[0].set_ylabel('concentration ($\mu$M)')
axs[0].legend(loc='upper left', bbox_to_anchor=(0, -0.2), borderaxespad=0)
axs[1].set_title('phase plane')
axs[1].plot(u,p)
axs[1].scatter(u[0], p[0], facecolors='none', color='k', s=30,
label=f'initial (u,p): {u[0], p[0]}')
axs[1].axis('equal')
axs[1].set_xlabel('$U$ ($\mu$M)'); axs[1].set_ylabel('$P$ ($\mu$M)')
axs[1].legend(loc='upper left', bbox_to_anchor=(0, -0.2), borderaxespad=0)
plt.show()
Let's define some kinetic parameters. These are what are given in the lecture notes page:
In [ ]:
k1 = 0.01
k2 = 0.0
kaplus = 0.5
kaminus = 0.20
kd = 0.08
K = 0.2
n = 2
params = [k1,k2,kaplus,kaminus,kd,K]
Now define initial conditions and solve the ODE:
In [ ]:
init_up = (0.1, 0.1)
T,u,p = solve_ode(init_up, params, n, tspan=(0, 100), dt=0.001)
Plot $U$ and $P$ through time (left) and phase plane (right)
In [ ]:
plot_PhoP(params, n, T, u, p)
The example above has $k_2 = 0$, meaning there is no inducible promoter. Both the phosphorylated and unphosphorylated forms remain low, and there is some change in their concentrations via $U \leftarrow \rightarrow P$ conversion.
Now let's set $k_2 = 0.3$ and turn on the inducible promoter:
In [ ]:
k1 = 0.01
k2 = 0.3
kaplus = 0.5
kaminus = 0.20
kd = 0.08
K = 0.2
n = 2
params = [k1,k2,kaplus,kaminus,kd,K]
init_up = (0.1, 0.1)
T,u,p = solve_ode(init_up, params, n, tspan=(0, 100), dt=0.001)
plot_PhoP(params, n, T, u, p)
Great, now that the inducible promoter is on, $P$ activates the inducible promoter, and the concentrations of both forms of the protein increase. We have recovered the example on the lecture notes page.
Let's run multiple integrations using different initial conditions. The plot below is a combined phase plane for multiple trajectories, each starting at different initial concentrations of $U$ and $P$ (open circles). The final $U$ and $P$ is marked with a solid black circle:
In [ ]:
n_trajs = 10
plt.figure()
for i in range(n_trajs):
# make sure we have a small set and a large set of initial concentrations
if i == 0:
init_up = (0.1, 0.1)
elif i == 1:
init_up = (4, 4)
elif i == 2:
init_up = (2.2, 1)
else:
init_up = np.random.uniform(0, 10, 2)
T,u,p = solve_ode(init_up, params, n=2, tspan=(0, 100), dt=0.001)
plt.scatter(u, p, s=2)
plt.scatter(u[0], p[0], facecolors='none', color='k', s=30)
plt.scatter(u[-1], p[-1], color='k', s=30)
plt.xlabel('$U$ ($\mu$M)')
plt.ylabel('$P$ ($\mu$P)')
us = np.linspace(0, 15, 100)
ps = kaplus/(kd + kaminus) * us
#plt.plot(us, ps, c='k', label=r'nulcline: $P = \frac{k_a^+}{k_d + k_a^-}U$', ls='--')
#plt.legend(loc='upper left', bbox_to_anchor=(1.03, 1), borderaxespad=0)
#plt.xlim(2,3)
#plt.ylim(0.5, 1)
plt.show()
In [ ]:
With these model parameters, there is one stable fixed point -- all the trajectories settle to the same final concentration.
Now, let's try $n=3$:
In [ ]:
n_trajs = 10
plt.figure()
for i in range(n_trajs-3):
# make sure we have a small set and a large set of initial concentrations
if i == 0:
init_up = (0.1, 0.1)
elif i == 1:
init_up = (4, 4)
elif i == 2:
init_up = (0.5, 0.5)
else:
init_up = np.random.uniform(0, 10, 2)
T,u,p = solve_ode(init_up, params, n=3, tspan=(0, 100), dt=0.001)
plt.plot(u, p)
plt.scatter(u[0], p[0], facecolors='none', color='k', s=30)
plt.scatter(u[-1], p[-1], color='k', s=30)
plt.xlabel('$U$ ($\mu$M)')
plt.ylabel('$P$ ($\mu$P)')
us = np.linspace(0, 15, 100)
ps = kaplus/(kd + kaminus) * us
#plt.plot(us, ps, c='k', label=r'nulcline: $P = \frac{k_a^+}{k_d + k_a^-}U$', ls='--')
#plt.legend(loc='upper left', bbox_to_anchor=(1.03, 1), borderaxespad=0)
plt.show()
--------------------------------------------------------------------------- OverflowError Traceback (most recent call last) <ipython-input-61-60c18706097b> in <module>() 14 else: 15 init_up = np.random.uniform(0, 10, 2) ---> 16 T,u,p = solve_ode(init_up, params, n=3, tspan=(0, 100), dt=0.001) 17 18 plt.plot(u, p) <ipython-input-55-f34549ea1760> in solve_ode(sol0, params, n, tspan, dt) 23 # (euler's method) 24 for t in T[0:-1]: ---> 25 u.append(u[-1] + dt * dudt(u[-1],p[-1])) 26 p.append(p[-1] + dt * dpdt(u[-1],p[-1])) 27 <ipython-input-55-f34549ea1760> in dudt(u, p) 7 # our differential equations 8 def dudt(u, p): ----> 9 return k1 + k2*p**n/(K**n+p**n) + kaminus*p - (kd + kaplus)*u 10 def dpdt(u, p): 11 return kaplus*u**2 - (kd + kaminus)*p OverflowError: (34, 'Numerical result out of range')
Positive feedback can produce bistability¶
With $n=3$, there is suddenly a second stable steady state that appears at low concentrations!
Let's make a quiver plot:
In [ ]:
init_up = (0.1, 0.1)
T,u,p = solve_ode(init_up, params, n=3, tspan=(0, 100), dt=0.001)
small_steady = u[-1], p[-1]
init_up = (0.4, 0.4)
T,u,p = solve_ode(init_up, params, n=3, tspan=(0, 100), dt=0.001)
large_steady = u[-1], p[-1]
small_steady, large_steady
Out[Â ]:
((0.10730117351711857, 0.026718498138532512), (3.0465959645539717, 0.7586127212686212))
In [ ]:
def phase_diagram(params, n, xlims=[0,1], ylims=[0,1], dx=0.1, dy=0.1):
# retrieve our kinetic rates
[k1, k2, kaplus, kaminus, kd, K] = params
# set up x and y arrays
x = np.arange(xlims[0],xlims[-1]+dx,dx)
y = np.arange(ylims[0],ylims[-1]+dy,dy)
Y,X = np.mgrid[ylims[0]:ylims[-1]:dy,xlims[0]:xlims[-1]:dx]
# get directions in U, P
U = k1 + k2*Y**n/(K**n+Y**n) + kaminus*Y - (kd + kaplus)*X
P = kaplus*X-(kd + kaminus)*Y
# quiver or streamplot
speed = np.sqrt(U**2 + P**2)
lw = 5*speed / speed.max()
#plt.quiver(X,Y,U,P) # alternative method by Carlos
plt.streamplot(X,Y,U,P,density=1,linewidth=lw)
plt.xlabel("U")
plt.ylabel("P")
Phase plane:
In [ ]:
plt.figure()
phase_diagram(params,n=3,xlims=[0,10],ylims=[0,5],dx=0.2,dy=0.1)
plt.scatter(*small_steady, c='k')
plt.scatter(*large_steady, c='k')
plt.show()
Zooming in to the small steady state:
In [ ]:
plt.figure()
phase_diagram(params,n=3,xlims=[0,0.2],ylims=[0,0.1],dx=0.01,dy=0.005)
plt.scatter(*small_steady, c='k')
#plt.scatter(*large_steady, c='k')
plt.show()
