4.4.8. Cooperative Binding Kinetics#

4.4.8.1. Motivation#

Some processes involving enzymes do not follow Michaelis-Menten kinetics. One such example is cooperative binding in which the binding of one molecule affects the binding of another molecule. In this notebook we will look at the ability of an enzyme to bind multiple copies of the same molecule.

4.4.8.2. Learning Goals#

After working through this notebook, you will be able to

  1. Write out the cooperative binding mechanism

  2. Derive the cooperative binding rate law

  3. Identify cooperative binding behavior in a plot of \(v_0\) vs \([S]_0\)

  4. Identify negative and postive cooperative binding behavior.

4.4.8.3. Coding Concepts#

The following coding concepts are used in this notebook:

  1. Variables

  2. Functions

  3. Plotting with matplotlib

4.4.8.4. Cooperative Binding#

One well-known example of a process that does not follow Michaelis-Menten kinetics is the binding of oxygen to hemoglobin. The binding of one oxygen molecule to Hemoglobin enhances the binding rate of subsequent oxygen molecules. Hemoglobins are actually heterotetramers can bind up to four oxygen molecules.

An example mechanism for this process is the simultaneous binding of \(n\) substrate molecules given by

(4.629)#\[\begin{align} E + nS &\overset{k_1}{\underset{k_{-1}}{\overset{\Longrightarrow}{\Longleftarrow}}} ES_n \\ ES_n &\overset{k_2}{\underset{k_{-2}}{\overset{\Longrightarrow}{\Longleftarrow}}} E + nP \end{align}\]

where \(n\) is the number of substrates and \(ES_n\) denotes the enzyme substrate complex with \(n\) substrates bound.

This type of behavior can be observed in a \(v_0\) vs \([S]_0\) plot by distinctly sigmoidal behavior and poor ability to fit the data with the Michaelis-Menten equation.

Hide code cell source
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
# setup plot parameters
fontsize=16
fig = plt.figure(figsize=(8,8), dpi= 80, facecolor='w', edgecolor='k')
ax = plt.subplot(111)
ax.grid(which='major', axis='both', color='#808080', linestyle='--')
ax.set_ylabel("$v_0$ ($\mu$M/s)",size=fontsize)
ax.set_xlabel("$[S]_0$ ($\mu$M)",size=fontsize)
plt.tick_params(axis='both',labelsize=fontsize)
def mm(S0,vmax,Km):  
    return vmax*S0/(Km + S0)
def hill(S0,vmax,Km,n):
    return vmax*S0**n/(Km + S0**n)
vmax = 0.001
Km = 10
S0 = np.arange(0,10,0.1)
ax.plot(S0,mm(S0,vmax,Km),'--',lw=3,label="Michaelis-Menten")
ax.plot(S0,hill(S0,vmax,Km,2),'-',lw=3,label="Non-Michaelis-Menten")
plt.legend(fontsize=fontsize);
../../_images/1f8b51ceb158fb4c7edb9aa655f278f85b4596543d42aa24c97761609c7b2820.png

4.4.8.5. Derivation of the Hill Equation#

The rate law stemming from the above mechanism can be derived in a similar fashion to the derivation for Michaelis-Menten. To make the derivation slightly easier, we will assume the second step is irreversible

(4.630)#\[\begin{align} E + nS &\overset{k_1}{\underset{k_{-1}}{\overset{\Longrightarrow}{\Longleftarrow}}} ES_n \\ ES_n &\overset{k_2}{\Longrightarrow} E + nP \end{align}\]

We start by writing the differential expressions for substrate and enzyme-substrate complex

(4.631)#\[\begin{align} -\frac{1}{n}\frac{d[S]}{dt} &= k_1[E][S]^n - k_{-1}[ES_n] \\ -\frac{d[ES_n]}{dt} &= (k_2+k_{-1})[ES_n] - k_1[E][S]^n \end{align}\]

we can again plug in \([E] = [E]_0 - [ES_n]\) to get

