<jupyter_start><jupyter_text>Bistable Double Negative Switch  Integrating individual Trajectory  Defining the system<jupyter_code>#Symbolic system
x,y,fxy, fyx,  kyx, kxy, bx, by, dx, dy, n = sm.symbols(
    'x,y,fyx, fxy, kyx, kxy, bx, by, dx, dy, n'
)


#System of equations
xdot = bx + fyx*kyx**n/(kyx**n + y**n) - dx*x
ydot = by + fxy*kxy**n/(kxy**n + x**n) - dy*y

#lambdify expression
dxdt = sm.lambdify((x,y, fyx, fxy,  kyx, kxy, bx, by, dx, dy, n), xdot)
dydt = sm.lambdify((x,y,fyx, fxy,  kyx, kxy, bx, by, dx, dy, n), ydot)

#Define system in single function
def double_neg(z0, t, *p):
    x,y = z0
    del_x = dxdt(x,y,*p)
    del_y = dydt(x,y,*p)
    return [del_x, del_y]<jupyter_output><empty_output><jupyter_text>Integrating potential path along said system<jupyter_code>#Setting ODE Parameters
p = (2, 2, 0.7, 0.5, 0.2, 0.2, 1, 1, 4)#fyx, fxy, kyx, kxy, b_x, b_y, d_x, d_y, n

#Integration Parameters
initials = [50, 50] #x0, y0
t0 = 0 #initial time
v0 = 0 #inital potential value
tol = 1e-5 #convergence criterion
tstep = 0.05 #time step



xy_path, v_path, times = integrate_potential_path(double_neg, initials, t0, v0, tol, tstep, *p)<jupyter_output><empty_output><jupyter_text>Plottting Results<jupyter_code>#Plot the trajectories with respect to time
x_path, y_path = zip(*xy_path)
plt.plot(times, np.log1p(x_path), label = "X")
plt.plot(times, np.log1p(y_path), label = "Y")
# plt.plot(times, v_path, label = "V")
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Approximating the Potential Surface<jupyter_code>#Setting up grid parameters
xy_min = 0.0
xy_max = 15.0
xy_step = 0.5

#integration parameters
int_params = (t0, v0, tol, tstep)

#Calculating and aligning Quasipotential trajectories with different inits along a 2d grid
x_paths, y_paths, v_paths, attractors = QuasiPotential2D(double_neg, p, int_params, xy_min, 
                                                         xy_max, xy_step,attractor_tol = 1e-2,  n_jobs = -1)
print("num attractors = ", len(attractors))<jupyter_output>num attractors =  2<jupyter_text>Interpolate through trajectories to map surface<jupyter_code>xs = []
ys = []
Vs = []

for x, y, v in zip(x_paths, y_paths, v_paths):
    xs += list(x)
    ys += list(y)
    Vs += list(v)
    
gridlines = 100

x_grid = np.linspace(min(xs), max(xs), gridlines)
y_grid = np.linspace(min(ys), max(ys), gridlines)

X_grid,Y_grid = np.meshgrid(x_grid, y_grid)
V_grid = griddata((xs,ys),Vs, (X_grid, Y_grid))
fig = plt.figure(dpi = 150)
ax = fig.add_subplot(111, projection='3d')

ax.plot_surface(X_grid, Y_grid, V_grid, cmap="cividis", lw=0.5, rstride=1, cstride=1)
ax.contour(X_grid, Y_grid, V_grid, 30, linestyles="solid")

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('V')

ax.view_init(30, -140)

plt.show()<jupyter_output><empty_output><jupyter_text>Stochastic Simulations  Simulate and plot one trajectory<jupyter_code># Code adapted from Sueskind: https://github.com/sueskind/gillespie
random.seed(a = 42)

x0 = [10, 10]
fold_yx, fold_xy,  k_yx, k_xy, b_x, b_y, d_x, d_y, n = (
    2, 2, 0.7, 0.5, 0.2, 0.2, 1, 1, 4
)

# Propensities
#  - A list of functions that, given the values of your current states, compute the transition rate of a reaction:
def x_birth(x,y):
    return b_x + fold_yx*k_yx**n/(k_yx**n + y**n)
def y_birth(x,y):
    return b_y + fold_xy*k_xy**n/(k_xy**n + x**n)
def x_death(x,y):
    return d_x*x
def y_death(x,y):
    return d_y*y

propensities = [
    x_birth,
    y_birth,
    x_death,
    y_death
]

# Stoichiometries
# Positive terms in the differential equations are considered as "birth" (+1 in stoich)
# Negative terms in the differential equations are considered as "death" (-1 in stoich)
stoichiometry = [
    [1,0],
    [0,1],
    [-1,0],
    [0,-1]
         ]
#Initialize Number of runs
runs = 10000
t_runs = {}
X_runs = {}
Y_runs = {}
V_runs = {}

#Define parallelizable gillespie function for asynchronous simulations
def gillespie_sim(r, initials, propoensities, stoichiometry, duration):
    times, XY_stoch = gillespie.simulate(initials, propensities, stoichiometry, duration)
    return (r, times, XY_stoch)

#Run simulations in parallel
pool = mp.Pool(mp.cpu_count())
results = []
results = pool.starmap_async(gillespie_sim, [(r, initials, propensities, stoichiometry, 50) for r in range(runs)]).get()

