Merge pull request #1 from Becksteinlab/add_leak

Added leak models and mean first arrival time simulator

.gitignore

@@ -1,2 +1,3 @@
 *__pycache__*
-*.ipynb_checkpoint*
+*.ipynb_checkpoint*
+effective_rates/target/

full_antiporter_example/equil.py

@@ -100,7 +100,8 @@ def run(scanner : MultibindScanner, basepath : Union[str, Path, None] = None) ->
              ('IFNA', 'OFNA'),
              ('OFNA', 'OF0'),
              ('OF0', 'OFH'),
-             ('OFH', 'IFH')]
+             ('OFH', 'IFH'),
+             ('OF0', 'IF0')]
 
     concentrations['H+'] = [8]
     scanner.run(concentrations, svd=False)
@@ -193,4 +194,4 @@ def plot_msp(pH, Na, scanner, outdir):
 
 
 if __name__ == "__main__":
-    run()
+    pass

full_antiporter_example/generate_paper_tables.py

@@ -1,3 +1,5 @@
+#!/usr/bin/env python
+
 from multibind.nonequilibrium import rate_matrix
 import math
 import pathlib
@@ -18,16 +20,16 @@
             (ofn, ifn),
             (ifn, if0),
             (if0, ifh),
-            ]
+            (of0, if0)]
 
 
-def dG2pKa(dG : float, pH : float = 0.0) -> float:
+def dG2pKa(dG: float, pH: float = 0.0) -> float:
     '''Convert Delta G to pKa given the pH.
     '''
     return pH - dG / math.log(10)
 
 
-def format_name(name : str):
+def format_name(name: str):
     '''Take in a state name and format it for latex.
     '''
     if "NA" in name:
@@ -55,7 +57,10 @@ def table_from_entries(entries, bars=True, dG_err=None) -> None:
                 entry = e
                 break
         if not entry:
+            if s1 == of0 and s2 == if0:
+                continue
             raise ValueError
+
         ordered_entries.append(entry)
 
     if dG_err:
@@ -68,7 +73,7 @@ def table_from_entries(entries, bars=True, dG_err=None) -> None:
         dG = f"{-math.log(k / bk):0.3f}"
         dG_std = math.sqrt(var / k**2 + bvar / bk**2)
 
-        if dG_err:
+        if dG_err:  # only used with the corrected rates
             dG_std = dG_err[i]
 
         append_value = f"{format_name(s1)} & {format_name(s2)} & ${k:0.2f} \\pm {var**0.5:0.2f}$ & ${bk:0.2f} \\pm {bvar**0.5:0.2f}$  & ${dG} \\pm {dG_std:0.3f}$ \\\\"
@@ -86,7 +91,7 @@ def table_from_entries(entries, bars=True, dG_err=None) -> None:
     print(table)
 
 
-def raw_rates_table(rate_file : Union[str, pathlib.Path]) -> None:
+def raw_rates_table(rate_file: Union[str, pathlib.Path]) -> None:
     '''Generate latex table from the raw rates file and print to screen.
     '''
     entries = []
@@ -97,11 +102,10 @@ def raw_rates_table(rate_file : Union[str, pathlib.Path]) -> None:
                 continue
             s1, s2, v, sigma = line
             entries.append((s1, s2, float(v), float(sigma)))
-
     table_from_entries(entries, bars=True)
 
 
-def corrected_rates_table(rate_file : Union[str, pathlib.Path]) -> None:
+def corrected_rates_table(rate_file: Union[str, pathlib.Path]):
     '''Generate latex table from the multibind corrected rates and print to screen.
     '''
     pH = 8
@@ -117,7 +121,6 @@ def corrected_rates_table(rate_file : Union[str, pathlib.Path]) -> None:
 
     for index, data in new_graph.iterrows():
         state1, state2, value, variance, ligand, std = data
-        # print(state1, state2, value, variance, ligand, std)
 
         if (state1[-1] == "H" and state2[-1] == "0") or (state1[-1] == "A" and state2[-1] == "0"):
             # backwards proton reaction
@@ -138,7 +141,7 @@ def corrected_rates_table(rate_file : Union[str, pathlib.Path]) -> None:
             new_graph.at[index, 'ligand'] = "Na+"
             new_graph.at[index, 'value'] = new_graph.value[index] + math.log(Na)
 
-    dG_err = []
+    dG_err = []  # collect errors from Cramer rao
     entries = []
     for s1, s2 in ordering:
         s1_idx = list(filter(lambda x: x[1][0] == s1, enumerate(states)))[0][0]
@@ -150,6 +153,35 @@ def corrected_rates_table(rate_file : Union[str, pathlib.Path]) -> None:
 
     table_from_entries(entries, bars=False, dG_err=dG_err)
 