(4.632)#\[\begin{align} -\frac{1}{n}\frac{d[S]}{dt} &= k_1[E]_0[S]^n - (k_{-1}+k_1[S]^n)[ES_n] \\ -\frac{d[ES_n]}{dt} &= (k_1[S]^n + k_2+k_{-1})[ES_n] - k_1[E]_0[S]^n \end{align}\]

Using the steady-state approximation for \([ES_n]\) we get

(4.633)#\[\begin{equation} [ES_n] \overset{s.s.}{=} \frac{k_1[E]_0[S]^n}{k_1[S]^n + k_2+k_{-1}} \end{equation}\]

Plugging this into the rate equation

(4.634)#\[\begin{align} v_0 = -\frac{1}{n}\frac{d[S]_0}{dt} &= k_1[E]_0[S]_0^n - (k_{-1}+k_1[S]_0^n)\frac{k_1[E]_0[S]_0^n}{k_1[S]_0^n + k_2+k_{-1}}\\ &= \frac{k_1[E]_0[S]_0^n\cdot(k_1[S]_0^n + k_2+k_{-1})-(k_{-1}+k_1[S]_0^n)k_1[E]_0[S]_0^n}{k_1[S]_0^n + k_2+k_{-1}} \\ &= \frac{k_1[E]_0[S]_0^n\cdot(k_1[S]_0^n + k_2+k_{-1})-k_{-1}k_1[E]_0[S]_0^n - k_1^2[E]_0[S]_0^{2n}}{k_1[S]_0^n + k_2+k_{-1}} \\ &= \frac{k_1k_2[E]_0[S]_0^n}{k_1[S]_0^n + k_2+k_{-1}} \\ &= \frac{k_2[E]_0[S]_0^n}{[S]_0^n + K_m} \end{align}\]

where \(K_m = \frac{k_2+k_{-1}}{k_1}\) is defined similar to the Michaelis-Menten constant.

We see that the resulting initial rate law is very similar to the Michaelis-Menten rate law. Indeed, they are identical for \(n=1\) as they should be (the mechanism is then the same).

The final Hill equation is given as

(4.635)#\[\begin{equation} v_0 = \frac{v_{max}[S]_0^n}{[S]_0^n + K_m} \end{equation}\]

where \(n\) is referred to as the Hill number. While it is appealing to interpret \(n\) as the number of substrates/ligands that bind to the enzyme, the Hill mechanism is so unphysical (simultaneous binding) that this interpretation is not necessarily true. Nonetheless, this mechanism and equation are used to fit enzyme kinetics and determine whether an enzyme displays cooperativity (negative or positive) or not.

Hide code cell source
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
# setup plot parameters
fontsize=16
fig = plt.figure(figsize=(8,8), dpi= 80, facecolor='w', edgecolor='k')
ax = plt.subplot(111)
ax.grid(which='major', axis='both', color='#808080', linestyle='--')
ax.set_ylabel("$v_0$ ($\mu$M/s)",size=fontsize)
ax.set_xlabel("$[S]_0$ ($\mu$M)",size=fontsize)
plt.tick_params(axis='both',labelsize=fontsize)
def hill(S0,vmax,Km,n):
    return vmax*S0**n/(Km + S0**n)
vmax = 0.001
Km = 2
S0 = np.arange(0,10,0.1)
ax.plot(S0,hill(S0,vmax,Km,1),'-',lw=3,label="n=1")
ax.plot(S0,hill(S0,vmax,Km,2),'-',lw=3,label="n=2")
ax.plot(S0,hill(S0,vmax,Km,3),'-',lw=3,label="n=3")
ax.plot(S0,hill(S0,vmax,Km,4),'-',lw=3,label="n=4")
ax.plot(S0,hill(S0,vmax,Km,5),'-',lw=3,label="n=5")
ax.plot(S0,hill(S0,vmax,Km,6),'-',lw=3,label="n=6")
plt.legend(fontsize=fontsize);
../../_images/8a84e79cabc2f9f7ba81e2e4ca1ac41d481b5d0320a06628f1f0541f80dfcf26.png

4.4.8.6. Positive and Negative Cooperativity#

