In [1]:
import numpy as np
import matplotlib.pyplot as plt
def find_intersection(lineA, lineB, nA, nB):
dAX = lineA[nA + 1, 0] - lineA[nA, 0]
dAY = lineA[nA + 1, 1] - lineA[nA, 1]
dBX = lineB[nB + 1, 0] - lineB[nB, 0]
dBY = lineB[nB + 1, 1] - lineB[nB, 1]
dABX = lineB[nB, 0] - lineA[nA, 0]
dABY = lineB[nB, 1] - lineA[nA, 1]
d = dAX * dBY - dAY * dBX
alpha = (dABX * dBY - dABY * dBX) / d
beta = (dABX * dAY - dABY * dAX) / d
# check whether there is an intersection in between
if (0 <= alpha <= 1) and (0 <= beta <= 1):
intersection = np.array([(1 - alpha) * lineA[nA, 0] + alpha * lineA[nA + 1, 0],
(1 - alpha) * lineA[nA, 1] + alpha * lineA[nA + 1, 1],
1])
# 1st and 2nd elements: X & Y
# 3rd element: 1 means found, 0 means not found
else:
intersection = np.array([alpha, beta, 0])
# 1st and 2nd elements: alpha & beta
return intersection
# test example with three points in each line. They intersect by the first segments.
lineA = np.array([[0, 0], [1, 0], [2, 0]])
lineB = np.array([[0.5, 1], [0.5, -1], [0.5, -3]])
plt.figure()
plt.plot(lineA[:, 0], lineA[:, 1], color='red', marker='o', label='Line A', linewidth=2)
plt.plot(lineB[:, 0], lineB[:, 1], color='blue', marker='o',label='Line B', linewidth=2)
plt.xlabel("X")
plt.ylabel("Y")
plt.legend(loc='upper right')
plt.show()
print(find_intersection(lineA, lineB, 0, 0)) # Intersect at (0.5, 0)
print(find_intersection(lineA, lineB, 0, 1)) # No intersection
print(find_intersection(lineA, lineB, 1, 1)) # No intersection
[0.5 0. 1. ] [ 0.5 -0.5 0. ] [-0.5 -0.5 0. ]
Find all intersections of two nullclines in Part 3A¶
In [2]:
def find_intersection_all(lineA, lineB):
lineA_ind = np.arange(lineA.shape[0] - 1)
lineB_ind = np.arange(lineB.shape[0] - 1)
lines_all = np.array(np.meshgrid(lineA_ind, lineB_ind)).T.reshape(-1, 2)
# all combination of indices for lineA and lineB
def find_intersection_single(inds):
return find_intersection(lineA, lineB, inds[0], inds[1])
results = np.apply_along_axis(find_intersection_single, 1, lines_all)
intersections = results[results[:, 2] == 1, :2]
return intersections
# Load nullcline data from Part 3A
null1 = np.loadtxt("null1_03A.csv", delimiter=',')
null2 = np.loadtxt("null2_03A.csv", delimiter=',')
# Find intersections
ss_all = find_intersection_all(null1, null2)
print(ss_all)
[[542.37185133 35.27209876] [189.96276765 126.64623313] [ 52.90805706 361.58503595]]
Find all intersections (improved performance)¶
In [3]:
# Same derivative function from Part 3A
def hill_inh(X, X_th, n):
# inhibitory Hill function
# X_th: Hill threshold, n: Hill coefficient
a = (X / X_th)**n
return 1 / (1 + a)
def derivs_ts(t, Xs):
# Calculate derivative function for a toggle switch circuit
X = Xs[0]
Y = Xs[1]
dxdt = 5 + 50 * hill_inh(Y, 100, 4) - 0.1 * X
dydt = 4 + 40 * hill_inh(X, 150, 4) - 0.12 * Y
return [dxdt, dydt]
# Calculate dX/dt along the nullcline dY/dt = 0
dxdt_y_null = np.apply_along_axis(lambda Xs: derivs_ts(0, Xs)[0], 1, null2)
# Calculate dY/dt along the nullcline dX/dt = 0
dydt_x_null = np.apply_along_axis(lambda Xs: derivs_ts(0, Xs)[1], 1, null1)
def find_intersection_all_fast(lineA, lineB, dxdt, dydt):
small = 1e-3 # used to allow some tolerance of numerical errors in detecting sign flipping
lineA_ind = np.where(dydt * np.concatenate([dydt[1:], [np.nan]]) <= small)[0]
lineB_ind = np.where(dxdt * np.concatenate([dxdt[1:], [np.nan]]) <= small)[0]
lines_all = np.array(np.meshgrid(lineA_ind, lineB_ind)).T.reshape(-1, 2)
def find_intersection_single(inds):
return find_intersection(lineA, lineB, inds[0], inds[1])
results = np.apply_along_axis(find_intersection_single, 1, lines_all)
intersections = results[results[:, 2] == 1, :2]
return intersections
# Find intersections
ss_all = find_intersection_all_fast(null1, null2, dxdt_y_null, dydt_x_null)
print(ss_all)
[[542.37185133 35.27209876] [189.96276765 126.64623313] [ 52.90805706 361.58503595]]
Stability¶
In [4]:
from scipy.linalg import eig
import pandas as pd
def stability_2D(derivs, ss, **kwargs):
delta = 0.001
f_current = derivs(0, ss, **kwargs)
f_plus_dx = derivs(0, ss + np.array([delta, 0]), **kwargs)
f_plus_dy = derivs(0, ss + np.array([0, delta]), **kwargs)
# Finite difference to approximate Jacobian
dfxdx = (f_plus_dx[0] - f_current[0]) / delta
dfxdy = (f_plus_dy[0] - f_current[0]) / delta
dfydx = (f_plus_dx[1] - f_current[1]) / delta
dfydy = (f_plus_dy[1] - f_current[1]) / delta
jacobian = np.array([[dfxdx, dfxdy], [dfydx, dfydy]])
lambdas = eig(jacobian)[0]
if np.any(np.iscomplex(lambdas)):
if np.real(lambdas[0]) < 0:
stability = 4 # stable spiral
else:
stability = 5 # unstable spiral
else:
if (lambdas[0] < 0) and (lambdas[1] < 0):
stability = 1 # stable
elif (lambdas[0] >= 0) and (lambdas[1] >= 0):
stability = 2 # unstable
else:
stability = 3 # saddle
return np.concatenate([ss, [stability]])
ss_with_stability = np.apply_along_axis(lambda ss: stability_2D(derivs_ts, ss), 1, ss_all)
results = pd.DataFrame(ss_with_stability, columns=["X", "Y", "Stability"])
print(results)
X Y Stability 0 542.371851 35.272099 1.0 1 189.962768 126.646233 3.0 2 52.908057 361.585036 1.0