Питон / Биомолекулярная физика - Попытка закодировать простое стохастическое моделирование системы, проявляющей условное поведение!

Question

* отредактировано 6/17/10

Я пытаюсь понять, как улучшить мой код (сделать его более питоническим). Кроме того, я заинтересован в написании более интуитивных «условий», которые описывают сценарии, которые являются обычным явлением в биохимии. Условные критерии в приведенной ниже программе я объяснил в ответе № 2, но код меня не устраивает - он работает нормально, но неочевиден и не так прост для реализации в более сложных условных сценариях. Идеи приветствуются. Комментарии / критика приветствуются. Первая публикация опыта @ stackoverflow - пожалуйста, прокомментируйте этикет, если это необходимо.

Код генерирует список значений, которые являются решением для следующего упражнения:

"На выбранном вами языке программирования реализуйте алгоритм первой реакции Гиллеспи, чтобы изучить временное поведение реакции A ---> B, в которой переход от A к B может иметь место только в том случае, если другое соединение, C, является присутствует и где C динамически взаимопревращается с D, как смоделировано в сети Петри ниже. Предположим, что в начале реакции присутствуют 100 молекул A, 1 C и нет B или D. В качестве значения kAB установлено значение 0,1 с. -1 и kCD и kDC до 1,0 с-1. Имитация поведения системы в течение 100 с. "

def sim():
    # Set the rate constants for all transitions
    kAB = 0.1
    kCD = 1.0
    kDC = 1.0

    # Set up the initial state
    A = 100
    B = 0
    C = 1
    D = 0

    # Set the start and end times
    t = 0.0
    tEnd = 100.0

    print "Time\t", "Transition\t", "A\t", "B\t", "C\t", "D"

    # Compute the first interval
    transition, interval = transitionData(A, B, C, D, kAB, kCD, kDC)
    # Loop until the end time is exceded or no transition can fire any more
    while t <= tEnd and transition >= 0:
        print t, '\t', transition, '\t', A, '\t', B, '\t', C, '\t', D
        t += interval
        if transition == 0:
            A -= 1
            B += 1
        if transition == 1:
            C -= 1
            D += 1
        if transition == 2:
            C += 1
            D -= 1
        transition, interval = transitionData(A, B, C, D, kAB, kCD, kDC)


def transitionData(A, B, C, D, kAB, kCD, kDC):
    """ Returns nTransition, the number of the firing transition (0: A->B,
    1: C->D, 2: D->C), and interval, the interval between the time of
    the previous transition and that of the current one. """
    RAB = kAB * A * C
    RCD = kCD * C
    RDC = kDC * D
    dt = [-1.0, -1.0, -1.0]
    if RAB > 0.0:
        dt[0] = -math.log(1.0 - random.random())/RAB
    if RCD > 0.0:
        dt[1] = -math.log(1.0 - random.random())/RCD
    if RDC > 0.0:
        dt[2] = -math.log(1.0 - random.random())/RDC
    interval = 1e36
    transition = -1
    for n in range(len(dt)):
        if dt[n] > 0.0 and dt[n] < interval:
            interval = dt[n]
            transition = n
    return transition, interval       


if __name__ == '__main__':
    sim()

Answer 1

не уверен, что вы видели это.

http://stompy.sourceforge.net/html/userguide_doc.html

Я работаю над подобными вещами, и я случайно обнаружил это недавно.

Answer 2

Информация о математике, стоящей за простым стохастическим моделированием химических rxns:

Как правило, подобные процессы имитируются как дискретные события, причем каждое событие происходит с вероятностью «P» при заданной константе скорости «k» иколичество возможных событий «n» в интервале времени «dt»: P = 1-e ** (- k dt n).Здесь мы пренебрегаем фактическим временем каждого события (~ 0) и вместо этого сосредотачиваемся на интервале времени, в котором происходит событие.Любой, кто знаком с N, выберет K задач / испытаний Бернулли, оценит наличие 1 / e, например, когда N = K и N-> oo, вероятность того, что не будет выбран определенный элемент из N, приближается к 1 / e.Следовательно, в стохастической химической реакции (первый порядок) вероятность того, что молекула не подвергнется реакции (не будет выбрана), равна некоторой степени 1 / e ... эта мощность зависит от временного интервала и константы скорости, а также отколичество молекул и константа скорости в вопросе.И наоборот, 1- (1 / e) ^ xyz дает вероятность того, что любая конкретная молекула будет реагировать (будет выбрана).