Both positive and negative cooperativity can be observed. Negative cooperativity simply means that binding of one ligand negatively impacts the binding of the second ligand. There is still cooperativity but the impact is to negate the binding of the second ligand. The Hill equation can be used to fit this type of behavior and it will be observed for values of \(n<1\). Let’s look at a plot of different \(n\) values.

Hide code cell source
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
# setup plot parameters
fontsize=16
fig = plt.figure(figsize=(8,8), dpi= 80, facecolor='w', edgecolor='k')
ax = plt.subplot(111)
ax.grid(which='major', axis='both', color='#808080', linestyle='--')
ax.set_ylabel("$v_0$ ($\mu$M/s)",size=fontsize)
ax.set_xlabel("$[S]_0$ ($\mu$M)",size=fontsize)
plt.tick_params(axis='both',labelsize=fontsize)
def hill(S0,vmax,Km,n):
    return vmax*S0**n/(Km + S0**n)
vmax = 0.001
Km = 2
S0 = np.arange(0,1000000,0.1)
ax.plot(S0,hill(S0,vmax,Km,1),'-',lw=3,label="n=1")
ax.plot(S0,hill(S0,vmax,Km,2),'-',lw=3,label="n=2")
ax.plot(S0,hill(S0,vmax,Km,3),'-',lw=3,label="n=3")
ax.plot(S0,hill(S0,vmax,Km,4),'-',lw=3,label="n=4")
ax.plot(S0,hill(S0,vmax,Km,0.1),'-',lw=3,label="n=0.1")
ax.plot(S0,hill(S0,vmax,Km,0.2),'-',lw=3,label="n=0.2")
ax.plot(S0,hill(S0,vmax,Km,0.3),'-',lw=3,label="n=0.3")
ax.plot(S0,hill(S0,vmax,Km,0.4),'-',lw=3,label="n=0.4")
ax.plot(S0,hill(S0,vmax,Km,0.5),'-',lw=3,label="n=0.5")
plt.legend(fontsize=fontsize);
/opt/anaconda3/lib/python3.13/site-packages/IPython/core/events.py:82: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  func(*args, **kwargs)
Error in callback <function flush_figures at 0x10e584220> (for post_execute), with arguments args (),kwargs {}:
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
File /opt/anaconda3/lib/python3.13/site-packages/matplotlib_inline/backend_inline.py:126, in flush_figures()
    123 if InlineBackend.instance().close_figures:
    124     # ignore the tracking, just draw and close all figures
    125     try:
--> 126         return show(True)
    127     except Exception as e:
    128         # safely show traceback if in IPython, else raise
    129         ip = get_ipython()

File /opt/anaconda3/lib/python3.13/site-packages/matplotlib_inline/backend_inline.py:90, in show(close, block)
     88 try:
     89     for figure_manager in Gcf.get_all_fig_managers():
---> 90         display(
     91             figure_manager.canvas.figure,
     92             metadata=_fetch_figure_metadata(figure_manager.canvas.figure),
     93         )
     94 finally:
     95     show._to_draw = []

File /opt/anaconda3/lib/python3.13/site-packages/IPython/core/display_functions.py:278, in display(include, exclude, metadata, transient, display_id, raw, clear, *objs, **kwargs)
    276     publish_display_data(data=obj, metadata=metadata, **kwargs)
    277 else:
--> 278     format_dict, md_dict = format(obj, include=include, exclude=exclude)
    279     if not format_dict:
    280         # nothing to display (e.g. _ipython_display_ took over)
    281         continue

File /opt/anaconda3/lib/python3.13/site-packages/IPython/core/formatters.py:238, in DisplayFormatter.format(self, obj, include, exclude)
    236 md = None
    237 try:
--> 238     data = formatter(obj)
    239 except:
    240     # FIXME: log the exception
    241     raise

File /opt/anaconda3/lib/python3.13/site-packages/decorator.py:235, in decorate.<locals>.fun(*args, **kw)
    233 if not kwsyntax:
    234     args, kw = fix(args, kw, sig)
