diff --git a/.gitignore b/.gitignore
old mode 100644
new mode 100755
diff --git a/LICENSE.md b/LICENSE.md
old mode 100644
new mode 100755
diff --git a/MANIFEST.in b/MANIFEST.in
old mode 100644
new mode 100755
diff --git a/README.md b/README.md
old mode 100644
new mode 100755
diff --git a/biocircuits/.gitignore b/biocircuits/.gitignore
old mode 100644
new mode 100755
diff --git a/biocircuits/__init__.py b/biocircuits/__init__.py
old mode 100644
new mode 100755
index 0522dfb..fbfa916
--- a/biocircuits/__init__.py
+++ b/biocircuits/__init__.py
@@ -15,4 +15,4 @@
 
 __author__ = """Justin Bois"""
 __email__ = "bois@caltech.edu"
-__version__ = "0.1.8"
+__version__ = "0.1.9"
diff --git a/biocircuits/apps/.gitignore b/biocircuits/apps/.gitignore
old mode 100644
new mode 100755
diff --git a/biocircuits/apps/__init__.py b/biocircuits/apps/__init__.py
old mode 100644
new mode 100755
diff --git a/biocircuits/apps/ffl.py b/biocircuits/apps/ffl.py
old mode 100644
new mode 100755
diff --git a/biocircuits/apps/promiscuous.py b/biocircuits/apps/promiscuous.py
old mode 100644
new mode 100755
diff --git a/biocircuits/dynsys.py b/biocircuits/dynsys.py
old mode 100644
new mode 100755
diff --git a/biocircuits/gillespie.py b/biocircuits/gillespie.py
old mode 100644
new mode 100755
index 252a16b..8bd43c9
--- a/biocircuits/gillespie.py
+++ b/biocircuits/gillespie.py
@@ -49,6 +49,10 @@ def _gillespie_draw(propensity_func, propensities, population, t, args):
     # Sum of propensities
     props_sum = _sum(propensities)
 
+    # Bail if the sum of propensities is zero (no moves to make)
+    if props_sum == 0.0:
+        return -1, -1.0
+
     # Compute time
     time = _draw_time(props_sum)
 
@@ -85,12 +89,16 @@ def _copy_population(population_previous, population):
             # draw the event and time step
             event, dt = draw_fun(propensity_func, propensities, population, t, args)
 
-            # Update the population
-            _copy_population(population_previous, population)
-            population += update[event, :]
+            if event == -1:
+                # Skip to the end of the simulation with the same population
+                t = time_points[-1]
+            else:
+                # Update the population
+                _copy_population(population_previous, population)
+                population += update[event, :]
 
-            # Increment time
-            t += dt
+                # Increment time
+                t += dt
 
         # Update the index (Have to be careful about types for Numba)
         j = np.searchsorted((time_points > t).astype(np.int64), 1)
@@ -105,7 +113,6 @@ def _copy_population(population_previous, population):
     return pop_out, None
 
 
-@numba.njit
 def _gillespie_trajectory_report_time_points(
     propensity_func, update, population_0, time_points, draw_fun, args=()
 ):
