initial commit of the kinetic module

README.md

-# kinetic
+# Kinetic
 
+`kinetic` is a Python module collecting functions useful for simulation of chemical and enzymatic kinetics.
+
+# Functions
+The main functionality of the module is covered by:
+Function `reaction(rxn_str)`, which takes one argument:
+* `rxn_str` text represenation of the reaction. The reaction specification syntax is simple, with two reaction types available: **reversible** by using `<=>` symbol  or **irreversible** `->`. The reagents are separated by `+` sign and spaces. The rate constants are named at the end after the `;` symbol.
+
+Function `concentration(rxns, c0, rates, ts)` takes four arguments:
+* `rxns` *(type: list)* a list of reactions, created by the **reaction** function
+* `c0` *(type: dict)* a dictionary of initial conditions. Only species with non-zero concentrations should be included
+* `rates` *(type: dict)* a dictionary holding the values of all defined kinetic rate constants
+* `ts` *(type: np.array)* an array with discrete time points, for which the transient concentration will be calculated
+
+Simple example:
+```
+rxns = [kn.reaction('E + S <=> ES ; k1on k1off'),
+        kn.reaction('ES -> E + P  ; k2')]
+
+c0 = {'E': 1.0e-3, 'S': 1.0}
+kin = {'k1on': 1.0e3, 'k1off': 200.0, 'k2': 250.0}
+ts = np.linspace(0.0, 20, 21)
+
+spcs, ct = kn.concentration(rxns, c0, kin, ts)
+```

example1.py

+import matplotlib.pyplot as plt
+import numpy as np
+
+import kinetic as kn
+
+# E + S <=> ES
+# ES -> E + P
+
+rxns = [kn.reaction('E + S <=> ES ; k1on k1off'),
+        kn.reaction('ES -> E + P  ; k2')]
+
+c0 = {'E': 1.0e-3, 'S': 1.0}
+kin = {'k1on': 1.0e3, 'k1off': 200.0, 'k2': 250.0}
+ts = np.linspace(0.0, 20, 21)
+
+spcs, ct = kn.concentration(rxns, c0, kin, ts)
+
+print(spcs)
+for spc in spcs:
+    plt.plot(ts, ct[:, spcs[spc]], '.-', label=spc)
+
+plt.grid()
+plt.legend()
+plt.show()
+
+for i in range(len(ts)):
+    print('%8.3f %12.6e' % (ts[i], ct[i, spcs['S']]))

kinetic.py