--> 235 return caller(func, *(extras + args), **kw)

File /opt/anaconda3/lib/python3.13/site-packages/IPython/core/formatters.py:282, in catch_format_error(method, self, *args, **kwargs)
    280 """show traceback on failed format call"""
    281 try:
--> 282     r = method(self, *args, **kwargs)
    283 except NotImplementedError:
    284     # don't warn on NotImplementedErrors
    285     return self._check_return(None, args[0])

File /opt/anaconda3/lib/python3.13/site-packages/IPython/core/formatters.py:402, in BaseFormatter.__call__(self, obj)
    400     pass
    401 else:
--> 402     return printer(obj)
    403 # Finally look for special method names
    404 method = get_real_method(obj, self.print_method)

File /opt/anaconda3/lib/python3.13/site-packages/IPython/core/pylabtools.py:170, in print_figure(fig, fmt, bbox_inches, base64, **kwargs)
    167     from matplotlib.backend_bases import FigureCanvasBase
    168     FigureCanvasBase(fig)
--> 170 fig.canvas.print_figure(bytes_io, **kw)
    171 data = bytes_io.getvalue()
    172 if fmt == 'svg':

File /opt/anaconda3/lib/python3.13/site-packages/matplotlib/backend_bases.py:2157, in FigureCanvasBase.print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
   2154     # we do this instead of `self.figure.draw_without_rendering`
   2155     # so that we can inject the orientation
   2156     with getattr(renderer, "_draw_disabled", nullcontext)():
-> 2157         self.figure.draw(renderer)
   2158 if bbox_inches:
   2159     if bbox_inches == "tight":

File /opt/anaconda3/lib/python3.13/site-packages/matplotlib/artist.py:94, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
     92 @wraps(draw)
     93 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 94     result = draw(artist, renderer, *args, **kwargs)
     95     if renderer._rasterizing:
     96         renderer.stop_rasterizing()

File /opt/anaconda3/lib/python3.13/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     68     if artist.get_agg_filter() is not None:
     69         renderer.start_filter()
---> 71     return draw(artist, renderer)
     72 finally:
     73     if artist.get_agg_filter() is not None:

File /opt/anaconda3/lib/python3.13/site-packages/matplotlib/figure.py:3257, in Figure.draw(self, renderer)
   3254             # ValueError can occur when resizing a window.
   3256     self.patch.draw(renderer)
-> 3257     mimage._draw_list_compositing_images(
   3258         renderer, self, artists, self.suppressComposite)
   3260     renderer.close_group('figure')
   3261 finally:

File /opt/anaconda3/lib/python3.13/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    132 if not_composite or not has_images:
    133     for a in artists:
--> 134         a.draw(renderer)
    135 else:
    136     # Composite any adjacent images together
    137     image_group = []

File /opt/anaconda3/lib/python3.13/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     68     if artist.get_agg_filter() is not None:
     69         renderer.start_filter()
---> 71     return draw(artist, renderer)
     72 finally:
     73     if artist.get_agg_filter() is not None:

File /opt/anaconda3/lib/python3.13/site-packages/matplotlib/axes/_base.py:3226, in _AxesBase.draw(self, renderer)
   3223 if artists_rasterized:
   3224     _draw_rasterized(self.get_figure(root=True), artists_rasterized, renderer)
-> 3226 mimage._draw_list_compositing_images(
   3227     renderer, self, artists, self.get_figure(root=True).suppressComposite)
   3229 renderer.close_group('axes')
   3230 self.stale = False

File /opt/anaconda3/lib/python3.13/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    132 if not_composite or not has_images:
    133     for a in artists:
--> 134         a.draw(renderer)
    135 else:
    136     # Composite any adjacent images together
    137     image_group = []

File /opt/anaconda3/lib/python3.13/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     68     if artist.get_agg_filter() is not None:
     69         renderer.start_filter()
---> 71     return draw(artist, renderer)
     72 finally:
     73     if artist.get_agg_filter() is not None:

File /opt/anaconda3/lib/python3.13/site-packages/matplotlib/legend.py:753, in Legend.draw(self, renderer)
    749     self._legend_box.set_width(self.get_bbox_to_anchor().width - pad)
    751 # update the location and size of the legend. This needs to
    752 # be done in any case to clip the figure right.
--> 753 bbox = self._legend_box.get_window_extent(renderer)
    754 self.legendPatch.set_bounds(bbox.bounds)
    755 self.legendPatch.set_mutation_scale(fontsize)

File /opt/anaconda3/lib/python3.13/site-packages/matplotlib/offsetbox.py:369, in OffsetBox.get_window_extent(self, renderer)
    367 bbox = self.get_bbox(renderer)
    368 try:  # Some subclasses redefine get_offset to take no args.
--> 369     px, py = self.get_offset(bbox, renderer)
    370 except TypeError:
    371     px, py = self.get_offset()

File /opt/anaconda3/lib/python3.13/site-packages/matplotlib/offsetbox.py:60, in _compat_get_offset.<locals>.get_offset(self, *args, **kwargs)
     56 params = _api.select_matching_signature(sigs, self, *args, **kwargs)
     57 bbox = (params["bbox"] if "bbox" in params else
     58         Bbox.from_bounds(-params["xdescent"], -params["ydescent"],
     59                          params["width"], params["height"]))
---> 60 return meth(params["self"], bbox, params["renderer"])

File /opt/anaconda3/lib/python3.13/site-packages/matplotlib/offsetbox.py:306, in OffsetBox.get_offset(self, bbox, renderer)
    291 @_compat_get_offset
    292 def get_offset(self, bbox, renderer):
    293     """
    294     Return the offset as a tuple (x, y).
    295 
   (...)    303     renderer : `.RendererBase` subclass
    304     """
    305     return (
--> 306         self._offset(bbox.width, bbox.height, -bbox.x0, -bbox.y0, renderer)
    307         if callable(self._offset)
    308         else self._offset)

File /opt/anaconda3/lib/python3.13/site-packages/matplotlib/legend.py:722, in Legend._findoffset(self, width, height, xdescent, ydescent, renderer)
    719 """Helper function to locate the legend."""
    721 if self._loc == 0:  # "best".
--> 722     x, y = self._find_best_position(width, height, renderer)
    723 elif self._loc in Legend.codes.values():  # Fixed location.
    724     bbox = Bbox.from_bounds(0, 0, width, height)

File /opt/anaconda3/lib/python3.13/site-packages/matplotlib/legend.py:1170, in Legend._find_best_position(self, width, height, renderer)
   1163 legendBox = Bbox.from_bounds(l, b, width, height)
   1164 # XXX TODO: If markers are present, it would be good to take them
   1165 # into account when checking vertex overlaps in the next line.
   1166 badness = (sum(legendBox.count_contains(line.vertices)
   1167                for line in lines)
   1168            + legendBox.count_contains(offsets)
   1169            + legendBox.count_overlaps(bboxes)
-> 1170            + sum(line.intersects_bbox(legendBox, filled=False)
   1171                  for line in lines))
   1172 # Include the index to favor lower codes in case of a tie.
   1173 candidates.append((badness, idx, (l, b)))

File /opt/anaconda3/lib/python3.13/site-packages/matplotlib/legend.py:1170, in <genexpr>(.0)
   1163 legendBox = Bbox.from_bounds(l, b, width, height)
   1164 # XXX TODO: If markers are present, it would be good to take them
   1165 # into account when checking vertex overlaps in the next line.
   1166 badness = (sum(legendBox.count_contains(line.vertices)
   1167                for line in lines)
   1168            + legendBox.count_contains(offsets)
   1169            + legendBox.count_overlaps(bboxes)
-> 1170            + sum(line.intersects_bbox(legendBox, filled=False)
   1171                  for line in lines))
   1172 # Include the index to favor lower codes in case of a tie.
   1173 candidates.append((badness, idx, (l, b)))