@@ -129,12 +136,20 @@ def _gillespie_trajectory_report_time_points(
             # draw the event and time step
             event, dt = draw_fun(propensity_func, propensities, population, t, args)
 
+            if event == -1:
+                # Skip to the end of the simulation with the same population
+                t = time_points[-1]
+            else:
+                # New population
+                population += update[event, :]
+
+                # Increment time
+                t += dt
+
             # Update the population
-            population += update[event, :]
             pop[i, :] = population
 
-            # Increment time
-            t += dt
+            # Updated time
             tp[i] = t
 
             # Increment indexes
@@ -226,6 +241,10 @@ def _draw(propensities, population, t):
             # Sum of propensities
             props_sum = np.sum(propensities)
 
+            # Bail if the sum of propensities is zero
+            if props_sum == 0.0:
+                return -1, -1.0
+
             # Compute time
             time = np.random.exponential(1 / props_sum)
 
@@ -258,12 +277,20 @@ def _traj():
                         # draw the event and time step
                         event, dt = _draw(propensities, population, t)
 
+                        if event == -1:
+                            # Skip to the end of the simulation with the same population
+                            t = time_points[-1]
+                        else:
+                            # New population
+                            population += update[event, :]
+
+                            # Increment time
+                            t += dt
+
                         # Update the population
-                        population += update[event, :]
                         pop[i, :] = population
 
-                        # Increment time
-                        t += dt
+                        # Updated time
                         tp[i] = t
 
                         # Increment indexes
@@ -299,12 +326,16 @@ def _traj():
                         # draw the event and time step
                         event, dt = _draw(propensities, population, t)
 
-                        # Update the population
-                        _copy_population(population_previous, population)
-                        population += update[event, :]
+                        if event == -1:
+                            # Skip to the end of the simulation with the same population
+                            t = time_points[-1]
+                        else:
+                            # Update the population
+                            _copy_population(population_previous, population)
+                            population += update[event, :]
 
-                        # Increment time
-                        t += dt
+                            # Increment time
+                            t += dt
 
                     # Update the index (Be careful about types for Numba)
                     j = np.searchsorted((time_points > t).astype(np.int64), 1)
@@ -321,7 +352,7 @@ def _traj():
     else:
         if return_time_points:
 
-            def traj():
+            def _traj():
                 return _gillespie_trajectory_report_time_points(
                     propensity_func,
                     update,
@@ -371,7 +402,7 @@ def _traj():
 
 
 def _gillespie_multi_fn(args):
-    """Convenient function for multithreading."""
+    """Convenience function for multithreading."""
     return _gillespie_ssa(*args)
 
 
diff --git a/biocircuits/jsfunctions.py b/biocircuits/jsfunctions.py
old mode 100644
new mode 100755
index 03ed70b..509f2b5
--- a/biocircuits/jsfunctions.py
+++ b/biocircuits/jsfunctions.py
@@ -1616,6 +1616,34 @@
 	cds.change.emit();
 }
 
+""",
+    "simple_binding_sensitivity": """
+function sensitivity(a0, b0, Kd) {
+	let b = a0 - b0 - Kd;
+	let discrim = b * b + 4 * a0 * Kd;
+	let sqrtDiscrim = Math.sqrt(discrim);
+
+	return a0 * (1 - b / sqrtDiscrim) / (sqrtDiscrim - b);
+}
+
+
+function callback() {
+    let b0 = Math.pow(10, b0Slider.value);
+    let KdVals = [0.001, 0.01, 0.1, 1.0, 10.0];        
+    let a0 = cds.data['a0'];
+    
+    for (let i = 0; i < KdVals.length; i++) {
+    	let ind = 's' + i.toString();
+    	for (let j = 0; j < a0.length; j++) {
+    		cds.data[ind][j] = sensitivity(a0[j], b0, KdVals[i]);
+    	}
+    }
+
+    span.location = b0;
+
+    cds.change.emit();
+}
+
 """,
     "toggle_nullclines": """
 function f(x, beta, n) {return beta / (1.0 + Math.pow(x, n));}
@@ -1758,4 +1786,4 @@
     cdsUnstable.change.emit();
 }
 """,
-}
\ No newline at end of file
+}
diff --git a/biocircuits/jsplots.py b/biocircuits/jsplots.py
old mode 100644
new mode 100755
index f2462aa..63ff23d
--- a/biocircuits/jsplots.py
+++ b/biocircuits/jsplots.py
@@ -718,10 +718,7 @@ def _solve_repressilator(beta, n, t_max):
     colors = colorcet.b_glasbey_category10[:3]
 
     p_rep = bokeh.plotting.figure(
-        frame_width=550,
-        frame_height=200,
-        x_axis_label="t",
-        x_range=[0, 40.0],
+        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]))
@@ -753,14 +750,11 @@ def _solve_repressilator(beta, n, t_max):
         + jsfuns["utils"]
         + jsfuns["linalg"]
         + jsfuns["proteinRepressilator"]
-        + 'callback()'
+        + "callback()"
     )
     callback = bokeh.models.CustomJS(
         args=dict(
-            cds=cds,
-            xRange=p_rep.x_range,
-            betaSlider=beta_slider,
-            nSlider=n_slider,
+            cds=cds, xRange=p_rep.x_range, betaSlider=beta_slider, nSlider=n_slider,
         ),
         code=js_code,
     )
@@ -794,7 +788,9 @@ def repressilator():
         end=4,
         step=0.1,
         value=1,
-        format=bokeh.models.FuncTickFormatter(code="return Math.pow(10, tick).toFixed(2)"),
+        format=bokeh.models.FuncTickFormatter(
+            code="return Math.pow(10, tick).toFixed(2)"
+        ),
     )
     gamma_slider = bokeh.models.Slider(
         title="γ",
@@ -802,7 +798,9 @@ def repressilator():
         end=0,
         step=0.1,
         value=0,
-        format=bokeh.models.FuncTickFormatter(code="return Math.pow(10, tick).toFixed(3)"),
+        format=bokeh.models.FuncTickFormatter(
+            code="return Math.pow(10, tick).toFixed(3)"
+        ),
     )
     rho_slider = bokeh.models.Slider(
         title="ρ",
@@ -810,11 +808,12 @@ def repressilator():
         end=0,
         step=0.1,
         value=-3,
-        format=bokeh.models.FuncTickFormatter(code="return Math.pow(10, tick).toFixed(6)"),
+        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 6-array of (dm_1/dt, dm_2/dt, dm_3/dt, dx_1/dt, dx_2/dt, dx_3/dt)
@@ -831,7 +830,6 @@ def repressilator_rhs(mx, t, beta, gamma, rho, n):
             ]
         )
 
-
     # Initial condiations
     x0 = np.array([0, 0, 0, 1, 1.1, 1.2])
 
@@ -848,13 +846,8 @@ def _solve_repressilator(log_beta, log_gamma, log_rho, n, t_max):
         m1, m2, m3, x1, x2, x3 = x.transpose()
         return t, m1, m2, m3, x1, x2, x3
 
-
     t, m1, m2, m3, x1, x2, x3 = _solve_repressilator(
-        beta_slider.value,
-        gamma_slider.value,
-        rho_slider.value,
-        n_slider.value,
-        40.0,
+        beta_slider.value, gamma_slider.value, rho_slider.value, n_slider.value, 40.0,
     )
 
     cds = bokeh.models.ColumnDataSource(
@@ -862,10 +855,7 @@ def _solve_repressilator(log_beta, log_gamma, log_rho, n, t_max):
     )
 
     p = bokeh.plotting.figure(
-        frame_width=500,
-        frame_height=200,
-        x_axis_label="t",
-        x_range=[0, 40.0],
+        frame_width=500, frame_height=200, x_axis_label="t", x_range=[0, 40.0],
     )
 
     colors = bokeh.palettes.d3["Category20"][6]
@@ -885,19 +875,13 @@ def _solve_repressilator(log_beta, log_gamma, log_rho, n, t_max):
         ("x₃", [x3_line]),
     ]
     legend = bokeh.models.Legend(items=legend_items)
-    legend.click_policy = 'hide'
+    legend.click_policy = "hide"
 
     p.add_layout(legend, "right")
 
     # Build the layout
     layout = bokeh.layouts.column(
-        bokeh.layouts.row(
-            beta_slider,
-            gamma_slider,
-            rho_slider,
-            n_slider,
-            width=575,
-        ),
+        bokeh.layouts.row(beta_slider, gamma_slider, rho_slider, n_slider, width=575,),
         bokeh.layouts.Spacer(height=10),
         p,
     )
@@ -910,7 +894,7 @@ def _solve_repressilator(log_beta, log_gamma, log_rho, n, t_max):
         + jsfuns["utils"]
         + jsfuns["linalg"]
         + jsfuns["repressilator"]
-        + 'callback()'
+        + "callback()"
     )
     callback = bokeh.models.CustomJS(
         args=dict(
@@ -933,6 +917,80 @@ def _solve_repressilator(log_beta, log_gamma, log_rho, n, t_max):
     return layout
 
 
+def simple_binding_sensitivity():
+    """Make a simple plot of sensitivity for binding of A and B.
+    """
+
+    def sensitivity(a_tot, b_tot, Kd):
+        b = a_tot - b_tot - Kd
+        discrim = b ** 2 + 4 * a_tot * Kd
+
+        return a_tot * (1 - b / np.sqrt(discrim)) / (-b + np.sqrt(discrim))
+
+    a_tot = np.logspace(-2, 2, 1001)
+    b_tot = 1.0
+    Kd_vals = np.logspace(-3, 1, 5)
+
+    colors = bokeh.palettes.Blues9[3:8]
+
+    b0_slider = bokeh.models.Slider(
+        title="total B conc. (µM)",
+        start=-1,
+        end=1,
+        value=np.log10(b_tot),
+        step=0.01,
+        width=200,
+        format=bokeh.models.FuncTickFormatter(
+            code="return Math.pow(10, tick).toFixed(2)"
+        ),
+    )
+
+    p = bokeh.plotting.figure(
+        frame_width=400,
+        height=325,
+        x_axis_label="total A conc. (µM)",
+        y_axis_label="sensitivity",
+        x_axis_type="log",
+        y_axis_type="log",
+        x_range=[a_tot.min(), a_tot.max()],
+    )
+
+    cds = bokeh.models.ColumnDataSource(
+        {
+            **{"a0": a_tot},
+            **{f"s{i}": sensitivity(a_tot, b_tot, Kd) for i, Kd in enumerate(Kd_vals)},
+        }
+    )
+
+    for i, (color, Kd) in enumerate(zip(colors, Kd_vals)):
+        p.line(
+            source=cds,
+            x="a0",
+            y=f"s{i}",
+            line_width=2,
+            color=color,
+            legend_label=str(Kd) + " µM",
+        )
+
+    span = bokeh.models.Span(location=10 ** b0_slider.value, dimension="height")
+    p.add_layout(span)
+
+    p.legend.location = "bottom_right"
+    p.legend.title = "Kd"
+
+    # Set up callbacks
+    js_code = jsfuns["simple_binding_sensitivity"] + "callback()"
+    callback = bokeh.models.CustomJS(
+        args=dict(cds=cds, b0Slider=b0_slider, span=span), code=js_code,
+    )
+
+    b0_slider.js_on_change("value", callback)
+
+    return bokeh.layouts.column(
+        bokeh.layouts.row(bokeh.models.Spacer(width=153), b0_slider), p
+    )
+
+
 def turing_dispersion_relation():
     """Plot of dispersion relation for Turing patterns.
 