+'''
+Kinetic reaction modelling in Python
+'''
+
+import unittest
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.integrate import odeint
+
+
+def parse_reagents(reagents_desc):
+    '''
+    creates reagents dictionary from the string
+
+    Examples:
+      A + B    produces {'A':1, 'B':1}
+    '''
+    reagents = {}
+
+    for rgnt in reagents_desc:
+        rgnt_label = ''
+        rgnt_coeff = 0
+        if rgnt[0].isdigit():
+            rgnt_label = rgnt[1:]
+            rgnt_coeff = int(rgnt[0])
+        else:
+            rgnt_label = rgnt
+            rgnt_coeff = 1
+        # assign substrate or increase coeff
+        if rgnt_label in reagents.keys():
+            reagents[rgnt_label] += rgnt_coeff
+        else:
+            reagents[rgnt_label] = rgnt_coeff
+
+    return reagents
+
+
+def reaction(rxn_str):
+    '''
+    creates reaction dictionary from a string
+
+    Examples:
+      A -> B    ; k1
+      2A <=> A2   ; k2on k2off
+    '''
+    reaction_type = ''
+    rateconst = {'forward': None, 'backward': None}
+
+    rxn = rxn_str.split()
+
+    # remove pluses
+    while '+' in rxn:
+        rxn.remove('+')
+
+    # determine reaction type
+    if rxn.count('->') == 1:
+        reaction_type = 'irreversible'
+        rateconst['forward'] = rxn.pop()
+        rxn.pop()  # pop ';' character
+        index = rxn.index('->')
+    elif rxn.count('<=>') == 1:
+        reaction_type = 'reversible'
+        rateconst['backward'] = rxn.pop()
+        rateconst['forward'] = rxn.pop()
+        rxn.pop()  # pop ';' character
+        index = rxn.index('<=>')
+    else:
+        raise Exception('sth\ss wrong here %s' % rxn_str)
+
+    left_side = rxn[:index]
+    right_side = rxn[index + 1:]
+
+    # parse substrates
+    substrates = parse_reagents(left_side)
+
+    # parse products
+    products = parse_reagents(right_side)
+
+    return {'type': reaction_type,
+            'substrates': substrates,
+            'products': products,
+            'rateconst': rateconst}
+
+
+def get_species(rxns):
+    species = {}
+    i = 0
+    for rxn in rxns:
+        for substrate in rxn['substrates'].keys():
+            if not substrate in species.keys():
+                species[substrate] = i
+                i += 1
+        for product in rxn['products'].keys():
+            if not product in species.keys():
+                species[product] = i
+                i += 1
+    return species
+
+
+def make_conc_dict(species, ct):
+    cnc_dict = {}
+    for spc, spc_id in species.items():
+        cnc_dict[spc] = ct[spc_id]
+    return cnc_dict
+
+
+def print_kinetic_equations(rxns):
+    species = get_species(rxns)
+    print(species)
+
+    rates = [['', ''] for i in range(len(species))]
+    for spc_label, spc_id in species.items():
+        rates[spc_id][0] = 'd[%s]/dt' % spc_label
+
+    for rxn in rxns:
+        on_rate = rxn['rateconst']['forward']
+        for substrate_label, substrate_coeff in rxn['substrates'].items():
+            on_rate += '[' + substrate_label + ']'
+            if substrate_coeff != 1:
+                on_rate += '^%d' % substrate_coeff
+
+        off_rate = ''
+        if rxn['type'] == 'reversible':
+            off_rate = rxn['rateconst']['backward']
+            for product_label, product_coeff in rxn['products'].items():
+                off_rate += '[' + product_label + ']'
+                if product_coeff != 1:
+                    off_rate += '^%d' % product_coeff
+
+        # Substrates
+        for substrate_label, substrate_coeff in rxn['substrates'].items():
+            if substrate_coeff != 1:
+                rates[species[substrate_label]
+                      ][1] += '  -%d*' % substrate_coeff + on_rate
+            else:
+                rates[species[substrate_label]][1] += '  -' + on_rate
+
+            if rxn['type'] == 'reversible':
+                if substrate_coeff != 1:
+                    rates[species[substrate_label]
+                          ][1] += '  +%d*' % substrate_coeff + off_rate
+                else:
+                    rates[species[substrate_label]][1] += '  +' + off_rate
+
+        # Products
+        for product_label, product_coeff in rxn['products'].items():
+            if product_coeff != 1:
+                rates[species[product_label]
+                      ][1] += '  +%d*' % product_coeff + on_rate
+            else:
+                rates[species[product_label]][1] += '  +' + on_rate
+
+            if rxn['type'] == 'reversible':
+                #rates[species[product_label]][1] += '  -' + off_rate
+                if product_coeff != 1:
+                    rates[species[product_label]
+                          ][1] += '  -%d*' % product_coeff + off_rate
+                else:
+                    rates[species[product_label]][1] += '  -' + off_rate
+
+    # print rates
+    print
+    for i in range(len(rates)):
+        print('%15s =   %-s' % tuple(rates[i]))
+
+    return
+
+
+def concentration(rxns, initial, kinCoeff, timepoints, what='all'):
+    '''
+    solves diff equations and returns the transient concentrations
+    '''
+    # get all the species
+    species = {}
+    i = 0
+    for rxn in rxns:
+        for substrate in rxn['substrates'].keys():
+            if not substrate in species.keys():
+                species[substrate] = i
+                i += 1
+        for product in rxn['products'].keys():
+            if not product in species.keys():
+                species[product] = i
+                i += 1
+
+    # print species
+
+    # set initial conct
+    c0 = {}
+    for spc in species:
+        c0[spc] = 0.0
+    for label in initial:
+        c0[label] = initial[label]
+
+    # print c0
+
+    # set y0
+    y0 = np.zeros(len(species))
+    for spc in species:
+        y0[species[spc]] = c0[spc]
+
+    # print 'y0:', y0
+
+    # construct the dy function
+
+    def dy(t, y):
+        rates = np.zeros(len(species))
+
+        for rxn in rxns:
+            # print rxn
+            # on_rate = str(rxn['rateconst']['forward']) + '*' + str(rxn['substrates'].keys())
+            # off_rate = ''
+            on_rate = kinCoeff[rxn['rateconst']['forward']]
+            for substrate_label, substrate_coeff in rxn['substrates'].items():
+                on_rate *= y[species[substrate_label]] ** substrate_coeff
+
+            off_rate = 0.0
+            if rxn['type'] == 'reversible':
+                off_rate = kinCoeff[rxn['rateconst']['backward']]
+                # off_rate = str(rxn['rateconst']['backward']) + '*' + str(rxn['products'].keys())
+                for product_label, product_coeff in rxn['products'].items():
+                    off_rate *= y[species[product_label]] ** product_coeff
+
+            # Substrates
+            # for substrate in rxn['substrates']:
+            for substrate_label, substrate_coeff in rxn['substrates'].items():
+                # print substrate,
+                # rates[species[substrate]] = rates[species[substrate]] + '  -' + on_rate
+                rates[species[substrate_label]] -= substrate_coeff * on_rate
+
+                if rxn['type'] == 'reversible':
+                    # rates[species[substrate]] = rates[species[substrate]] + '  +' + off_rate
+                    rates[species[substrate_label]
+                          ] += substrate_coeff * off_rate
+
+            # Products
+            # for product in rxn['products']:
+            for product_label, product_coeff in rxn['products'].items():
+                # print product,
+                # rates[species[product]] = rates[species[product]] + '  +' + on_rate
+                rates[species[product_label]] += product_coeff * on_rate
+                if rxn['type'] == 'reversible':
+                    # rates[species[product]] = rates[species[product]] + '  -' + off_rate
+                    rates[species[product_label]] -= product_coeff * off_rate
+        # print rates
+        return rates
+
+    # TODO: test it separately
+    # dy(y0, [1.0])
+    # return
+
+    ct = odeint(dy, y0, timepoints, tfirst=True)
+
+    return (species, ct)
+
+
+# ## Tests
+
+class reaction_test(unittest.TestCase):
+    def test_simple_irreversible(self):
+        r1 = reaction('A -> B ; kon')
+        # print r1
+        assert r1['type'] == 'irreversible'
+        assert 'A' in r1['substrates'].keys()
+        assert not 'B' in r1['substrates'].keys()
+        assert 'B' in r1['products'].keys()
+        assert not 'A' in r1['products'].keys()
+        assert r1['rateconst']['forward'] == 'kon'
+        assert r1['rateconst']['backward'] == None
+
+    def test_simple_reversible(self):
+        r2 = reaction('2A + B <=> A2B ; k2f k2b')
+        # print r2
+        assert r2['type'] == 'reversible'
+        assert 'A' in r2['substrates'].keys()
+        assert 'B' in r2['substrates'].keys()
+        assert not '+' in r2['substrates'].keys()
+        assert 'A2B' in r2['products'].keys()
+        assert not 'A' in r2['products'].keys()
+        assert r2['rateconst']['forward'] == 'k2f'
+        assert r2['rateconst']['backward'] == 'k2b'
+
+    def test_simple_dimerization(self):
+        r3 = reaction('A + 3A <=> B + B + 6B ; k3f k3b')
+        print(r3)
+        print_kinetic_equations([r3])
+        assert r3['type'] == 'reversible'
+        assert 'A' in r3['substrates'].keys()
+        assert not 'B' in r3['substrates'].keys()
+        assert 'B' in r3['products'].keys()
+        assert not 'A' in r3['products'].keys()
+        assert r3['rateconst']['forward'] == 'k3f'
+        assert r3['rateconst']['backward'] == 'k3b'
+
+    def test_enzymatic_concentration(self):
+        '''
+        test E + S <=> ES -> E + P reaction scheme
+        '''
+        rxns = [reaction('E + S <=> ES ; k1on k1off'),
+                reaction('ES -> E + P ; k2')]
+        c0 = {'E': 1.0e-3, 'S': 0.1}
+        kin = {'k1on': 1.0e6, 'k1off': 100, 'k2': 5.0}
+        ts = np.linspace(0.0, 30, 100)
+
+        spcs, ct = concentration(rxns, c0, kin, ts)
+        print(spcs)
+        for spc in spcs:
+            plt.plot(ts, ct[:, spcs[spc]], '.-', label=spc)
+
+        # plt.plot(ts, ct[:, 0], 'k.-', label='S')
+        # plt.plot(ts, ct[:, 2], 'b.-', label='ES')
+        # plt.plot(ts, ct[:, 3], 'r.-', label='P')
+        plt.grid()
+        plt.legend()
+        plt.show()
+
+        for i in range(len(ts)):
+            print('%8.3f %12.6e' % (ts[i], ct[i, spcs['S']]))
+        # test some concentrations
+
+        assert ct[-1, spcs['S']] <= 1.0e-9
+
+    def test_dimerization_concentration(self):
+        '''
+        test A + A <=> B reaction scheme
+        '''
+        rxns = [reaction('A + A <=> B ; k1on k1off'), ]
+        c0 = {'A': 1.0, 'B': 0.0}
+        kin = {'k1on': 1.0, 'k1off': 1}
+        ts = np.linspace(0.0, 5, 10)
+
+        spcs, ct = concentration(rxns, c0, kin, ts)
+        print(spcs)
+        print(ct)
+        for spc in spcs:
+            plt.plot(ts, ct[:, spcs[spc]], '.-', label=spc)
+        # plt.plot(ts, ct[:,0], 'k.-', label='S')
+        # plt.plot(ts, ct[:,2], 'b.-', label='ES')
+        # plt.plot(ts, ct[:,3], 'r.-', label='P')
+        plt.grid()
+        plt.legend()
+        plt.show()
+        for i in range(len(ts)):
+            print('%8.3f %12.6e' % (ts[i], ct[i, 1]))
+        # test some concentrations
+
+        assert ct[spcs['A']][-1] <= 1.0e-9
+
+    def test_psensor_s1s2_concentration(self):
+        #  1: E + S1 <=> ES1; Kd_ES1
+        #  2: E + S2 <=> ES2; Kd_ES2
+        #  3: ES2 + S1 <=> ES1S2; Kd_ES1
+        #  4: ES1 + S2 <=> ES1S2; Kd_ES2
+        #  5: ES1S2 --> E + 2P + P2; kcat
+        #  6: R + P <=> RP; Kd_RP
+
+        rxns = [reaction('E + S1 <=> ES1           ; k1on k1off'),
+                reaction('E + S2 <=> ES2           ; k2on k2off'),
+                reaction('ES2 + S1 <=> ES1S2       ; k12on k12off'),
+                # reaction('ES1 + S1 <=> ES1S2       ; k21on k21off'), # TODO: that causes an interesting bug
+                reaction('ES1 + S2 <=> ES1S2       ; k21on k21off'),
+                reaction('ES1S2 -> E + 2P + P2     ; k3cat'),
+                reaction('R + P <=> RP             ; k4on k4off')]
+
+        for rxn in rxns:
+            print(rxn)
+
+        print
+
+        print_kinetic_equations(rxns)
+
+        print
+
+        # from B1   5.0e-4  1.00E-07    5.00E-04    5.00E-07
+        # convert to uM
+        c0 = {'E': 1.00E-01, 'S1': 5.1e2, 'S2': 5.00E2,
+              'R': 5.00E-01}  # , 'P' : 0.0}
+        Kd_ES1 = 1.0e3  # again in uM compatibility
+        Kd_ES2 = 1.0e3
+        Kd_RP = 2.7e0
+        kcat = 50.0
+        kr = 1000.0
+
+        kin = {'k1on':  kr / Kd_ES1, 'k1off': kr,
+               'k12on': kr / Kd_ES1, 'k12off': kr,
+               'k2on':  kr / Kd_ES2, 'k2off': kr,
+               'k21on': kr / Kd_ES2, 'k21off': kr,
+               'k3cat': kcat,
+               'k4on': kr / Kd_RP, 'k4off': kr}
+
+        ts = np.linspace(0.0, 100, 100)
+
+        species, ct = concentration(rxns, c0, kin, ts)
+
+        print(species)
+        print(ct)
+        for spc in species:  # ['P', 'X', 'RP']:
+            plt.plot(ts, ct[:, species[spc]], '.-', label=spc)
+
+        plt.grid()
+        plt.legend()
+        plt.show()
+
+        # for i in range(len(ts)):
+        #     print '%8.3f %12.6e' % (ts[i], ct[i, 1])
+
+        # test some concentrations
+        assert ct[-1, species['E']] >= 1.0e-9
+
+
+if __name__ == '__main__':
+    unittest.main()