File /opt/anaconda3/lib/python3.13/site-packages/matplotlib/path.py:686, in Path.intersects_bbox(self, bbox, filled)
    677 def intersects_bbox(self, bbox, filled=True):
    678     """
    679     Return whether this path intersects a given `~.transforms.Bbox`.
    680 
   (...)    684     The bounding box is always considered filled.
    685     """
--> 686     return _path.path_intersects_rectangle(
    687         self, bbox.x0, bbox.y0, bbox.x1, bbox.y1, filled)

KeyboardInterrupt: 

From the above plot we observe a number of differences between curves displaying positive cooperativity (\(n>1\)) and negative cooperativity ($n<1).

  1. Negative cooperativity displays a more rapid initial increase in \(v_0\) as compared to MM or positive cooperative behavior

  2. The rate of all value is equivalent at \([S]_0 = 1\) \(\mu\)M.

  3. The rate of negative cooperative enzymes is lower than MM or positive cooperative behvaior for \([S]_0 > 1\)

  4. Postive cooperative and MM curves all approach \(v_{max}\) asymptotically.

  5. Negative cooperative curves do not approach \(v_{max}\) asymptotically.

4.4.8.7. Fitting the Hill Equation#

When fitting the Hill equation, we must consider the Michaelis-Menten like parameters (\(v_{max}\) and \(K_m\)) as well as the Hill number, \(n\). To fit these values, we have options:

  1. Use non-linear fitting of the equation

  2. Linearize in a Lineweaver-Burk style equation

Consider the second case, we must write the reciprocal Hill equation

(4.636)#\[\begin{equation} \frac{1}{v_0} = \frac{K_m}{v_{max}}\frac{1}{[S]_0^n} + \frac{1}{v_{max}} \end{equation}\]

Observe that this equation suggests that \(\frac{1}{v_0}\) vs \(\frac{1}{[S]_0^n}\) will be linear with slope \(\frac{K_m}{v_{max}}\) and intercept \(\frac{1}{v_{max}}\). To use this equation in a linear least-squares fitting procedure, we will need to plot \(\frac{1}{v_0}\) vs \(\frac{1}{[S]_0}\), \(\frac{1}{v_0}\) vs \(\frac{1}{[S]_0^2}\), \(\frac{1}{v_0}\) vs \(\frac{1}{[S]_0^3}\), \(...\), fit lines to each one and determine the best fit by comparing \(R^2\) values. This is not ideal, and, additionally, \(n\) can take on non-integer values. We see that non-linear fitting is the appropriate way to go here.

Let’s see how it is done.

4.4.8.7.1. Example: Fitting the Hill Equation for \(n=1\)#

In this example, we will use the same data as the example from the Michaelis-Menten notes. Thus, we expect to get a fit with \(n\approx1\). The data is as follows

\([S]_0\) (mM)

\(v_0\) (\(\mu\)M/s)

1

2.5

2

4.0

5

6.3

10

7.6

20

9.0

# Put the datat into numpy arrays
s0 = np.array([1,2,5,10,20.0])
v0 = np.array([2.5,4.0,6.3,7.6,9.0])

Perform the non-linear fit and get the parameters:

Hide code cell source
# perform non-linear fit
# import least squares function from scipy library
from scipy.optimize import least_squares
# define Michaelis-Menten function
def hill(x,s):  
    return x[0]*s**x[2]/(x[1] + s**x[2])
# define residual function (difference between function and data)
def loss(x,s,data):
    return hill(x,s) - data