С точки зрения симуляции, было бы логично разделить наш общий временной интервал на все меньшие интервалы и использовать генератор случайных чисел, чтобы предсказать, произошло ли событие в данном временном интервале - например, если мы разделили dt дляединое целое даже на 10 меньших интервалов, число от 0 до 0,1 будет означать, что произошло событие, а число от 0,1 до 1,0 будет указывать, что это не произошло.Однако существует неопределенность относительно того, когда именно произошло событие, - поэтому мы должны сократить интервалы - это быстро превращается в проигрышную битву, поскольку неопределенность сохраняется с этим методом.

Решение этой проблемы состоит в том, чтобы взять натуральный логарифм (здесь 'ln', log () по умолчанию в py) обеих сторон вышеприведенного уравнения и решить для dt, что дает dt = (-ln(1-Р)) / (К * п).Вероятность P затем генерируется случайным образом, давая определенное значение dt для каждого события.

Answer 3

**** Редактировать **** Я изначально объяснил это неправильно !!!!Следующее верно, хотя - Джастин, эта программа использует умные критерии для «взвешивания» каждого события.Все значения RAB, RCD и RDC получают параметр истина / ложь путем умножения kAB, kCD и kDC на C или D, которые в этом случае могут быть равны единице или нулю.Нулевое значение для D, и, следовательно, RDC предотвратит отрисовку dt [2] в

для n в диапазоне (len (dt)): если dt [n]> 0,0 и dt [n] <интервал: </p>

оператор.Кроме того, следующее -

    if transition == 1:
        C -= 1
        D += 1
    if transition == 2:
        C += 1
        D -= 1

диктует, что когда происходит событие C-> D (переход 1), следующее событие обязательно должно быть D-> C (переход 2), поскольку из трех значений вdt [], только dt [1] отлично от нуля и, таким образом, соответствует вышеупомянутым критериям.Итак, как мы оцениваем вероятность того, что переход 0 или переход 1 произойдет?Это немного сложно, но присуще следующим строкам:

interval = 1e36
transition = -1
for n in range(len(dt)):  
    if dt[n] > 0.0 and dt[n] < interval:            
        interval = dt[n]            
        transition = n            
return transition, interval   

"для n в диапазоне (len (dt)):" возвращает все значения списка dt [].Следующая строка определяет критерии, которые должны быть выполнены, а именно, что каждое значение должно быть больше 0 и меньше интервала.Для перехода 0 интервал равен 1е36 (который должен представлять бесконечность).Суть в том, что интервал затем устанавливается на переход 0, поэтому для второго значения в dt [], переход 1, критерии утверждают, что он должен быть меньше, чем dt для перехода 0, чтобы происходить ... или другими словамичто это должно было произойти быстрее, чтобы произошло вообще, что согласуется с химической логикой.Больше всего меня беспокоит то, что накопленные значения t, определяемые линией «t + = интервал», могут быть не совсем справедливыми ... потому что, поскольку запуск t1 не зависит от запуска t0, запуск t0 и принятие, скажем, 0,1 секунды, не должныисключить t1 от использования того же .1 сек для стрельбы ... но я работаю над исправлением для этого ... предложения приветствуются!Это подробная распечатка из сценария, включая запуск перехода 1 и 2:

Переход времени ABCD

dt0= 0.0350693547214
dt1= 2.26710773787
interval= 1e+36
dt= 0.0350693547214
transition= 0

0.0 0 100 0 1 0

dt0= 0.000339596342313
dt1= 0.21083283004
interval= 1e+36
dt= 0.000339596342313 
transition= 0

