SHOW:
|
|
- or go back to the newest paste.
| 1 | import numpy as np | |
| 2 | import matplotlib.pyplot as plt | |
| 3 | from matplotlib import cm | |
| 4 | from tqdm import tqdm | |
| 5 | import os | |
| 6 | from scipy.integrate import solve_ivp | |
| 7 | ||
| 8 | # Parameters | |
| 9 | SM = 0.0 | |
| 10 | INT = 5.0 | |
| 11 | PROD = -2.0 | |
| 12 | ||
| 13 | TS1_values = np.linspace(5, 30, 40) | |
| 14 | TS2_values = np.linspace(5, 40, 40) | |
| 15 | ||
| 16 | R = 8.314e-3 | |
| 17 | T = 298.0 | |
| 18 | RT = R * T | |
| 19 | A = 1.0 | |
| 20 | ||
| 21 | # ODE system | |
| 22 | def reaction_odes(t, y, k1f, k1r, k2f, k2r): | |
| 23 | SM_conc, INT_conc, PROD_conc = y | |
| 24 | dSMdt = -k1f * SM_conc + k1r * INT_conc | |
| 25 | dINTdt = k1f * SM_conc - k1r * INT_conc - k2f * INT_conc + k2r * PROD_conc | |
| 26 | dPRODdt = k2f * INT_conc - k2r * PROD_conc | |
| 27 | return [dSMdt, dINTdt, dPRODdt] | |
| 28 | ||
| 29 | y0 = [1.0, 0.0, 0.0] | |
| 30 | t_span = (0, 1000.0) | |
| 31 | t_eval = [1000.0] | |
| 32 | ||
| 33 | TS1_grid, TS2_grid = np.meshgrid(TS1_values, TS2_values) | |
| 34 | rate_grid = np.zeros_like(TS1_grid) | |
| 35 | ||
| 36 | # Calculate rates | |
| 37 | for i in tqdm(range(len(TS2_values)), desc="Rows"): | |
| 38 | for j in range(len(TS1_values)): | |
| 39 | TS1_val = TS1_grid[i, j] | |
| 40 | TS2_val = TS2_grid[i, j] | |
| 41 | ||
| 42 | - | if (TS1_val >= 5) and (TS2_val >= 5) and ((TS2_val - INT) > (TS1_val - SM)): |
| 42 | + | if (TS1_val >= 5) and (TS2_val >= 5) and ((TS2_val - INT) < (TS1_val - SM)): |
| 43 | k1f = A * np.exp(-(TS1_val - SM) / RT) | |
| 44 | k1r = A * np.exp(-(TS1_val - INT) / RT) | |
| 45 | k2f = A * np.exp(-(TS2_val - INT) / RT) | |
| 46 | k2r = A * np.exp(-(TS2_val - PROD) / RT) | |
| 47 | ||
| 48 | sol = solve_ivp(reaction_odes, t_span, y0, t_eval=t_eval, args=(k1f, k1r, k2f, k2r)) | |
| 49 | if sol.success: | |
| 50 | SM_final, INT_final, PROD_final = sol.y[:, -1] | |
| 51 | rate = (k2f * INT_final) - (k2r * PROD_final) | |
| 52 | else: | |
| 53 | rate = np.nan | |
| 54 | else: | |
| 55 | rate = np.nan | |
| 56 | ||
| 57 | rate_grid[i, j] = rate | |
| 58 | ||
| 59 | # Stats | |
| 60 | valid_rates = rate_grid[~np.isnan(rate_grid)] | |
| 61 | mean_rate = np.mean(valid_rates) if valid_rates.size > 0 else np.nan | |
| 62 | std_rate = np.std(valid_rates) if valid_rates.size > 0 else np.nan | |
| 63 | max_rate = np.nanmax(rate_grid) | |
| 64 | min_rate = np.nanmin(rate_grid) | |
| 65 | ||
| 66 | print("Summary Statistics:")
| |
| 67 | print(f"Mean rate: {mean_rate:.4e}")
| |
| 68 | print(f"Std dev of rate: {std_rate:.4e}")
| |
| 69 | print(f"Max rate: {max_rate:.4e}")
| |
| 70 | print(f"Min rate: {min_rate:.4e}")
| |
| 71 | ||
| 72 | # Plot | |
| 73 | fig = plt.figure(figsize=(10, 7)) | |
| 74 | ax = fig.add_subplot(111, projection='3d') | |
| 75 | ||
| 76 | masked_rate = np.ma.masked_invalid(rate_grid) | |
| 77 | surf = ax.plot_surface(TS1_grid, TS2_grid, masked_rate, cmap=cm.viridis, edgecolor='none') | |
| 78 | fig.colorbar(surf, shrink=0.5, aspect=5) | |
| 79 | ||
| 80 | ax.set_xlabel('TS1 Free Energy')
| |
| 81 | ax.set_ylabel('TS2 Free Energy')
| |
| 82 | ax.set_zlabel('Net Rate')
| |
| 83 | ax.set_title('Rate Dependence (Assuming A=1)')
| |
| 84 | ||
| 85 | ax.invert_xaxis() | |
| 86 | ax.invert_yaxis() | |
| 87 | ax.view_init(elev=30, azim=210) | |
| 88 | ||
| 89 | plt.tight_layout() | |
| 90 | ||
| 91 | output_path = "3d_rate_plot_odes_reversible.png" | |
| 92 | plt.savefig(output_path, dpi=300) | |
| 93 | print(f"Plot saved as {os.path.abspath(output_path)}")
| |
| 94 | ||
| 95 | plt.show() | |
| 96 |