# make an initial guess of parameters
x0 = np.array([1.0,1.0,1.0])
res_lsq = least_squares(loss, x0, args=(s0, v0))
print("v_max = ", np.round(res_lsq.x[0],1), "micro M/s")
print("Km = ", np.round(res_lsq.x[1],1), "mM")
print("n = ", np.round(res_lsq.x[2],1), "unitless")
v_max =  10.8 micro M/s
Km =  3.2 mM
n =  0.9 unitless
Hide code cell source
# plot data
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
# setup plot parameters
fontsize=16
fig = plt.figure(figsize=(8,8), dpi= 80, facecolor='w', edgecolor='k')
ax = plt.subplot(111)
ax.grid(which='major', axis='both', color='#808080', linestyle='--')
ax.set_ylabel("$v_0$ ($\mu$M/s)",size=fontsize)
ax.set_xlabel("$[S]_0$ (mM)",size=fontsize)
plt.tick_params(axis='both',labelsize=fontsize)
ax.plot(s0,v0,'o',lw=2)
s = np.arange(np.amin(s0),np.amax(s0),0.01)
ax.plot(s,hill(res_lsq.x,s),lw=3,label="fit")
plt.legend(fontsize=fontsize);
../../_images/aa18dabf0e7798e8d4118433876d9efa8d97d15e72a39ba7e2df71d9c806cf5d.png

4.4.8.7.2. Example: Fitting the Hill Equation for \(n\) Different than 1#

Consider the following data, fit the Hill equation

\([S]_0\) (mM)

\(v_0\) (\(\mu\)M/s)

0.05

0.10

1

1.26

2

4.55

5

7.18

10

7.52

20

7.42

# Put the data into numpy arrays
import numpy as np
s0 = np.array([0.5,1,2,5,10,20.0])
v0 = np.array([0.19,1.26,4.55,7.18,7.52,7.42])
Hide code cell source
# plot data
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
# setup plot parameters
fontsize=16
fig = plt.figure(figsize=(8,8), dpi= 80, facecolor='w', edgecolor='k')
ax = plt.subplot(111)
ax.grid(which='major', axis='both', color='#808080', linestyle='--')
ax.set_ylabel("$v_0$ ($\mu$M/s)",size=fontsize)
ax.set_xlabel("$[S]_0$ (mM)",size=fontsize)
plt.tick_params(axis='both',labelsize=fontsize)
ax.plot(s0,v0,'o',lw=2);
../../_images/3ce475a3c39e8a061a19de93170a8b5998a539f669ef856fb8bd3c751c51f3d2.png
Hide code cell source
# plot data
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
# setup plot parameters
fontsize=16
fig = plt.figure(figsize=(8,8), dpi= 80, facecolor='w', edgecolor='k')
ax = plt.subplot(111)
ax.grid(which='major', axis='both', color='#808080', linestyle='--')
ax.set_ylabel("$v_0$ ($\mu$M/s)",size=fontsize)
ax.set_xlabel("$[S]_0$ (mM)",size=fontsize)
plt.tick_params(axis='both',labelsize=fontsize)
ax.plot(s0,v0,'o',lw=2)
s = np.arange(np.amin(s0),np.amax(s0),0.01)
# define Michaelis-Menten function
def hill(s,vmax,Km,n):  
    return vmax*s**n/(Km + s**n)
plt.plot(s,hill(s,7.5,5,1),lw=2,label="n=1")
plt.plot(s,hill(s,7.5,5,2),lw=2,label="n=2")
plt.plot(s,hill(s,7.5,5,3),lw=2,label="n=3")
plt.plot(s,hill(s,7.5,5,4),lw=2,label="n=4")
plt.legend(fontsize=fontsize);
../../_images/1da5b578388b6fdbe679d7944519893b41515be28fe238e5392fe093e7ca0731.png

Perform the non-linear fit to the Hill equation and get the parameters:

Hide code cell source
# perform non-linear fit
# import least squares function from scipy library
from scipy.optimize import curve_fit
# define Michaelis-Menten function
def hill(s,vmax,Km,n):  
    return vmax*s**n/(Km + s**n)
# make an initial guess of parameters
x0 = np.array([1.0,1.0])
popt, pcov = curve_fit(hill, s0, v0)
err = np.sqrt(np.diag(pcov))
print("v_max = ", np.round(popt[0],2),"+/-", np.round(err[0],2), "muM/s")
print("Km = ", np.round(popt[1],1),"+/-", np.round(err[1],1), "mM")
print("n = ", np.round(popt[2],1),"+/-", np.round(err[2],1))
v_max =  7.49 +/- 0.04 muM/s
Km =  5.0 +/- 0.2 mM
n =  2.9 +/- 0.1