+    return entries
+
+
+def corrected_rates_table_to_tuple(entries):
+
+    rows = ""
+
+    def corrected_ligand_name(name):
+        if name == "H":
+            return "H"
+        elif name == "NA":
+            return "Na"
+        else:
+            return "Empty"
+
+    for s1, s2, rate, sigma in entries:
+        conf1 = s1[0:2]
+        ligand1 = corrected_ligand_name(s1[2:])
+        conf2 = s2[0:2]
+        ligand2 = corrected_ligand_name(s2[2:])
+
+        rows += f"    (State({conf1}, {ligand1}), State({conf2}, {ligand2}), {rate}, {sigma**0.5}),\n"
+
+    command = f"""
+rates = (
+{rows})
+"""
+    print(command)
+
 
 def main():
     from sys import argv
@@ -164,8 +196,13 @@ def main():
     print('=======================================\n\n')
 
     print('=========== CORRECTED RATES ===========\n')
-    corrected_rates_table(rate_file)
+    entries = corrected_rates_table(rate_file)
+    print('=======================================\n\n')
+
+    print('========== CORRECTED RATES* ===========\n')
+    corrected_rates_table_to_tuple(entries)
     print('=======================================\n\n')
 
+
 if __name__ == "__main__":
     main()

full_antiporter_example/inputs/diffusion_rates.csv

@@ -10,4 +10,6 @@ OF0,OFNA,320000000,2890000000000
 OF0,OFH,200,225
 OFH,OF0,6000,10000
 OFH,IFH,5000,10000
-IFH,OFH,8000,10000
+IFH,OFH,8000,10000
+IF0,OF0,135,1000
+OF0,IF0,100,1000

full_antiporter_example/inputs/diffusion_rates_large.csv

@@ -0,0 +1,15 @@
+state1,state2,k,var
+IFH,IF0,8000,160000
+IF0,IFH,200,625
+IF0,IFNA,320000000,2890000000000
+IFNA,IF0,70000000,90250000000000
+IFNA,OFNA,5000,10000
+OFNA,IFNA,8000,10000
+OFNA,OF0,170000000,100000000000000
+OF0,OFNA,320000000,2890000000000
+OF0,OFH,200,225
+OFH,OF0,6000,10000
+OFH,IFH,5000,10000
+IFH,OFH,8000,10000
+IF0,OF0,1350,1000
+OF0,IF0,1000,1000

full_antiporter_example/inputs/diffusion_rates_none.csv

@@ -0,0 +1,13 @@
+state1,state2,k,var
+IFH,IF0,8000,160000
+IF0,IFH,200,625
+IF0,IFNA,320000000,2890000000000
+IFNA,IF0,70000000,90250000000000
+IFNA,OFNA,5000,10000
+OFNA,IFNA,8000,10000
+OFNA,OF0,170000000,100000000000000
+OF0,OFNA,320000000,2890000000000
+OF0,OFH,200,225
+OFH,OF0,6000,10000
+OFH,IFH,5000,10000
+IFH,OFH,8000,10000

full_antiporter_example/pka_kd_values.py

@@ -60,24 +60,58 @@ def __str__(self) -> str:
 Na = "Na+"
 Empty = "Empty"
 