diff --git a/biocircuits/rd.py b/biocircuits/rd.py
old mode 100644
new mode 100755
diff --git a/biocircuits/reg.py b/biocircuits/reg.py
old mode 100644
new mode 100755
diff --git a/biocircuits/utils.py b/biocircuits/utils.py
old mode 100644
new mode 100755
diff --git a/biocircuits/viz.py b/biocircuits/viz.py
old mode 100644
new mode 100755
index 08205e3..1a37455
--- a/biocircuits/viz.py
+++ b/biocircuits/viz.py
@@ -267,7 +267,7 @@ def interactive_xy_plot(
     """
     warnings.warn(
         DeprecationWarning,
-        "`interactive_xy_plot() is deprecated. Either custom build an interactive plot with base Bokeh or use Panel.",
+        "`interactive_xy_plot() is deprecated.",
     )
 
     def _plot_app(doc):
@@ -1451,3 +1451,55 @@ def _display_clicks(div, attributes=[], style="float:left;clear:left;font_size=0
     """
         % (attributes, style),
     )
+
+
+def phase_portrait2(
+    du_dt,
+    dv_dt,
+    u_range,
+    v_range,
+    args_u=(),
+    args_v=(),
+    log=False,
+    p=None,
+    zoomable=False,
+    **kwargs,
+):
+    """
+    Plots the phase portrait for a 2D dynamical system in the u-v plane.
+
+    Parameters
+    ----------
+    du_dt : function
+        A function to compute the right hand side of du/dt. Must have
+        call signature `du_dt(u, v, *args_u)`. Note that there can be
+        no explicit time dependence.
+    dv_dt : function
+        A function to compute the right hand side of dv/dt. Must have
+        call signature `dv_dt(u, v, *args_v)`. Note that there can be
+        no explicit time dependence.
+    u_range : array_like, shape (2,)
+        Minimum and maximum values of u to consider.
+    v_range : array_like, shape (2,)
+        Minimum and maximum values of v to consider.
+    args_u : tuple, default ()
+        Tuple of extra arguments to be passed to `du_dt`.
+    args_v : tuple, default ()
+        Tuple of extra arguments to be passed to `dv_dt`.
+    log : bool, default False
+        If True, plot u and v on a logarithmic scale.
+    p : bokeh.plotting.Figure instance, default None
+        Figure to use for the phase portrait. If None, a new one is
+        created according to `streamplot()`.
+    zoomable : bool, default False
+        If True, generates a zoomable plot where the streamlines will be
+        redrawn to correspond for the level so zoom. This requires
+        Python to be running for the plot to render.
+    kwargs :
+        All other kwargs are passed to `streamplot`.
+
+    Returns
+    -------
+    output : bokeh.plotting.Figure instance populated with streamplot
+    """
+
diff --git a/js/.gitignore b/js/.gitignore
old mode 100644
new mode 100755
diff --git a/js/circuits.js b/js/circuits.js
old mode 100644
new mode 100755
diff --git a/js/generate_jsfunctions.py b/js/generate_jsfunctions.py
old mode 100644
new mode 100755
diff --git a/js/linalg.js b/js/linalg.js
old mode 100644
new mode 100755
diff --git a/js/ode.js b/js/ode.js
old mode 100644
new mode 100755
index 1c9f2b1..9182284
--- a/js/ode.js
+++ b/js/ode.js
@@ -10,8 +10,12 @@ function rkf45(
     sBounds,
     hMin,
     enforceNonnegative,
+    event,
+    eventArgs,
     debugMode,
 ) {
+    // event must be a function with call signature event(y, t, ...eventArgs)
+
     // Set up return variables
     let tSol = [timePoints[0]];
     let t = timePoints[0];
@@ -37,11 +41,18 @@ function rkf45(
     if (hMin === undefined) hMin = 0.0;
     if (enforceNonnegative === undefined) enforceNonnegative = true;
     if (maxDeadSteps === undefined) maxDeadSteps = 10;
+    if (event === undefined) {
+        event = (function(y, t) {return true;});
+        eventArgs = [];
+    }
+    else if (eventArgs === undefined) {
+        eventArgs = [];
+    }
     if (debugMode === undefined) debugMode = false;
 
-    while (i < iMax && nDeadSteps < maxDeadSteps) {
+    while (i < iMax && nDeadSteps < maxDeadSteps && event(y, t, ...eventArgs)) {
         nDeadSteps = 0;
-        while (t < timePoints[i] && nDeadSteps < maxDeadSteps) {
+        while (t < timePoints[i] && nDeadSteps < maxDeadSteps && event(y, t, ...eventArgs)) {
             [y0, t, h, deadStep] = rkf45Step(
                 f,
                 y0,
@@ -465,7 +476,7 @@ function vsimex(
 
             // DEBUG
 		    nSteps += 1;
-		    // END EDEBUG
+		    // END DEBUG
         }
         if (t > tSol[tSol.length - 1]) {
             y.push(y0);
diff --git a/js/plot_specific/autoactivator_fixed_points.js b/js/plot_specific/autoactivator_fixed_points.js
old mode 100644
new mode 100755
diff --git a/js/plot_specific/autorepressor_response_to_pulse.js b/js/plot_specific/autorepressor_response_to_pulse.js
old mode 100644
new mode 100755
diff --git a/js/plot_specific/gaussian_pulse.js b/js/plot_specific/gaussian_pulse.js
old mode 100644
new mode 100755
diff --git a/js/plot_specific/michaelis_menten_approx.js b/js/plot_specific/michaelis_menten_approx.js
old mode 100644
new mode 100755
diff --git a/js/plot_specific/proteinRepressilator.js b/js/plot_specific/proteinRepressilator.js
old mode 100644
new mode 100755
diff --git a/js/plot_specific/repressilator.js b/js/plot_specific/repressilator.js
old mode 100644
new mode 100755
diff --git a/js/plot_specific/simple_binding_sensitivity.js b/js/plot_specific/simple_binding_sensitivity.js
new file mode 100755
index 0000000..c657f5e
--- /dev/null
+++ b/js/plot_specific/simple_binding_sensitivity.js
@@ -0,0 +1,25 @@
+function sensitivity(a0, b0, Kd) {
+	let b = a0 - b0 - Kd;
+	let discrim = b * b + 4 * a0 * Kd;
+	let sqrtDiscrim = Math.sqrt(discrim);
+
+	return a0 * (1 - b / sqrtDiscrim) / (sqrtDiscrim - b);
+}
+
+
+function callback() {
+    let b0 = Math.pow(10, b0Slider.value);
+    let KdVals = [0.001, 0.01, 0.1, 1.0, 10.0];        
+    let a0 = cds.data['a0'];
+    
+    for (let i = 0; i < KdVals.length; i++) {
+    	let ind = 's' + i.toString();
+    	for (let j = 0; j < a0.length; j++) {
+    		cds.data[ind][j] = sensitivity(a0[j], b0, KdVals[i]);
+    	}
+    }
+
+    span.location = b0;
+
+    cds.change.emit();
+}
diff --git a/js/plot_specific/toggle_nullclines.js b/js/plot_specific/toggle_nullclines.js
old mode 100644
new mode 100755
diff --git a/js/reg.js b/js/reg.js
old mode 100644
new mode 100755
diff --git a/js/rootfinding.js b/js/rootfinding.js
old mode 100644
new mode 100755
diff --git a/js/utils.js b/js/utils.js
old mode 100644
new mode 100755
diff --git a/requirements.txt b/requirements.txt
old mode 100644
new mode 100755
diff --git a/setup.cfg b/setup.cfg
old mode 100644
new mode 100755
diff --git a/setup.py b/setup.py
old mode 100644
new mode 100755
index 6e13485..9baf179
--- a/setup.py
+++ b/setup.py
@@ -2,7 +2,7 @@
 from codecs import open
 from os import path
 
-__version__ = '0.1.8'
+__version__ = '0.1.9'
 
 here = path.abspath(path.dirname(__file__))