Perform a non-linear fit the MM equation for reference:

Hide code cell source
# perform non-linear fit
# import least squares function from scipy library
from scipy.optimize import curve_fit
# define Michaelis-Menten function
def mm(s,vmax,Km):  
    return vmax*s/(Km + s)
# make an initial guess of parameters
x0 = np.array([1.0,1.0])
popt_mm, pcov = curve_fit(mm, s0, v0)
err_mm = np.sqrt(np.diag(pcov))
print("v_max = ", np.round(popt_mm[0],1),"+/-", np.round(err_mm[0],1), "muM/s")
print("Km = ", np.round(popt_mm[1],1),"+/-", np.round(err_mm[1],1), "mM")
v_max =  9.4 +/- 1.5 muM/s
Km =  2.8 +/- 1.4 mM

Plot the results:

Hide code cell source
# plot data
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
# setup plot parameters
fontsize=16
fig = plt.figure(figsize=(8,8), dpi= 80, facecolor='w', edgecolor='k')
ax = plt.subplot(111)
ax.grid(which='major', axis='both', color='#808080', linestyle='--')
ax.set_ylabel("$v_0$ ($\mu$M/s)",size=fontsize)
ax.set_xlabel("$[S]_0$ (mM)",size=fontsize)
plt.tick_params(axis='both',labelsize=fontsize)
ax.plot(s0,v0,'o',lw=2)
s = np.arange(0.5,20,0.01)
ax.plot(s,hill(s,popt[0],popt[1],popt[2]),lw=2,label="Hill")
ax.plot(s,mm(s,popt_mm[0],popt_mm[1]),lw=2,label="Michaelis-Menten")
plt.legend(fontsize=fontsize);
<matplotlib.legend.Legend at 0x7fd400bc9eb0>
../../_images/b67b6fdac8608549dddf9f5e40c79536bc39cd5e20203b39d018c9c18a1df341.png

4.4.8.8. Example: Multiple Data Sets#

Hide code cell source
from tabulate import tabulate
import numpy as np
s0 = np.array([1,2,5,10,20.0])
# Generate a data set
def hill(S0,vmax,Km,n):  
    return vmax*S0**n/(Km + S0**n)
vmax = 7.5
Km = 4.0
n = 0.5
truth = hill(s0,vmax,Km,n)
n_trials = 1
data = np.empty((truth.shape[0],n_trials))
s0_total = np.empty((s0.shape[0],n_trials))
for i in range(n_trials):
    # estimate error based on normal distribution 99.9% data within 7.5%
    error = np.random.normal(0,0.03,truth.shape[0])
    # estimate error from uniform distribution with maximum value of 5%
    #error = 0.1*(np.random.rand(truth.shape[0])-0.5)
    # generate data by adding error to truth
    data[:,i] = truth*(1+error)
    # keep flattened s0 array
    s0_total[:,i] = s0
combined_data = np.column_stack((s0,data))
print(tabulate(combined_data,headers=["[S]0","Trial 1", "Trial 2", "Trial 3", "Trial 4", "Trial 5"]))
  [S]0    Trial 1
------  ---------
     1    1.5454
     2    2.0766
     5    2.75409
    10    3.23132
    20    3.74453
# perform non-linear fit
# import least squares function from scipy library
from scipy.optimize import curve_fit
# define Michaelis-Menten function
def hill(s,vmax,Km,n):  
    return vmax*s**n/(Km + s**n)
# make an initial guess of parameters
popt, pcov = curve_fit(hill, s0_total.flatten(), data.flatten())
err = np.sqrt(np.diag(pcov))
print("v_max = ", np.round(popt[0],1),"+/-", np.round(err[0],1), "muM/s")
print("Km = ", np.round(popt[1],1),"+/-", np.round(err[1],1), "mM")
print("n = ", np.round(popt[2],1),"+/-", np.round(err[2],1))
v_max =  5.2 +/- 0.4 muM/s
Km =  2.3 +/- 0.2 mM
n =  0.6 +/- 0.1