genODE

Synopsis: Generate Python ODE script from a given model specification file.

Usage: python astools.py genODE [option]

where [option] can be

• modelfile: Name of model specification file in models folder. This assumes that the model file is not in models folder.
• type: Type of model specification file. Allowable types are 'ASM' (AdvanceSyn Model Specification). Default = 'ASM'.
• solver: Type of solver to use. Allowable types are 'Euler', 'Heun' (Runge-Kutta 2nd method or Trapezoidal), 'RK3' (third order Runge-Kutta), 'RK4' (fourth order Runge-Kutta), 'RK38' (fourth order Runge-Kutta method, 3/8 rule), 'CK4' (fourth order Cash-Karp), 'CK5' (fifth order Cash-Karp), 'RKF4' (fourth order Runge-Kutta-Fehlberg), 'RKF5' (fifth order Runge-Kutta-Fehlberg), 'DP4' (fourth order Dormand-Prince), and 'DP5' (fifth order Dormand-Prince). Default = 'RK4'.
• timestep: Time step interval for simulation. Default = 1.0.
• endtime: Time to end simulation - the simulation will run from 0 to end time. Default = 21600.
• lowerbound: Define lower boundary of objects. For example, "1;2" means that when the value of the object hits 1, it will be bounced back to 2. Default = 0;0; that is, when the value of the object goes to negative, it will be bounced back to zero.
• upperbound: Define upper boundary of objects. For example, "10;9" means that when the value of the object hits 1, it will be pushed down to 9. Default = 1e-3;1e-3; that is, when the value of the object above 1e-3, it will be pushed back to 1e-3.
• odefile: Python ODE script file to write out. This file will be written into odescript folder. Default = odescript.py.

For example:

python genODE \
    --modelfile=models/asm/glycolysis.modelspec \
    --mtype=ASM \
    --solver=RK4 \
    --timestep=1 \
    --endtime=21600 \
    --lowerbound=0;0 \
    --upperbound=1e-3;1e-3 \
    --odefile=glycolysis.py

Working example:

D:\Dropbox\MyProjects\astoolkit>python astools.py genODE --modelfile=models/asm/glycolysis.modelspec --mtype=ASM --solver=RK4 --timestep=1 --endtime=21600 --lowerbound=0;0 --upperbound=1e-3;1e-3 --odefile=glycolysis.py
''' -------------------------------------
Python ODE script generated by ASModeller
(a package in AdvanceSynToolKit)

Date Time: 2020-5-31 3:32:32:358011

name: glycolysis
author: Maurice Ling
------------------------------------- '''

def glucose(t, y):
    r1 = (1e-6 * y[19] * y[0])/(1e-6 + (y[19] * y[0]))
    return (0) - (r1)

def atp(t, y):
    r7 = (1e-6 * y[25] * y[13])/(1e-6 + (y[25] * y[13]))
    r10 = (1e-6 * y[28] * y[16])/(1e-6 + (y[28] * y[16]))
    r1 = (1e-6 * y[19] * y[0])/(1e-6 + (y[19] * y[0]))
    r3 = (1e-6 * y[21] * y[9])/(1e-6 + (y[21] * y[9]))
    return (r7 + r10) - (r1 + r3)

def adp(t, y):
    r1 = (1e-6 * y[19] * y[0])/(1e-6 + (y[19] * y[0]))
    r3 = (1e-6 * y[21] * y[9])/(1e-6 + (y[21] * y[9]))
    r7 = (1e-6 * y[25] * y[13])/(1e-6 + (y[25] * y[13]))
    r10 = (1e-6 * y[28] * y[16])/(1e-6 + (y[28] * y[16]))
    return (r1 + r3) - (r7 + r10)

def proton(t, y):
    r1 = (1e-6 * y[19] * y[0])/(1e-6 + (y[19] * y[0]))
    r3 = (1e-6 * y[21] * y[9])/(1e-6 + (y[21] * y[9]))
    r6 = (1e-6 * y[24] * y[11])/(1e-6 + (y[24] * y[11]))
    r10 = (1e-6 * y[28] * y[16])/(1e-6 + (y[28] * y[16]))
    return (r1 + r3 + r6) - (r10)

def nad(t, y):
    r6 = (1e-6 * y[24] * y[11])/(1e-6 + (y[24] * y[11]))
    return (0) - (r6)

def pi(t, y):
    r6 = (1e-6 * y[24] * y[11])/(1e-6 + (y[24] * y[11]))
    return (0) - (r6)

def nadh(t, y):
    r6 = (1e-6 * y[24] * y[11])/(1e-6 + (y[24] * y[11]))
    return (r6) - (0)

def water(t, y):
    r9 = (1e-6 * y[27] * y[15])/(1e-6 + (y[27] * y[15]))
    return (r9) - (0)