-rates = (
-    # OF
-    (State(OF, Empty), State(OF, H),           194.45,    3.58),
-    (State(OF, H),     State(OF, Empty),      6007.89,    4.14),
-    (State(OF, Empty), State(OF, Na),    320043930.33,  395.36),
-    (State(OF, Na),    State(OF, Empty), 167088412.30, 3259.91),
-    # IF
-    (State(IF, Empty), State(IF, H),           215.89,    4.12),
-    (State(IF, H),     State(IF, Empty),      7888.67,   31.22),
-    (State(IF, Empty), State(IF, Na),    319948241.25, 6504.61),
-    (State(IF, Na),    State(IF, Empty),  76739017.89, 2677.41),
-    # CONF
-    (State(IF, Na),    State(OF, Na),         4990.12,   18.53),
-    (State(IF, Na),    State(OF, Na),         8006.22,   12.17),
-    (State(IF, H),     State(OF, H),          8006.22,   12.17),
-    (State(IF, H),     State(OF, H),          4990.12,   18.53),
+none_rates = (
+    (State(IF, H), State(OF, H), 8006.223169408338, 12.57314756359655),
+    (State(OF, H), State(IF, H), 4990.015244309923, 19.140787642841556),
+    (State(OF, H), State(OF, Empty), 6007.98133761928, 4.139689625469715),
+    (State(OF, Empty), State(OF, H), 194.45150148101385, 3.5841640766354605),
+    (State(OF, Empty), State(OF, Na), 320043930.3290041, 395.36290105497517),
+    (State(OF, Na), State(OF, Empty), 167088412.30087355, 3259.9134983152007),
+    (State(OF, Na), State(IF, Na), 8006.223169408335, 12.169381784880281),
+    (State(IF, Na), State(OF, Na), 4990.01524430993, 18.526112996833973),
+    (State(IF, Na), State(IF, Empty), 76739017.89163141, 2677.4053171912883),
+    (State(IF, Empty), State(IF, Na), 319948241.25336146, 6504.605667410289),
+    (State(IF, Empty), State(IF, H), 215.89117545131967, 4.123301815264508),
+    (State(IF, H), State(IF, Empty), 7888.666130672948, 31.218756052842668),
+    (State(OF, Empty), State(IF, Empty), 0.0, 0.0),
+    (State(IF, Empty), State(OF, Empty), 0.0, 0.0),
 )
 
+small_rates = (
+    (State(IF, H), State(OF, H), 8006.205593777888, 12.569971591604132),
+    (State(OF, H), State(IF, H), 4990.043521784862, 19.13587689337721),
+    (State(OF, H), State(OF, Empty), 6007.959673036285, 4.120731730091776),
+    (State(OF, Empty), State(OF, H), 194.46702398836024, 3.567464723819803),
+    (State(OF, Empty), State(OF, Na), 320044062.77480483, 392.94891103734113),
+    (State(OF, Na), State(OF, Empty), 167079476.74613082, 3240.1849035495575),
+    (State(OF, Na), State(IF, Na), 8006.242380110126, 12.146920571454022),
+    (State(IF, Na), State(OF, Na), 4989.984335655457, 18.491999128412196),
+    (State(IF, Na), State(IF, Empty), 76760846.63180155, 2645.7557493785434),
+    (State(IF, Empty), State(IF, Na), 319948058.79832435, 6426.80159525659),
+    (State(IF, Empty), State(IF, H), 215.84542374567658, 4.076852122475119),
+    (State(IF, H), State(IF, Empty), 7889.0151036064735, 30.86985335989601),
+    (State(OF, Empty), State(IF, Empty), 99.70388815206579, 2.401037558616935),
+    (State(IF, Empty), State(OF, Empty), 135.2183394847137, 3.789998694486533),
+)
+
+large_rates = (
+    (State(IF, H), State(OF, H), 8006.002877514489, 12.533154300899445),
+    (State(OF, H), State(IF, H), 4990.369640303549, 19.07895672839883),
+    (State(OF, H), State(OF, Empty), 6007.709441691196, 3.8787346954443973),
+    (State(OF, Empty), State(OF, H), 194.64612022231438, 3.354863361763936),
+    (State(OF, Empty), State(OF, Na), 320045588.7058398, 361.3535690488545),
+    (State(OF, Na), State(OF, Empty), 166976457.88334376, 2981.519619855568),
+    (State(OF, Na), State(IF, Na), 8006.463912125178, 11.878496964270255),
+    (State(IF, Na), State(OF, Na), 4989.6278679430425, 18.08426456609654),
+    (State(IF, Na), State(IF, Empty), 77013041.25217529, 2147.63862169223),
+    (State(IF, Empty), State(IF, Na), 319945943.9022867, 5208.28444500144),
+    (State(IF, Empty), State(IF, H), 215.31805431858123, 3.3692238082462045),
+    (State(IF, H), State(IF, Empty), 7893.0254385417575, 25.53827745474294),
+    (State(OF, Empty), State(IF, Empty), 999.6291938053557, 5.654090509498371),
+    (State(IF, Empty), State(OF, Empty), 1350.2745135849646, 8.920341089700932),
+)
+
+all_rates = {"small": small_rates, "none": none_rates, "large": large_rates}
 
 c_H = 1e-8    # pH = 8
 c_Na = 0.100  # [Na^+] = 100 mM
@@ -125,6 +159,15 @@ def main():
 
     The forward and reverse rates are hardcoded globally.
     """
+    import sys
+
+    try:
+        target = sys.argv[1]
+    except:
+        target = "small"
+
+    rates = all_rates[target]
+
     IFNA = State(IF, Na)
     IFH = State(IF, H)
     IF0 = State(IF, Empty)

full_antiporter_example/transport.py

@@ -119,6 +119,8 @@ def transport(c, rates, h_counter, na_counter, **kwargs):
 
     states = c.states.values
 
+    target_flux = None
+
     for i, j in connections:
         i_i = np.argwhere(states == i)[0, 0]
         i_j = np.argwhere(states == j)[0, 0]
@@ -132,4 +134,10 @@ def transport(c, rates, h_counter, na_counter, **kwargs):
         drive.append(-np.log(Gp.T[i_i, i_j] / Gp.T[i_j, i_i]))
         connection_labels.append(f"{i}/{j}")
 
-    return net[0], drive, connection_labels, steady_state_populations, Gp.T
+        if (j == "IFH") and (i == "OFH"):
+            target_flux = net[-1]
+
+    if not target_flux:
+        raise RuntimeError("Didn't find operational edge")
+
+    return -target_flux, drive, connection_labels, steady_state_populations, Gp.T