#Close pool after finished
pool.close()<jupyter_output><empty_output><jupyter_text>Calculate potential for each point along trajectory<jupyter_code>x_grid_step = (max(xs) - min(xs)) / gridlines
print(x_grid_step)
y_grid_step = (max(ys) - min(ys)) / gridlines
print(y_grid_step)
#Collect results into dictionary
for result in results:
    r = result[0]
    t_runs[r] = result[1]
    X_stoch, Y_stoch = zip(*result[2])
    X_runs[r] = X_stoch
    Y_runs[r] = Y_stoch
    V_stoch = []
    for i in range(0, len(X_stoch)):
        x_idx = int(X_stoch[i] / x_grid_step)
        y_idx = int(Y_stoch[i] / y_grid_step)
        V_stoch.append(V_grid[y_idx, x_idx])
    V_runs[r] = V_stoch
max([max(Y_runs[r]) for r in range(runs)])<jupyter_output><empty_output><jupyter_text>Overlay stoch traj on potential surface<jupyter_code>fig = plt.figure(dpi = 150)
ax = fig.add_subplot(111, projection='3d')

ax.plot_surface(X_grid, Y_grid, V_grid, cmap="cividis", lw=0.5, rstride=1, cstride=1, alpha=0.3)
ax.contour(X_grid, Y_grid, V_grid, 30, linestyles="solid", alpha = 0.3)

ax.plot(X_stoch, Y_stoch, V_stoch, color="red")
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('V')

ax.view_init(30, -140)

plt.show()
#plt.savefig('Lscape_tristab_n10_By4-0_stochTraj.png', bbox_inches='tight')
coords = {"X":[X_runs[r][-1] for r in range(runs)], "Y":[Y_runs[r][-1] for r in range(runs)], "V":[V_runs[r][-1] for r in range(runs)]}
fig = plt.figure(dpi = 150)
ax = fig.add_subplot(111, projection='3d')

ax.plot_surface(X_grid, Y_grid, V_grid, cmap="cividis", lw=0.5, rstride=1, cstride=1, alpha=0.3)
ax.contour(X_grid, Y_grid, V_grid, 30, linestyles="solid", alpha = 0.3)

xy_sim = np.vstack([coords["X"], coords["Y"]])
density = gaussian_kde(xy_sim)(xy_sim)
ax.scatter(coords["X"], coords["Y"], coords["V"], c = density, cmap = "coolwarm")

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('V')

ax.view_init(30, -140)

plt.show()
#plt.savefig('Lscape_tristab_n10_By4-0_stochTraj.png', bbox_inches='tight')
coords = {"X":[X_runs[r][-1] for r in range(runs)], "Y":[Y_runs[r][-1] for r in range(runs)], "V":[V_runs[r][-1] for r in range(runs)]}
plt.figure(dpi = 150)
plt.title("50 Time Step Distribution")
plt.hist(coords["X"], bins = 100)
plt.xlabel("Log X Expression")
plt.ylabel("Counts")
plt.show()
plt.figure(dpi = 150)
plt.title("50 Time Step Distribution")
plt.hist(coords["Y"], bins = 100)
plt.xlabel("Log Y Expression")
plt.ylabel("Counts")
plt.show()
#Initialize Number of runs
runs = 10000
#Initialize run dictionaries for t, X, Y, and V
t_runs = {}
X_runs = {}
Y_runs = {}
V_runs = {}

#Run simulations in parallel
pool = mp.Pool(mp.cpu_count())
results = []
results = pool.starmap_async(gillespie_sim, [(r, initials, propensities, stoichiometry, 25) for r in range(runs)]).get()
pool.close()
# Collect all the results into dictionaries
for result in results:
    r = result[0]
    t_runs[r] = result[1]
    X_stoch, Y_stoch = zip(*result[2])
    X_runs[r] = X_stoch
    Y_runs[r] = Y_stoch
    V_stoch = []
    for i in range(0, len(X_stoch)):
        x_idx = int(X_stoch[i] / x_grid_step)
        y_idx = int(Y_stoch[i] / y_grid_step)
        V_stoch.append(V_grid[y_idx, x_idx])
    V_runs[r] = V_stoch
coords = {"X":[X_runs[r][-1] for r in range(runs)], "Y":[Y_runs[r][-1] for r in range(runs)], "V":[V_runs[r][-1] for r in range(runs)]}
fig = plt.figure(dpi = 150)
ax = fig.add_subplot(111, projection='3d')

ax.plot_surface(X_grid, Y_grid, V_grid, cmap="cividis", lw=0.5, rstride=1, cstride=1, alpha=0.3)
ax.contour(X_grid, Y_grid, V_grid, 30, linestyles="solid", alpha = 0.3)

xy_sim = np.vstack([coords["X"], coords["Y"]])
density = gaussian_kde(xy_sim)(xy_sim)
ax.scatter(coords["X"], coords["Y"], coords["V"], c = density, cmap = "coolwarm")

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('V')

ax.view_init(30, -140)

plt.show()
#plt.savefig('Lscape_tristab_n10_By4-0_stochTraj.png', bbox_inches='tight')
plt.figure(dpi = 150)
plt.title("25 Time Step Distribution")
plt.hist(coords["X"], bins = 100)
plt.xlabel("Log X Expression")
plt.ylabel("Counts")
plt.show()
plt.figure(dpi = 150)
plt.title("25 Time Step Distribution")
plt.hist(coords["Y"], bins = 100)
plt.xlabel("Log Y Expression")
plt.ylabel("Counts")
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16427822439706397, 'fraction_numerical': 0.028078760340207388, 'mean_word_length': 4.91592005513439, 'pattern_counts': {'":': 9, '<': 32, '<?xml version=': 0, '>': 32, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '100283', 'license': 'MIT', 'n_tokens_mistral': 3896, 'n_tokens_neox': 3651, 'n_words': 800, 'path': 'notebooks/Surface_Mapping_Bistable.ipynb', 'repo_name': 'omiethescientist/Waddington-Potential-Surface', 'size': '8583'}