def g6p(t, y):
    r1 = (1e-6 * y[19] * y[0])/(1e-6 + (y[19] * y[0]))
    r2 = (1e-6 * y[20] * y[8])/(1e-6 + (y[20] * y[8]))
    return (r1) - (r2)

def f6p(t, y):
    r2 = (1e-6 * y[20] * y[8])/(1e-6 + (y[20] * y[8]))
    r3 = (1e-6 * y[21] * y[9])/(1e-6 + (y[21] * y[9]))
    return (r2) - (r3)

def f16p(t, y):
    r3 = (1e-6 * y[21] * y[9])/(1e-6 + (y[21] * y[9]))
    r4 = (1e-6 * y[22] * y[10])/(1e-6 + (y[22] * y[10]))
    return (r3) - (r4)

def gadp(t, y):
    r4 = (1e-6 * y[22] * y[10])/(1e-6 + (y[22] * y[10]))
    r5 = (1e-6 * y[23] * y[12])/(1e-6 + (y[23] * y[12]))
    r6 = (1e-6 * y[24] * y[11])/(1e-6 + (y[24] * y[11]))
    return (r4 + r5) - (r6)

def dhap(t, y):
    r4 = (1e-6 * y[22] * y[10])/(1e-6 + (y[22] * y[10]))
    r5 = (1e-6 * y[23] * y[12])/(1e-6 + (y[23] * y[12]))
    return (r4) - (r5)

def bpg13(t, y):
    r6 = (1e-6 * y[24] * y[11])/(1e-6 + (y[24] * y[11]))
    r7 = (1e-6 * y[25] * y[13])/(1e-6 + (y[25] * y[13]))
    return (r6) - (r7)

def pg3(t, y):
    r7 = (1e-6 * y[25] * y[13])/(1e-6 + (y[25] * y[13]))
    r8 = (1e-6 * y[26] * y[14])/(1e-6 + (y[26] * y[14]))
    return (r7) - (r8)

def pg2(t, y):
    r8 = (1e-6 * y[26] * y[14])/(1e-6 + (y[26] * y[14]))
    r9 = (1e-6 * y[27] * y[15])/(1e-6 + (y[27] * y[15]))
    return (r8) - (r9)

def pep(t, y):
    r9 = (1e-6 * y[27] * y[15])/(1e-6 + (y[27] * y[15]))
    r10 = (1e-6 * y[28] * y[16])/(1e-6 + (y[28] * y[16]))
    return (r9) - (r10)

def pyr(t, y):
    r10 = (1e-6 * y[28] * y[16])/(1e-6 + (y[28] * y[16]))
    return (r10) - (0)

def hk_rna(t, y):
    e1_rna = 1e-7
    return (e1_rna) - (0)

def hk(t, y):
    e1 = 1e-5 * y[18]
    return (e1) - (0)

def pgi(t, y):
    return (0) - (0)

def pfk(t, y):
    return (0) - (0)

def aldo(t, y):
    return (0) - (0)

def tpi(t, y):
    return (0) - (0)

def gapdh(t, y):
    return (0) - (0)

def pkg(t, y):
    return (0) - (0)

def pgm(t, y):
    return (0) - (0)

def eno(t, y):
    return (0) - (0)

def pk(t, y):
    return (0) - (0)

def boundary_checker(y, boundary, type):
    '''
    Private function - called by ODE solvers to perform boundary checking
    of variable values and reset them to specific values if the original
    values fall out of the boundary values.

    Boundary parameter takes the form of a dictionary with variable number
    as key and a list of [<boundary value>, <value to set if boundary is
    exceeded>]. For example, the following dictionary for lower boundary
    (type = 'lower') {'1': [0.0, 0.0], '5': [2.0, 2.0]} will set the lower
    boundary of variable y[0] and [5] to 0.0 and 2.0 respectively. This
    also allows for setting to a different value - for example, {'1': [0.0,
    1.0]} will set variable y[0] to 2.0 if the original y[0] value is
    negative.

    @param y: values for variables
    @type y: list
    @param boundary: set of values for boundary of variables
    @type boundary: dictionary
    @param type: the type of boundary to be checked, either 'upper' (upper
    boundary) or 'lower' (lower boundary)
    '''
    for k in list(boundary.keys()):
        if y[int(k)] < boundary[k][0] and type == 'lower':
            y[int(k)] = boundary[k][1]
        if y[int(k)] > boundary[k][0] and type == 'upper':
            y[int(k)] = boundary[k][1]
    return y