0,0350693547214 0 99 1 1 0

dt0= 0.0510125874767
dt1= 1.26127048627
interval= 1e+36 
dt= 0.0510125874767
transition= 0

0,0354089510637 0 98 2 1 0

dt0= 0.0809691957218
dt1= 0.593246425076
interval= 1e+36
dt= 0.0809691957218
transition= 0

0,0864215385404 0 97 3 1 0

dt0= 0.00205040633531
dt1= 1.70623338677
interval= 1e+36
dt= 0.00205040633531
transition= 0

0,167390734262 0 964 1 0

dt0= 0.106140534256
dt1= 0.0915160378053
interval= 1e+36
dt= 0.106140534256
transition= 0
interval= 0.106140534256
dt= 0.0915160378053
transition= 1

0.169441140598 1 95 5 1 0

dt2= 0.0482892532952
interval= 1e+36
dt= 0.0482892532952
transition= 2

0.260957178403 2 95 5 0 1

dt0= 0.112545351421
dt1= 1.84936696832
interval= 1e+36
dt= 0.112545351421
transition= 0

0.309246431698 0 95 5 1 0

Джастин, я не уверен, что вы имеете в виду, когда dt [2] меньше 1e36, заставляя его «оставаться» на переходе 2?Этого не происходит из-за оператора

if transition == 2:
            C += 1
            D -= 1

.Кто-нибудь знает более прямой способ сделать это

Ха-ха, позвольте разгореться пламени - вы, ребята, просто потрясающие - я очень ценю обратную связь!Stackoverflow является оооочень законным.

Answer 4

Я не знаю алгоритм Гиллеспи, но я предполагаю, что вы проверили, что программа сходится к правильному равновесию. Поэтому я интерпретирую ваши вопросы как

«Вот работающая физическая программа, как я могу сделать ее более питонической»

Вероятно, было бы более питоничным сделать что-то вроде следующего

R = [ kAB * A * C,  kCD * C, kAB * A * C]
dt = [(-math.log(1-random.random())/x,i) for i,x in enumerate(R) if x > 0]
if dt:
     interval,transition = min(dt)
else:
     transition = None

Если вы хотите использовать python в физике, тогда я предлагаю вам изучить numpy. Потому что NumPy быстрее для многих проблем. Итак, вот некоторые непроверенные части решения для numpy. Добавьте следующее в шапку вашей программы

from numpy import log, array, isinf, zeros
from numpy.random import rand

Затем вы можете заменить внутреннюю TransitionData чем-то вроде следующего

R = array([ kAB * A * C,  kCD * C, kAB * A * C])
dt = -log(1-rand(3))/R
transition = dt.argmin()
interval = dt[transition]
if isinf(interval):
    transition = None

Я не знаю, было бы более питонным выдавать исключение StopIteration вместо возврата None, но это деталь.

Вы также должны хранить ваши концентрации в единой структуре данных. Если вы используете NumPy, то я предлагаю вам использовать массив. Аналогично, yoy может использовать массив dABCD для хранения изменений в концентрации (возможно, вы найдете более подходящие имена переменных). Добавьте следующий код вне цикла

ABCD = array([A,B,C,D])

dABCD = zeros(3,4)
dABCD[0,0] = -1#A decreases when transition == 0
dABCD[0,1] = 1 #B increases when transition == 0
dABCD[1,2] = -1#C decreases when transition == 1
dABCD[1,3] = 1 #D increases when transition == 1
.....      etc

Теперь вы можете заменить основной цикл на что-то вроде следующего

while t <= tEnd:
       print t, '\t', transition, '\t', ABCD
       transition, interval = transitionData(ABCD, kAB, kCD, kDC)
       if transition != None:
            t += interval
           ABCD += dABCD[transition,:]
       else:
           break;
else:
       print "Warning: Stopping due to timeout. The system did not equilibrate"

Вероятно, есть еще что сделать. В качестве примера dABCD, вероятно, должен быть разреженным массивом, но я надеюсь, что эти идеи могут быть началом.