## ~dhogan/rnap

e889dcf373785aaecda1ebff2aaa992cce41c3cc — Daniel Hogan 2 years ago
```add michaelis-menten kinetics
```
```1 files changed, 86 insertions(+), 0 deletions(-)

A mm.py
```
`A mm.py => mm.py +86 -0`
```@@ 0,0 1,86 @@
+#!/usr/bin/env python3
+
+import numpy as np
+from scipy.optimize import curve_fit as cf
+from scipy.stats import norm, chi2
+import matplotlib.pyplot as plt
+import warnings
+
+warnings.filterwarnings(
+    "ignore", message="Covariance of the parameters could not be estimated"
+)
+
+
+def summarize(data):
+    estimate = np.median(data)
+    edge = (100 - 68) / 2
+    return (
+        estimate,
+        np.percentile(data, edge) - estimate,
+        np.percentile(data, 100 - edge) - estimate,
+    )
+
+
+def mm(conc, vmax, km):
+    return vmax * conc / (km + conc)
+
+
+data = {}
+data[("hairpin", 10)] = ((250, 1000), (7.27, 3.07), (1.02, 0.16), (0.80, 0.14))
+data[("hairpin", 5)] = ((250, 1000), (10.93, 3.71), (1.76, 0.80), (1.33, 0.56))
+data[("no hairpin", 10)] = (
+    (25, 50, 250),
+    (6.91, 5.39, 3.17),
+    (0.89, 1.27, 0.16),
+    (0.71, 0.86, 0.14),
+)
+data[("termination", 10)] = (
+    (25, 50, 250, 1000),
+    (4.21, 4.41, 2.63, 1.51),
+    (1.01, 1.13, 0.48, 0.14),
+    (0.68, 1.13, 0.35, 0.12),
+)
+data[("termination", 5)] = ((250, 1000), (2.71, 1.74), (0.62, 0.62), (0.43, 0.36))
+
+pauses = ("hairpin", "no hairpin", "termination")
+forces = (5, 10)
+
+plot = False
+n = 10000
+
+output = []
+
+for pause in pauses:
+    for force in forces:
+        try:
+            x, y, plus, minus = data[(pause, force)]
+        except:
+            continue
+        print(pause, force)
+        y = np.array(y)
+        yerr = 1 / y - 1 / (y + plus)
+        tmins = []
+        kms = []
+        invyrands = []
+        for _ in range(n):
+            invyrand = 1 / y + yerr * norm.rvs()
+            invyrands.append(invyrand)
+            popt, pcov = cf(mm, x, invyrand)
+            tmins.append(1/popt[0])
+            kms.append(popt[1])
+        xs = 10 * np.arange(125)
+        print(summarize(tmins))
+        print(summarize(kms))
+        output.append(pause, force, summarize(tmins), summarize(kms))
+        print()
+        if not plot:
+            continue
+        plt.hist(vmaxes, bins="auto")
+        plt.show()
+        plt.hist(kms, bins="auto")
+        plt.show()
+        for vmax, km, invyrand in zip(vmaxes, kms, invyrands):
+            plt.plot(xs, mm(xs, vmax, km), alpha=1 / 100)
+            plt.scatter(x, invyrand, alpha=1 / 100)
+        plt.errorbar(x, 1 / y, yerr=yerr)
+        plt.show()

```