def RK4(funcs, x0, y0, step, xmax, nonODEfunc=None,
        lower_bound=None, upper_bound=None,
        overflow=1e100, zerodivision=1e100):
    '''
    Generator to integrate a system of ODEs, y' = f(x, y), using fourth
    order Runge-Kutta method.

    A function (as nonODEfunc parameter) can be included to modify one or
    more variables (y0 list). This function will not be an ODE (not a
    dy/dt). This can be used to consolidate the modification of one or
    more variables at each ODE solving step. For example, y[0] = y[1] / y[2]
    can be written as

    >>> def modifying_function(y, step):
            y[0] = y[1] / y[2]
            return y

    This function must take 'y' (variable list) and 'step' (time step) as
    parameters and must return 'y' (the modified variable list). This
    function will execute before boundary checking at each time step.

    Upper and lower boundaries of one or more variable can be set using
    upper_bound and lower_bound parameters respectively. These parameters
    takes the form of a dictionary with variable number as key and a list
    of [<boundary value>, <value to set if boundary is exceeded>]. For
    example, the following dictionary for lower boundary {'1': [0.0, 0.0],
    '5': [2.0, 2.0]} will set the lower boundary of variable y[0] and y[5]
    to 0.0 and 2.0 respectively. This also allows for setting to a different
    value - for example, {'1': [0.0, 1.0]} will set variable y[0] to 2.0 if
    the original y[0] value is negative.

    @param funcs: system of differential equations
    @type funcs: list
    @param x0: initial value of x-axis, which is usually starting time
    @type x0: float
    @param y0: initial values for variables
    @type y0: list
    @param step: step size on the x-axis (also known as step in calculus)
    @type step: float
    @param xmax: maximum value of x-axis, which is usually ending time
    @type xmax: float
    @param nonODEfunc: a function to modify the variable list (y0)
    @type nonODEfunc: function
    @param lower_bound: set of values for lower boundary of variables
    @type lower_bound: dictionary
    @param upper_bound: set of values for upper boundary of variables
    @type upper_bound: dictionary
    @param overflow: value (usually a large value) to assign in event of
    over flow error (usually caused by a large number) during integration.
    Default = 1e100.
    @type overflow: float
    @param zerodivision: value (usually a large value) to assign in event
    of zero division error, which results in positive infinity, during
    integration. Default = 1e100.
    @type zerodivision: float
    '''
    yield [x0] + y0
    def solver(funcs, x0, y0, step):
        n = len(funcs)
        f1, f2, f3, f4 = [0]*n, [0]*n, [0]*n, [0]*n
        y1 = [0]*n
        for i in range(n):
            try: f1[i] = funcs[i](x0, y0)
            except TypeError: pass
            except ZeroDivisionError: f1[i] = zerodivision
            except OverflowError: f1[i] = overflow
        for j in range(n):
            y1[j] = y0[j] + (0.5*step*f1[j])
        for i in range(n):
            try: f2[i] = funcs[i]((x0+(0.5*step)), y1)
            except TypeError: pass
            except ZeroDivisionError: f2[i] = zerodivision
            except OverflowError: f2[i] = overflow
        for j in range(n):
            y1[j] = y0[j] + (0.5*step*f2[j])
        for i in range(n):
            try: f3[i] = funcs[i]((x0+(0.5*step)), y1)
            except TypeError: pass
            except ZeroDivisionError: f3[i] = zerodivision
            except OverflowError: f3[i] = overflow
        for j in range(n):
            y1[j] = y0[j] + (step*f3[j])
        for i in range(n):
            try: f4[i] = funcs[i]((x0+step), y1)
            except TypeError: pass
            except ZeroDivisionError: f4[i] = zerodivision
            except OverflowError: f4[i] = overflow
        for i in range(n):
            try: y1[i] = y0[i] + (step * \
                    (f1[i] + (2.0*f2[i]) + (2.0*f3[i]) + f4[i]) / 6.0)
            except TypeError: pass
            except ZeroDivisionError: y1[i] = zerodivision
            except OverflowError: y1[i] = overflow
        return y1
    while x0 < xmax:
        y1 = solver(funcs, x0, y0, step)
        if nonODEfunc:
            y1 = nonODEfunc(y1, step)
        if lower_bound:
            y1 = boundary_checker(y1, lower_bound, 'lower')
        if upper_bound:
            y1 = boundary_checker(y1, upper_bound, 'upper')
        y0 = y1
        x0 = x0 + step
        yield [x0] + y0

ODE = list(range(29))
ODE[0] = glucose
ODE[1] = atp
ODE[2] = adp
ODE[3] = proton
ODE[4] = nad
ODE[5] = pi
ODE[6] = nadh
ODE[7] = water
ODE[8] = g6p
ODE[9] = f6p
ODE[10] = f16p
ODE[11] = gadp
ODE[12] = dhap
ODE[13] = bpg13
ODE[14] = pg3
ODE[15] = pg2
ODE[16] = pep
ODE[17] = pyr
ODE[18] = hk_rna
ODE[19] = hk
ODE[20] = pgi
ODE[21] = pfk
ODE[22] = aldo
ODE[23] = tpi
ODE[24] = gapdh
ODE[25] = pkg
ODE[26] = pgm
ODE[27] = eno
ODE[28] = pk

