|
8 | 8 | import os |
9 | 9 |
|
10 | 10 | import numpy as np |
11 | | -import matplotlib.pyplot as mpl |
| 11 | +import matplotlib.pyplot as plt |
12 | 12 | from scipy.integrate import odeint |
13 | 13 | from control import phase_plot, box_grid |
14 | 14 |
|
15 | 15 | # Simple model of a genetic switch |
16 | | -# |
| 16 | + |
17 | 17 | # This function implements the basic model of the genetic switch |
18 | 18 | # Parameters taken from Gardner, Cantor and Collins, Nature, 2000 |
19 | 19 | def genswitch(y, t, mu=4, n=2): |
20 | | - return (mu / (1 + y[1]**n) - y[0], mu / (1 + y[0]**n) - y[1]) |
| 20 | + return mu/(1 + y[1]**n) - y[0], mu/(1 + y[0]**n) - y[1] |
21 | 21 |
|
22 | 22 | # Run a simulation from an initial condition |
23 | 23 | tim1 = np.linspace(0, 10, 100) |
24 | 24 | sol1 = odeint(genswitch, [1, 5], tim1) |
25 | 25 |
|
26 | | -# Extract the equlibirum points |
27 | | -mu = 4; n = 2; # switch parameters |
28 | | -eqpt = np.empty(3); |
29 | | -eqpt[0] = sol1[0,-1] |
30 | | -eqpt[1] = sol1[1,-1] |
31 | | -eqpt[2] = 0; # fzero(@(x) mu/(1+x^2) - x, 2); |
| 26 | +# Extract the equilibrium points |
| 27 | +mu = 4; n = 2 # switch parameters |
| 28 | +eqpt = np.empty(3) |
| 29 | +eqpt[0] = sol1[0, -1] |
| 30 | +eqpt[1] = sol1[1, -1] |
| 31 | +eqpt[2] = 0 # fzero(@(x) mu/(1+x^2) - x, 2) |
32 | 32 |
|
33 | 33 | # Run another simulation showing switching behavior |
34 | | -tim2 = np.linspace(11, 25, 100); |
35 | | -sol2 = odeint(genswitch, sol1[-1,:] + [2, -2], tim2) |
| 34 | +tim2 = np.linspace(11, 25, 100) |
| 35 | +sol2 = odeint(genswitch, sol1[-1, :] + [2, -2], tim2) |
36 | 36 |
|
37 | 37 | # First plot out the curves that define the equilibria |
38 | 38 | u = np.linspace(0, 4.5, 46) |
39 | | -f = np.divide(mu, (1 + u**n)) # mu / (1 + u^n), elementwise |
| 39 | +f = np.divide(mu, (1 + u**n)) # mu/(1 + u^n), element-wise |
40 | 40 |
|
41 | | -mpl.figure(1); mpl.clf(); |
42 | | -mpl.axis([0, 5, 0, 5]); # box on; |
43 | | -mpl.plot(u, f, '-', f, u, '--') # 'LineWidth', AM_data_linewidth); |
44 | | -mpl.legend(('z1, f(z1)', 'z2, f(z2)')) # legend(lgh, 'boxoff'); |
45 | | -mpl.plot([0, 3], [0, 3], 'k-') # 'LineWidth', AM_ref_linewidth); |
46 | | -mpl.plot(eqpt[0], eqpt[1], 'k.', eqpt[1], eqpt[0], 'k.', |
47 | | - eqpt[2], eqpt[2], 'k.') # 'MarkerSize', AM_data_markersize*3); |
48 | | -mpl.xlabel('z1, f(z2)'); |
49 | | -mpl.ylabel('z2, f(z1)'); |
| 41 | +plt.figure(1); plt.clf() |
| 42 | +plt.axis([0, 5, 0, 5]) # box on; |
| 43 | +plt.plot(u, f, '-', f, u, '--') # 'LineWidth', AM_data_linewidth) |
| 44 | +plt.legend(('z1, f(z1)', 'z2, f(z2)')) # legend(lgh, 'boxoff') |
| 45 | +plt.plot([0, 3], [0, 3], 'k-') # 'LineWidth', AM_ref_linewidth) |
| 46 | +plt.plot(eqpt[0], eqpt[1], 'k.', eqpt[1], eqpt[0], 'k.', |
| 47 | + eqpt[2], eqpt[2], 'k.') # 'MarkerSize', AM_data_markersize*3) |
| 48 | +plt.xlabel('z1, f(z2)') |
| 49 | +plt.ylabel('z2, f(z1)') |
50 | 50 |
|
51 | 51 | # Time traces |
52 | | -mpl.figure(3); mpl.clf(); # subplot(221); |
53 | | -mpl.plot(tim1, sol1[:,0], 'b-', tim1, sol1[:,1], 'g--'); |
54 | | -# set(pl, 'LineWidth', AM_data_linewidth); |
55 | | -mpl.plot([tim1[-1], tim1[-1]+1], |
56 | | - [sol1[-1,0], sol2[0,1]], 'ko:', |
57 | | - [tim1[-1], tim1[-1]+1], [sol1[-1,1], sol2[0,0]], 'ko:'); |
58 | | -# set(pl, 'LineWidth', AM_data_linewidth, 'MarkerSize', AM_data_markersize); |
59 | | -mpl.plot(tim2, sol2[:,0], 'b-', tim2, sol2[:,1], 'g--'); |
60 | | -# set(pl, 'LineWidth', AM_data_linewidth); |
61 | | -mpl.axis([0, 25, 0, 5]); |
| 52 | +plt.figure(3); plt.clf() # subplot(221) |
| 53 | +plt.plot(tim1, sol1[:, 0], 'b-', tim1, sol1[:, 1], 'g--') |
| 54 | +# set(pl, 'LineWidth', AM_data_linewidth) |
| 55 | +plt.plot([tim1[-1], tim1[-1] + 1], |
| 56 | + [sol1[-1, 0], sol2[0, 1]], 'ko:', |
| 57 | + [tim1[-1], tim1[-1] + 1], [sol1[-1, 1], sol2[0, 0]], 'ko:') |
| 58 | +# set(pl, 'LineWidth', AM_data_linewidth, 'MarkerSize', AM_data_markersize) |
| 59 | +plt.plot(tim2, sol2[:, 0], 'b-', tim2, sol2[:, 1], 'g--') |
| 60 | +# set(pl, 'LineWidth', AM_data_linewidth) |
| 61 | +plt.axis([0, 25, 0, 5]) |
62 | 62 |
|
63 | | -mpl.xlabel('Time {\itt} [scaled]'); |
64 | | -mpl.ylabel('Protein concentrations [scaled]'); |
65 | | -mpl.legend(('z1 (A)', 'z2 (B)')) # 'Orientation', 'horizontal'); |
66 | | -# legend(legh, 'boxoff'); |
| 63 | +plt.xlabel('Time {\itt} [scaled]') |
| 64 | +plt.ylabel('Protein concentrations [scaled]') |
| 65 | +plt.legend(('z1 (A)', 'z2 (B)')) # 'Orientation', 'horizontal') |
| 66 | +# legend(legh, 'boxoff') |
67 | 67 |
|
68 | 68 | # Phase portrait |
69 | | -mpl.figure(2); mpl.clf(); # subplot(221); |
70 | | -mpl.axis([0, 5, 0, 5]); # set(gca, 'DataAspectRatio', [1, 1, 1]); |
71 | | -phase_plot(genswitch, X0 = box_grid([0, 5, 6], [0, 5, 6]), T = 10, |
72 | | - timepts = [0.2, 0.6, 1.2]) |
| 69 | +plt.figure(2) |
| 70 | +plt.clf() # subplot(221) |
| 71 | +plt.axis([0, 5, 0, 5]) # set(gca, 'DataAspectRatio', [1, 1, 1]) |
| 72 | +phase_plot(genswitch, X0=box_grid([0, 5, 6], [0, 5, 6]), T=10, |
| 73 | + timepts=[0.2, 0.6, 1.2]) |
73 | 74 |
|
74 | 75 | # Add the stable equilibrium points |
75 | | -mpl.plot(eqpt[0], eqpt[1], 'k.', eqpt[1], eqpt[0], 'k.', |
76 | | - eqpt[2], eqpt[2], 'k.') # 'MarkerSize', AM_data_markersize*3); |
| 76 | +plt.plot(eqpt[0], eqpt[1], 'k.', eqpt[1], eqpt[0], 'k.', |
| 77 | + eqpt[2], eqpt[2], 'k.') # 'MarkerSize', AM_data_markersize*3) |
77 | 78 |
|
78 | | -mpl.xlabel('Protein A [scaled]'); |
79 | | -mpl.ylabel('Protein B [scaled]'); # 'Rotation', 90); |
| 79 | +plt.xlabel('Protein A [scaled]') |
| 80 | +plt.ylabel('Protein B [scaled]') # 'Rotation', 90) |
80 | 81 |
|
81 | 82 | if 'PYCONTROL_TEST_EXAMPLES' not in os.environ: |
82 | | - mpl.show() |
| 83 | + plt.show() |
0 commit comments