diff --git a/biocircuits/apps/ffl.py b/biocircuits/apps/ffl.py index bcc9066..52bedeb 100644 --- a/biocircuits/apps/ffl.py +++ b/biocircuits/apps/ffl.py @@ -259,7 +259,7 @@ def ffl_app(): import biocircuits.apps import bokeh.plotting - app = biocircuits.apps.promiscuous_222_app() + app = biocircuits.apps.ffl_app() app(bokeh.plotting.curdoc()) ``` diff --git a/biocircuits/jsfunctions.py b/biocircuits/jsfunctions.py index 0295770..03ed70b 100644 --- a/biocircuits/jsfunctions.py +++ b/biocircuits/jsfunctions.py @@ -1248,6 +1248,68 @@ return LUPSolve(LU, p, b); } +""", + "michaelis_menten_approx": """ +function michaelisMentenRHS(c, t, kappa, zeta){ + let [cs, ces, cp] = c; + + return [ + (-(1.0 - ces) * cs + (1.0 - kappa) * ces) / kappa, + ((1.0 - ces) * cs - ces) / kappa / zeta, + ces + ]; +} + + +function approxMichaelisMentenRHS(c, t, kappa, zeta) { + let cs = c[0]; + return [-cs / (1.0 + cs)]; +} + + +function callback() { + let xRangeMax = xRange.end; + let dt = 0.01; + let cs0 = Math.pow(10.0, cs0Slider.value); + let cp0 = 0.0; // No product to start + let ces0 = 0.0; + let c0 = [cs0, 0.0, 0.0]; + let c0Approx = [cs0]; + let kappa = kappaSlider.value; + let zeta = Math.pow(10.0, zetaSlider.value); + + let t = linspace(0.0, xRangeMax, cds.data['t'].length); + let args = [kappa, zeta]; + + // Integrate ODES + let cSolve = rkf45(michaelisMentenRHS, c0, t, args, dt); + let csApprox = rkf45(approxMichaelisMentenRHS, c0Approx, t, args)[0]; + + // Compute product and enzyme conc from substrate in approximate regime + let cesApprox = []; + let cpApprox = []; + for (let i = 0; i < csApprox.length; i++) { + cesApprox[i] = csApprox[i] / (1.0 + csApprox[i]); + cpApprox[i] = cs0 + cp0 - csApprox[i] + zeta * (ces0 - cesApprox[i]); + } + + // Update data + cds.data['t'] = t; + cds.data['cs'] = cSolve[0]; + cds.data['ces'] = cSolve[1]; + cds.data['cp'] = cSolve[2]; + cds.data['cs_approx'] = csApprox; + cds.data['ces_approx'] = cesApprox; + cds.data['cp_approx'] = cpApprox; + + // Update y-range + yRange.start = -0.02 * cs0; + yRange.end = 1.02 * cs0; + + cds.change.emit(); + yRange.change.emit(); +} + """, "autorepressor_response_to_pulse": """ function negAutoRHS(x, t, beta0, gamma, k, n, ks, ns, sFun, sArgs=[]) { @@ -1476,6 +1538,84 @@ cds_fp_unstable.change.emit(); } +""", + "proteinRepressilator": """ +function proteinRepressilator(x, t, beta, n) { + // Unpack + var [x1, x2, x3] = x; + + return [ + beta * rep_hill(x3, n) - x1, + beta * rep_hill(x1, n) - x2, + beta * rep_hill(x2, n) - x3 + ]; +} + + +function callback() { + let xRangeMax = xRange.end; + let dt = 0.01; + let x0 = [1.0, 1.0, 1.2]; + let beta = betaSlider.value; + let n = nSlider.value; + + let t = linspace(0.0, xRangeMax, cds.data['t'].length); + let args = [beta, n]; + + // Integrate ODES + let xSolve = rkf45(proteinRepressilator, x0, t, args, t[1] - t[0], 1e-7, 1e-3, 100); + + cds.data['t'] = t; + cds.data['x1'] = xSolve[0]; + cds.data['x2'] = xSolve[1]; + cds.data['x3'] = xSolve[2]; + + cds.change.emit(); +} + +""", + "repressilator": """ +function repressilator(x, t, beta, rho, gamma, n) { + // Unpack + let [m1, m2, m3, x1, x2, x3] = x; + + return [ + beta * (rho + rep_hill(x3, n)) - m1, + beta * (rho + rep_hill(x1, n)) - m2, + beta * (rho + rep_hill(x2, n)) - m3, + gamma * (m1 - x1), + gamma * (m2 - x2), + gamma * (m3 - x3) + ]; +} + + +function callback() { + let xRangeMax = xRange.end; + let dt = 0.01; + let x0 = [0.0, 0.0, 0.0, 1.0, 1.0, 1.2]; + let beta = Math.pow(10, betaSlider.value); + let gamma = Math.pow(10, gammaSlider.value); + let rho = Math.pow(10, rhoSlider.value); + let n = nSlider.value; + + let t = linspace(0.0, xRangeMax, cds.data['t'].length); + let args = [beta, rho, gamma, n]; + + // Integrate ODES + let xSolve = rkf45(repressilator, x0, t, args, t[1] - t[0], 1e-7, 1e-3, 100); + + cds.data['t'] = t; + cds.data['m1'] = xSolve[0]; + cds.data['m2'] = xSolve[1]; + cds.data['m3'] = xSolve[2]; + cds.data['x1'] = xSolve[3]; + cds.data['x2'] = xSolve[4]; + cds.data['x3'] = xSolve[5]; + + cds.change.emit(); +} + """, "toggle_nullclines": """ function f(x, beta, n) {return beta / (1.0 + Math.pow(x, n));} @@ -1618,4 +1758,4 @@ cdsUnstable.change.emit(); } """, -} +} \ No newline at end of file diff --git a/biocircuits/jsplots.py b/biocircuits/jsplots.py index 0eeff15..f2462aa 100644 --- a/biocircuits/jsplots.py +++ b/biocircuits/jsplots.py @@ -2,6 +2,7 @@ import numpy as np import scipy.integrate +import scipy.special import matplotlib.streamplot @@ -57,8 +58,167 @@ def _sin_plot(): return layout +def michaelis_menten_approx(): + def michaelis_menten_rhs(c, t, kappa, zeta): + cs, ces, cp = c + return np.array( + [ + (-(1 - ces) * cs + (1 - kappa) * ces) / kappa, + ((1 - ces) * cs - ces) / kappa / zeta, + ces, + ] + ) + + def approx_michaelis_menten(c0, t, kappa, zeta): + """Analytical solution to the Michaelis-Menten equation.""" + cs0, ces0, cp0 = c0 + cs = scipy.special.lambertw(cs0 * np.exp(cs0 - t)).real + ces = cs / (1 + cs) + cp = cs0 + cp0 - cs - zeta * (ces0 + ces) + + return cs, ces, cp + + kappa_slider = bokeh.models.Slider( + title="κ", start=0.01, end=1, value=0.5, step=0.01, width=100 + ) + zeta_slider = bokeh.models.Slider( + title="ζ", + start=-2, + end=2, + value=-2, + step=0.05, + width=100, + format=bokeh.models.FuncTickFormatter( + code="return Math.pow(10, tick).toFixed(2)" + ), + ) + cs0_slider = bokeh.models.Slider( + title="init. substr. conc.", + start=-1, + end=1, + value=0.0, + step=0.01, + width=100, + format=bokeh.models.FuncTickFormatter( + code="return Math.pow(10, tick).toFixed(2)" + ), + ) + + def solve_mm(kappa, zeta, cs0): + # Initial condition + c0 = np.array([cs0, 0.0, 0.0]) + + # Time points + t = np.linspace(0, 10, 400) + + # Solve the full system + c = scipy.integrate.odeint(michaelis_menten_rhs, c0, t, args=(kappa, zeta)) + cs, ces, cp = c.transpose() + + # Solve the approximate system + cs_approx, ces_approx, cp_approx = approx_michaelis_menten(c0, t, kappa, zeta) + + return t, cs, ces, cp, cs_approx, ces_approx, cp_approx + + # Get solution for initial glyphs + t, cs, ces, cp, cs_approx, ces_approx, cp_approx = solve_mm( + kappa_slider.value, 10 ** zeta_slider.value, 10 ** cs0_slider.value + ) + + # Set up ColumnDataSource for plot + cds = bokeh.models.ColumnDataSource( + dict( + t=t, + cs=cs, + ces=ces, + cp=cp, + cs_approx=cs_approx, + ces_approx=ces_approx, + cp_approx=cp_approx, + ) + ) + + # Make the plot + p = bokeh.plotting.figure( + plot_width=500, + plot_height=250, + x_axis_label="dimensionless time", + y_axis_label="dimensionless concentration", + x_range=[0, 10], + y_range=[-0.02, 1.02], + ) + + colors = colorcet.b_glasbey_category10 + + # Populate glyphs + p.line( + source=cds, x="t", y="ces", line_width=2, color=colors[0], legend_label="ES", + ) + p.line( + source=cds, x="t", y="cs", line_width=2, color=colors[1], legend_label="S", + ) + p.line( + source=cds, x="t", y="cp", line_width=2, color=colors[2], legend_label="P", + ) + p.line( + source=cds, x="t", y="ces_approx", line_width=4, color=colors[0], alpha=0.3, + ) + p.line( + source=cds, x="t", y="cs_approx", line_width=4, color=colors[1], alpha=0.3, + ) + p.line( + source=cds, x="t", y="cp_approx", line_width=4, color=colors[2], alpha=0.3, + ) + + p.legend.location = "center_right" + + # JavaScript callback + js_code = ( + jsfuns["michaelis_menten_approx"] + + jsfuns["utils"] + + jsfuns["linalg"] + + jsfuns["ode"] + + "callback()" + ) + callback = bokeh.models.CustomJS( + args=dict( + cds=cds, + kappaSlider=kappa_slider, + zetaSlider=zeta_slider, + cs0Slider=cs0_slider, + xRange=p.x_range, + yRange=p.y_range, + ), + code=js_code, + ) + + # Link sliders + kappa_slider.js_on_change("value", callback) + zeta_slider.js_on_change("value", callback) + cs0_slider.js_on_change("value", callback) + + # Also trigger if x_range changes + p.x_range.js_on_change("end", callback) + + # Build layout + layout = bokeh.layouts.column( + bokeh.layouts.row( + bokeh.models.Spacer(width=40), + kappa_slider, + bokeh.models.Spacer(width=10), + zeta_slider, + bokeh.models.Spacer(width=10), + cs0_slider, + ), + bokeh.models.Spacer(width=10), + p, + ) + + return layout + + def gaussian_pulse(): - """Make a plot of a Gaussian pulse/ + """Make a plot of a Gaussian pulse. """ # t/s data for plotting t_0 = 4.0 @@ -514,14 +674,14 @@ def toggle_nullclines(): ) -def repressilator(): - """Replaces the plot of the protein-only repressilator. Replaces - Python code: +def protein_repressilator(): + """Plot the dynamics of a protein-only repressilator circuit. + """ - def repressilator_rhs(x, t, beta, n): - ''' + def protein_repressilator_rhs(x, t, beta, n): + """ Returns 3-array of (dx_1/dt, dx_2/dt, dx_3/dt) - ''' + """ x_1, x_2, x_3 = x return np.array( @@ -532,7 +692,6 @@ def repressilator_rhs(x, t, beta, n): ] ) - # Initial condiations x0 = np.array([1, 1, 1.2]) @@ -540,26 +699,29 @@ def repressilator_rhs(x, t, beta, n): n_points = 1000 # Widgets for controlling parameters - beta_slider = bokeh.models.Slider(title="β", start=0, end=100, step=0.1, value=10) + beta_slider = bokeh.models.Slider( + title="β", start=0.01, end=100, step=0.01, value=10.0 + ) n_slider = bokeh.models.Slider(title="n", start=1, end=5, step=0.1, value=3) - t_max_slider = bokeh.models.Slider(title="t_max", start=1, end=100, step=1, value=40) # Solve for species concentrations def _solve_repressilator(beta, n, t_max): t = np.linspace(0, t_max, n_points) - x = scipy.integrate.odeint(repressilator_rhs, x0, t, args=(beta, n)) + x = scipy.integrate.odeint(protein_repressilator_rhs, x0, t, args=(beta, n)) return t, x.transpose() - # Obtain solution for plot - t, x = _solve_repressilator(beta_slider.value, n_slider.value, t_max_slider.value) + t, x = _solve_repressilator(beta_slider.value, n_slider.value, 40.0) # Build the plot colors = colorcet.b_glasbey_category10[:3] p_rep = bokeh.plotting.figure( - frame_width=550, frame_height=200, x_axis_label="t", x_range=[0, t_max_slider.value] + frame_width=550, + frame_height=200, + x_axis_label="t", + x_range=[0, 40.0], ) cds = bokeh.models.ColumnDataSource(data=dict(t=t, x1=x[0], x2=x[1], x3=x[2])) @@ -576,7 +738,6 @@ def _solve_repressilator(beta, n, t_max): p_rep.legend.location = "top_left" - # Set up plot p_phase = bokeh.plotting.figure( frame_width=200, frame_height=200, x_axis_label="x₁", y_axis_label="x₂", @@ -584,194 +745,190 @@ def _solve_repressilator(beta, n, t_max): p_phase.line(source=cds, x="x1", y="x2", line_width=2) + # Set up callbacks + js_code = ( + jsfuns["reg"] + + jsfuns["ode"] + + jsfuns["circuits"] + + jsfuns["utils"] + + jsfuns["linalg"] + + jsfuns["proteinRepressilator"] + + 'callback()' + ) + callback = bokeh.models.CustomJS( + args=dict( + cds=cds, + xRange=p_rep.x_range, + betaSlider=beta_slider, + nSlider=n_slider, + ), + code=js_code, + ) - if fully_interactive_plots: - # Set up callbacks - def _callback(attr, old, new): - t, x = _solve_repressilator( - beta_slider.value, n_slider.value, t_max_slider.value - ) - cds.data = dict(t=t, x1=x[0], x2=x[1], x3=x[2]) - p_rep.x_range.end = t_max_slider.value - - beta_slider.on_change("value", _callback) - n_slider.on_change("value", _callback) - t_max_slider.on_change("value", _callback) - - # Build layout - repressilator_layout = bokeh.layouts.column( - p_rep, - bokeh.layouts.Spacer(height=10), - bokeh.layouts.row( - p_phase, - bokeh.layouts.Spacer(width=70), - bokeh.layouts.column(beta_slider, n_slider, t_max_slider, width=150), - ), - ) + beta_slider.js_on_change("value", callback) + n_slider.js_on_change("value", callback) + p_rep.x_range.js_on_change("end", callback) - # Build the app - def repressilator_app(doc): - doc.add_root(repressilator_layout) - - bokeh.io.show(repressilator_app, notebook_url=notebook_url) - else: - beta_slider.disabled = True - n_slider.disabled = True - t_max_slider.disabled = True - - # Build layout - repressilator_layout = bokeh.layouts.column( - p_rep, - bokeh.layouts.Spacer(height=10), - bokeh.layouts.row( - p_phase, - bokeh.layouts.Spacer(width=70), - bokeh.layouts.column( - bokeh.layouts.column(beta_slider, n_slider, t_max_slider, width=150), - bokeh.models.Div( - text=''' -

Sliders are inactive. To get active sliders, re-run notebook with - fully_interactive_plots = True - in the first code cell.

- ''', - width=250, - ), - ), - ), - ) + # Build layout + layout = bokeh.layouts.column( + p_rep, + bokeh.layouts.Spacer(height=10), + bokeh.layouts.row( + p_phase, + bokeh.layouts.Spacer(width=70), + bokeh.layouts.column(beta_slider, n_slider, width=150), + ), + ) + + return layout - bokeh.io.show(repressilator_layout) +def repressilator(): + """Plot the dynamics of a repressilator circuit. """ - def repressilator_rhs(x, t, beta, n): + # Sliders + beta_slider = bokeh.models.Slider( + title="β", + start=0, + end=4, + step=0.1, + value=1, + format=bokeh.models.FuncTickFormatter(code="return Math.pow(10, tick).toFixed(2)"), + ) + gamma_slider = bokeh.models.Slider( + title="γ", + start=-3, + end=0, + step=0.1, + value=0, + format=bokeh.models.FuncTickFormatter(code="return Math.pow(10, tick).toFixed(3)"), + ) + rho_slider = bokeh.models.Slider( + title="ρ", + start=-6, + end=0, + step=0.1, + value=-3, + format=bokeh.models.FuncTickFormatter(code="return Math.pow(10, tick).toFixed(6)"), + ) + n_slider = bokeh.models.Slider(title="n", start=1, end=5, step=0.1, value=3) + + + def repressilator_rhs(mx, t, beta, gamma, rho, n): """ - Returns 3-array of (dx_1/dt, dx_2/dt, dx_3/dt) + Returns 6-array of (dm_1/dt, dm_2/dt, dm_3/dt, dx_1/dt, dx_2/dt, dx_3/dt) """ - x_1, x_2, x_3 = x - + m_1, m_2, m_3, x_1, x_2, x_3 = mx return np.array( [ - beta / (1 + x_3 ** n) - x_1, - beta / (1 + x_1 ** n) - x_2, - beta / (1 + x_2 ** n) - x_3, + beta * (rho + 1 / (1 + x_3 ** n)) - m_1, + beta * (rho + 1 / (1 + x_1 ** n)) - m_2, + beta * (rho + 1 / (1 + x_2 ** n)) - m_3, + gamma * (m_1 - x_1), + gamma * (m_2 - x_2), + gamma * (m_3 - x_3), ] ) + # Initial condiations - x0 = np.array([1, 1, 1.2]) + x0 = np.array([0, 0, 0, 1, 1.1, 1.2]) # Number of points to use in plots n_points = 1000 - # Widgets for controlling parameters - beta_slider = bokeh.models.Slider( - title="β", start=0.01, end=100, step=0.01, value=10.0 - ) - n_slider = bokeh.models.Slider(title="n", start=1, end=5, step=0.1, value=3) - t_max_slider = bokeh.models.Slider( - title="t_max", start=1, end=100, step=1, value=40 - ) - # Solve for species concentrations - def _solve_repressilator(beta, n, t_max): + def _solve_repressilator(log_beta, log_gamma, log_rho, n, t_max): + beta = 10 ** log_beta + gamma = 10 ** log_gamma + rho = 10 ** log_rho t = np.linspace(0, t_max, n_points) - x = scipy.integrate.odeint(repressilator_rhs, x0, t, args=(beta, n)) + x = scipy.integrate.odeint(repressilator_rhs, x0, t, args=(beta, gamma, rho, n)) + m1, m2, m3, x1, x2, x3 = x.transpose() + return t, m1, m2, m3, x1, x2, x3 - return t, x.transpose() - # Obtain solution for plot - t, x = _solve_repressilator(beta_slider.value, n_slider.value, t_max_slider.value) + t, m1, m2, m3, x1, x2, x3 = _solve_repressilator( + beta_slider.value, + gamma_slider.value, + rho_slider.value, + n_slider.value, + 40.0, + ) - # Build the plot - colors = colorcet.b_glasbey_category10[:3] + cds = bokeh.models.ColumnDataSource( + dict(t=t, m1=m1, m2=m2, m3=m3, x1=x1, x2=x2, x3=x3) + ) - p_rep = bokeh.plotting.figure( - frame_width=550, + p = bokeh.plotting.figure( + frame_width=500, frame_height=200, x_axis_label="t", - x_range=[0, t_max_slider.value], + x_range=[0, 40.0], ) - cds = bokeh.models.ColumnDataSource(data=dict(t=t, x1=x[0], x2=x[1], x3=x[2])) - labels = dict(x1="x₁", x2="x₂", x3="x₃") - for color, x_val in zip(colors, labels): - p_rep.line( - source=cds, - x="t", - y=x_val, - color=color, - legend_label=labels[x_val], - line_width=2, - ) + colors = bokeh.palettes.d3["Category20"][6] + m1_line = p.line(source=cds, x="t", y="m1", line_width=2, color=colors[1]) + x1_line = p.line(source=cds, x="t", y="x1", line_width=2, color=colors[0]) + m2_line = p.line(source=cds, x="t", y="m2", line_width=2, color=colors[3]) + x2_line = p.line(source=cds, x="t", y="x2", line_width=2, color=colors[2]) + m3_line = p.line(source=cds, x="t", y="m3", line_width=2, color=colors[5]) + x3_line = p.line(source=cds, x="t", y="x3", line_width=2, color=colors[4]) + + legend_items = [ + ("m₁", [m1_line]), + ("x₁", [x1_line]), + ("m₂", [m2_line]), + ("x₂", [x2_line]), + ("m₃", [m3_line]), + ("x₃", [x3_line]), + ] + legend = bokeh.models.Legend(items=legend_items) + legend.click_policy = 'hide' - p_rep.legend.location = "top_left" + p.add_layout(legend, "right") - # Set up plot - p_phase = bokeh.plotting.figure( - frame_width=200, frame_height=200, x_axis_label="x₁", y_axis_label="x₂", + # Build the layout + layout = bokeh.layouts.column( + bokeh.layouts.row( + beta_slider, + gamma_slider, + rho_slider, + n_slider, + width=575, + ), + bokeh.layouts.Spacer(height=10), + p, ) - p_phase.line(source=cds, x="x1", y="x2", line_width=2) - # Set up callbacks js_code = ( jsfuns["reg"] + jsfuns["ode"] + jsfuns["circuits"] - + """ -var t = cds.data['t']; - -var beta = beta_slider.value; -var n = n_slider.value; -var t_max = t_max_slider.value; - -t = linspace(0.0, t_max, t.length); - -var x = rkf45(repressilator, [1.0, 1.0, 1.2], t, [beta, n], t[1] - t[0], 1e-7, 1e-3, 100); - -cds.data['t'] = t; -cds.data['x1'] = x[0]; -cds.data['x2'] = x[1]; -cds.data['x3'] = x[2]; - -x_range.end = t_max; - -cds.change.emit(); -""" + + jsfuns["utils"] + + jsfuns["linalg"] + + jsfuns["repressilator"] + + 'callback()' ) callback = bokeh.models.CustomJS( args=dict( cds=cds, - x_range=p_rep.x_range, - beta_slider=beta_slider, - n_slider=n_slider, - t_max_slider=t_max_slider, + xRange=p.x_range, + betaSlider=beta_slider, + rhoSlider=rho_slider, + gammaSlider=gamma_slider, + nSlider=n_slider, ), code=js_code, ) - def _callback(attr, old, new): - t, x = _solve_repressilator( - beta_slider.value, n_slider.value, t_max_slider.value - ) - cds.data = dict(t=t, x1=x[0], x2=x[1], x3=x[2]) - p_rep.x_range.end = t_max_slider.value - beta_slider.js_on_change("value", callback) + gamma_slider.js_on_change("value", callback) + rho_slider.js_on_change("value", callback) n_slider.js_on_change("value", callback) - t_max_slider.js_on_change("value", callback) - - # Build layout - layout = bokeh.layouts.column( - p_rep, - bokeh.layouts.Spacer(height=10), - bokeh.layouts.row( - p_phase, - bokeh.layouts.Spacer(width=70), - bokeh.layouts.column(beta_slider, n_slider, t_max_slider, width=150), - ), - ) + p.x_range.js_on_change("end", callback) return layout diff --git a/js/plot_specific/michaelis_menten_approx.js b/js/plot_specific/michaelis_menten_approx.js new file mode 100644 index 0000000..9ed7534 --- /dev/null +++ b/js/plot_specific/michaelis_menten_approx.js @@ -0,0 +1,59 @@ +function michaelisMentenRHS(c, t, kappa, zeta){ + let [cs, ces, cp] = c; + + return [ + (-(1.0 - ces) * cs + (1.0 - kappa) * ces) / kappa, + ((1.0 - ces) * cs - ces) / kappa / zeta, + ces + ]; +} + + +function approxMichaelisMentenRHS(c, t, kappa, zeta) { + let cs = c[0]; + return [-cs / (1.0 + cs)]; +} + + +function callback() { + let xRangeMax = xRange.end; + let dt = 0.01; + let cs0 = Math.pow(10.0, cs0Slider.value); + let cp0 = 0.0; // No product to start + let ces0 = 0.0; + let c0 = [cs0, 0.0, 0.0]; + let c0Approx = [cs0]; + let kappa = kappaSlider.value; + let zeta = Math.pow(10.0, zetaSlider.value); + + let t = linspace(0.0, xRangeMax, cds.data['t'].length); + let args = [kappa, zeta]; + + // Integrate ODES + let cSolve = rkf45(michaelisMentenRHS, c0, t, args, dt); + let csApprox = rkf45(approxMichaelisMentenRHS, c0Approx, t, args)[0]; + + // Compute product and enzyme conc from substrate in approximate regime + let cesApprox = []; + let cpApprox = []; + for (let i = 0; i < csApprox.length; i++) { + cesApprox[i] = csApprox[i] / (1.0 + csApprox[i]); + cpApprox[i] = cs0 + cp0 - csApprox[i] + zeta * (ces0 - cesApprox[i]); + } + + // Update data + cds.data['t'] = t; + cds.data['cs'] = cSolve[0]; + cds.data['ces'] = cSolve[1]; + cds.data['cp'] = cSolve[2]; + cds.data['cs_approx'] = csApprox; + cds.data['ces_approx'] = cesApprox; + cds.data['cp_approx'] = cpApprox; + + // Update y-range + yRange.start = -0.02 * cs0; + yRange.end = 1.02 * cs0; + + cds.change.emit(); + yRange.change.emit(); +} diff --git a/js/plot_specific/proteinRepressilator.js b/js/plot_specific/proteinRepressilator.js new file mode 100644 index 0000000..8a9f154 --- /dev/null +++ b/js/plot_specific/proteinRepressilator.js @@ -0,0 +1,32 @@ +function proteinRepressilator(x, t, beta, n) { + // Unpack + var [x1, x2, x3] = x; + + return [ + beta * rep_hill(x3, n) - x1, + beta * rep_hill(x1, n) - x2, + beta * rep_hill(x2, n) - x3 + ]; +} + + +function callback() { + let xRangeMax = xRange.end; + let dt = 0.01; + let x0 = [1.0, 1.0, 1.2]; + let beta = betaSlider.value; + let n = nSlider.value; + + let t = linspace(0.0, xRangeMax, cds.data['t'].length); + let args = [beta, n]; + + // Integrate ODES + let xSolve = rkf45(proteinRepressilator, x0, t, args, t[1] - t[0], 1e-7, 1e-3, 100); + + cds.data['t'] = t; + cds.data['x1'] = xSolve[0]; + cds.data['x2'] = xSolve[1]; + cds.data['x3'] = xSolve[2]; + + cds.change.emit(); +} diff --git a/js/plot_specific/repressilator.js b/js/plot_specific/repressilator.js new file mode 100644 index 0000000..de93a84 --- /dev/null +++ b/js/plot_specific/repressilator.js @@ -0,0 +1,40 @@ +function repressilator(x, t, beta, rho, gamma, n) { + // Unpack + let [m1, m2, m3, x1, x2, x3] = x; + + return [ + beta * (rho + rep_hill(x3, n)) - m1, + beta * (rho + rep_hill(x1, n)) - m2, + beta * (rho + rep_hill(x2, n)) - m3, + gamma * (m1 - x1), + gamma * (m2 - x2), + gamma * (m3 - x3) + ]; +} + + +function callback() { + let xRangeMax = xRange.end; + let dt = 0.01; + let x0 = [0.0, 0.0, 0.0, 1.0, 1.0, 1.2]; + let beta = Math.pow(10, betaSlider.value); + let gamma = Math.pow(10, gammaSlider.value); + let rho = Math.pow(10, rhoSlider.value); + let n = nSlider.value; + + let t = linspace(0.0, xRangeMax, cds.data['t'].length); + let args = [beta, rho, gamma, n]; + + // Integrate ODES + let xSolve = rkf45(repressilator, x0, t, args, t[1] - t[0], 1e-7, 1e-3, 100); + + cds.data['t'] = t; + cds.data['m1'] = xSolve[0]; + cds.data['m2'] = xSolve[1]; + cds.data['m3'] = xSolve[2]; + cds.data['x1'] = xSolve[3]; + cds.data['x2'] = xSolve[4]; + cds.data['x3'] = xSolve[5]; + + cds.change.emit(); +}