y = list(range(29))
y[0] = 1e-5    # glucose : D-glucose
y[1] = 1e-5    # atp : adenosine-triphosphate
y[2] = 1e-5    # adp : adenosine-diphosphate
y[3] = 1e-5    # proton : proton
y[4] = 1e-5    # nad : NAD
y[5] = 1e-5    # pi : phosphate
y[6] = 1e-5    # nadh : NADH
y[7] = 1e-5    # water : water
y[8] = 1e-9    # g6p : a-D-Glucose-6-phosphate
y[9] = 1e-9    # f6p : b-D-Fructose-6-phosphate
y[10] = 1e-9    # f16p : b-D-Fructose-1,6-phosphate
y[11] = 1e-9    # gadp : D-glyceraldehyde 3-phosphate
y[12] = 1e-9    # dhap : Dihydroxyacetone phosphate
y[13] = 1e-9    # bpg13 : D-1,3-bisphosphoglycerate
y[14] = 1e-9    # pg3 : 3-phosphoglycerate
y[15] = 1e-9    # pg2 : 2-phosphoglycerate
y[16] = 1e-9    # pep : phosphoenolpyruvate
y[17] = 1e-9    # pyr : pyruvate
y[18] = 0    # hk.rna : hexokinase
y[19] = 0    # hk : hexokinase
y[20] = 1e-6    # pgi : Phosphoglucose isomerase
y[21] = 1e-6    # pfk : phosphofructokinase
y[22] = 1e-6    # aldo : fructose-bisphosphate aldolase
y[23] = 1e-6    # tpi : triosephosphate isomerase
y[24] = 1e-6    # gapdh : glyceraldehyde phosphate dehydrogenase
y[25] = 1e-6    # pkg : phosphoglycerate kinase
y[26] = 1e-6    # pgm : phosphoglycerate mutase
y[27] = 1e-6    # eno : enolase
y[28] = 1e-6    # pk : pyruvate kinase

labels = ['time', 'glucose', 'atp', 'adp', 'proton', 'nad', 'pi', 'nadh', 'water', 'g6p', 'f6p', 'f16p', 'gadp', 'dhap', 'bpg13', 'pg3', 'pg2', 'pep', 'pyr', 'hk_rna', 'hk', 'pgi', 'pfk', 'aldo', 'tpi', 'gapdh', 'pkg', 'pgm', 'eno', 'pk']

lowerbound = {'0': [0.0, 0.0], '1': [0.0, 0.0], '2': [0.0, 0.0], '3': [0.0, 0.0], '4': [0.0, 0.0], '5': [0.0, 0.0], '6': [0.0, 0.0], '7': [0.0, 0.0], '8': [0.0, 0.0], '9': [0.0, 0.0], '10': [0.0, 0.0], '11': [0.0, 0.0], '12': [0.0, 0.0], '13': [0.0, 0.0], '14': [0.0, 0.0], '15': [0.0, 0.0], '16': [0.0, 0.0], '17': [0.0, 0.0], '18': [0.0, 0.0], '19': [0.0, 0.0], '20': [0.0, 0.0], '21': [0.0, 0.0], '22': [0.0, 0.0], '23': [0.0, 0.0], '24': [0.0, 0.0], '25': [0.0, 0.0], '26': [0.0, 0.0], '27': [0.0, 0.0], '28': [0.0, 0.0], }

upperbound = {'0': [0.001, 0.001], '1': [0.001, 0.001], '2': [0.001, 0.001], '3': [0.001, 0.001], '4': [0.001, 0.001], '5': [0.001, 0.001], '6': [0.001, 0.001], '7': [0.001, 0.001], '8': [0.001, 0.001], '9': [0.001, 0.001], '10': [0.001, 0.001], '11': [0.001, 0.001], '12': [0.001, 0.001], '13': [0.001, 0.001], '14': [0.001, 0.001], '15': [0.001, 0.001], '16': [0.001, 0.001], '17': [0.001, 0.001], '18': [0.001, 0.001], '19': [0.001, 0.001], '20': [0.001, 0.001], '21': [0.001, 0.001], '22': [0.001, 0.001], '23': [0.001, 0.001], '24': [0.001, 0.001], '25': [0.001, 0.001], '26': [0.001, 0.001], '27': [0.001, 0.001], '28': [0.001, 0.001], }

model = RK4(ODE, 0.0, y, 1, 21600, None, lowerbound, upperbound)

(py37) D:\Dropbox\MyProjects\astoolkit>