Вы когда-нибудь задумывались, как действуют клетки — элементарные единицы живой материи? Я – программист, но одновременно увлекаюсь клеточной биологией. Поэтому я решил смоделировать работу простейшей клетки на TypeScript. Вообще, клетки невероятно сложны; по оценкам учёных, человеческая клетка в среднем содержит 100 триллионов атомов. По-прежнему очень мало известно о том, как все эти биомолекулы взаимодействуют в клетке, поэтому в точности смоделировать работу клетки невозможно.
Размышляя на эту тему, я нашёл статью Fundamental behaviors emerge from simulations of a living minimal cell (Фундаментальные виды поведения возникают на основе моделирования простейшей живой клетки). Опираясь на кинетические параметры, авторы статьи создали модель взаимодействия молекул и химических реакций между ними в простейшей известной клетке. Затем эта симуляция запускается, и на её основе можно наблюдать такие процессы как репликация ДНК, метаболизм и синтез белков.
Я хотел лучше понять, что именно сделали авторы этой статьи, поэтому воспроизвёл их работу в Typescript – и мой вариант симуляции можно легко выполнять в браузере, причём, для этого ничего не требуется устанавливать. В этой статье мы сосредоточимся на хорошо перемешиваемой модели, которая представляет целую клетку как целостную однородную смесь. В исходной статье хорошо перемешиваемая модель откалибрована и верифицирована на основе более продвинутой «пространственной» модели, использующей 3D-сетку, которая также учитывает движение молекул по клетке.
Минимальная клетка (JCVI-syn3A)
Чтобы смоделировать клетку, нужно узнать максимум информации о генах и белках у неё внутри, а также о том, какой метаболизм в ней происходит. В основном эти процессы до сих пор плохо изучены, особенно в крупных и сложных клетках, например, эукариотических. Поэтому учёные решили создать клетку с минимальной возможной сложностью, где сложность измеряется количеством генов.
Для этого взяли бактериальную клетку и стали удалять из неё гены до тех пор, пока клетка не утрачивала жизнеспособность. Таким образом было выяснено, какие гены абсолютно необходимы для выживания. При помощи этого метода удалось картировать жизненно важные и второстепенные гены. В конце концов, у них получилась минимальная клетка, которую они назвали JCVI-syn3A.
В такой максимально урезанной бактериальной клетке оказалось «всего» 493 гена. Сравните это число, например, с кишечной палочкой (E. Coli), у которой 4637 гена. Из 493 генов минимальной клетки у 20% функция не выяснена, но оставшейся части генов хватило, чтобы приступить к составлению модели этой клетки. Эта работа была описана в статье Kinetic Modeling of the Genetic Information Processes in a Minimal Cell. Полученные в итоге кинетические параметры легли в основу той симуляции, которую я построил.
Геномные данные
Все данные выложены в открытом доступе в репозитории на GitHub, который ведёт лаборатория Люти Шультен.
Весь геном сохранён в файле syn3A.gb. Это файл в формате GenBank, в котором хранится код ДНК и информация о генах и о тех белках, которые они кодируют. Например, видно, что ген, занимающий позиции 497075..498190, кодирует белок дигидрофолатсинтазу. Этот белок состоит из аминокислот MISV…NNKF, где M = метионин, I = изолейцин, т.д…
complement(497075..498190)
/gene="folC"
/locus_tag="JCVISYN3A_0823"
/product="Dihydrofolate synthase"
/protein_id="AVX54985.1"
/translation="MISVDQELFPINQRLNKEVVFDKVLDELNHPEEFLKVINVVGTN
GKGSTSFYLSKGLLKKYQKVGLFISPAFLYQNERIQINNTPISDNDLKSYLHKIDYLI
KKYQLLFFEIWTLIMILYFKDQKVDIVVCEAGIGGIKDTTSFLTNQLFSLCTSISYDH
MDILGNSIDEIIYNKINIAKPNTKLFISYDNLKYKDKINEQLINKNVELIYTDLYEDQ
IIYQQANKGLVKKVLEYLNIFQTNIFQLQPPLGRFTTIRTFPNHIIIDGAHNVDGINK
LIQSVKMLNKEFIVLYASITTKEYLKSLELLDQNFNEVYICEFNFIKSWSIDNIDHKN
KIKDWKKFLKNNTKNIIICGSLYFIPLVYNYLTNNKF"Кинетические параметры
Кинетические параметры хранятся в файлах SBTab, таких как этот, в котором описаны реакции, описывающие центральный метаболизм. В частности, здесь содержится информация об исходных концентрациях, реакциях между продуктами обмена, а также формулы и параметры, на основе которых вычислялась скорость протекания реакций.
Например, исходная концентрация АТФ (обозначенная как M_atp_c) составляет 3,6529 миллимоль (мM). Одна из реакций, в которой участвует АТФ, это M_atp_c + M_f6p_c <=> M_adp_c + M_fdp_c. В рамках этой реакции фруктозо-6-фосфат преобразуется во фруктозо-1,6-бифосфат, и для этого используется именно АТФ. Формула протекания скорости реакции приведена ниже:
( 0.001 / 1.0) * (( kcrg_R_PFK * ( ( ((keq_R_PFK)^(hco_R_PFK/2.)) * M_atp_c * M_f6p_c) - (((keq_R_PFK)^(-hco_R_PFK/2.)) * M_adp_c * M_fdp_c) ) )) / sqrt(kmc_R_PFK_M_atp_c * kmc_R_PFK_M_f6p_c * kmc_R_PFK_M_adp_c * kmc_R_PFK_M_fdp_c) / (( (1 + (M_atp_c/kmc_R_PFK_M_atp_c)) * (1 + (M_f6p_c/kmc_R_PFK_M_f6p_c)) ) + ( (1 + (M_adp_c/kmc_R_PFK_M_adp_c)) * (1 + (M_fdp_c/kmc_R_PFK_M_fdp_c)) ) - 1)Каждый из приведённых здесь параметров также определяется в конце файла. Например, kmc_R_PFK_M_atp_c = 0.117.
Основное химическое уравнение (CME)
Для моделирования генетических систем в статье используется основное химическое уравнение. Этот метод применяется при моделировании стохастических систем. В рамках уравнения CME моделируются следующие функции:
-
Транскрипция ДНК в мРНК
-
Трансляция мРНК в белки
-
Разложение мРНК
-
Репликация ДНК
-
Аминоацилирование тРНК
Термин «основное химическое уравнение» довольно громкий, но на самом деле оно сводится к достаточно простому алгоритму, который называется «алгоритм Гиллеспи». В любой момент времени t мы вычисляем «интенсивность» каждой реакции, а этот показатель вычисляется как произведение скорости реакции на концентрацию реагентов. В качестве шага развития модели мы выбираем произвольный промежуток времени, исходя из общей интенсивности. Чем выше общая интенсивность, тем короче в среднем будет такой шаг.
Δt = (1 / totalPropensity) * ln(1 / randomNumber)При выборе случайной реакции шанс на выбор конкретной реакции тем выше, чем выше её интенсивность. Затем эта реакция «выполняется», а концентрации реагентов соответствующим образом корректируются. Вот и всё, просто повторяем этот процесс, пока не выйдем на целевое время и не вернём результирующие значения концентрации.
Итак, на практике происходит следующее: мы постоянно выбираем случайные гены для последующей их трансляции в мРНК, а далее случайные нити мРНК для трансляции в белки. В результате от клетки к клетке возникает изменчивость. Некоторые клетки растут быстрее, либо в них быстрее начинается репликация ДНК. При многократных прогонах программы вы всякий раз получаете разные результаты, точно как происходило бы и в реальном мире.
Обыкновенное дифференциальное уравнение (ODE)
Метаболические реакции вычисляются детерминированно (а не случайным образом) при помощи обыкновенного дифференциального уравнения (ODE). Среди различных систем, к которым оно применимо – следующие:
-
Центральный метаболизм
-
Метаболизм липидов
-
Метаболизм нуклеотидов
-
Метаболизм аминокислот
-
Метаболизм кофакторов
-
Транспортные реакции
Дифференциальные уравнения часто используются в компьютерных программах. Например, создавая игру, мы зачастую имитируем движение путём постоянного продвижения позиции путём умножения скорости на Δt. В таком случае скорость является производной позиции.
В нашей биологической симуляции уровни концентрации увеличиваются путём умножения скоростей реакции на Δt. Для каждого метаболита учитываются все реакции, при которых он расходуется и при которых он производится. В таком случае производная функция очень велика, и при этом она очень «жёсткая». Функция называется жёсткой, если она очень нестабильна и быстро устремляется в бесконечность, если её как следует не интегрировать.
Чтобы решить это дифференциальное уравнение, необходимо воспользоваться специальным методом интегрирования, который называется LSODA. Это алгоритм, постоянно проверяющий ошибки и корректирующий длительность шага, чтобы поддерживать стабильность интегрирования. Исходная версия алгоритма LSODA была написана в 1983 году на языке программирования Fortran. Но, к счастью, существует проект libsoda, в котором содержится версия этого алгоритма на C.
При помощи инструмента Emscripten этот код можно скомпилировать в Web Assembly (WASM), так, чтобы его можно было выполнять в браузере. Здесь возникает ещё одно преимущество – выполнение идёт очень быстро, поскольку Web Assembly работает непосредственно на ЦП и, соответственно, сильно превосходит по скорости JavaScript.
Коммуникация между CME и ODE
Итак, теперь у нас две отдельные системы: CME для описания генетических процессов и ODE для моделирования метаболизма. Чтобы скомбинировать их, между ними нужно наладить обмен информацией в каком-либо виде.
Наша симуляция основана на CME, эта процедура инициализируется вместе с исходными концентрациями. Спустя каждую секунду симуляции, частицы CME преобразуются в значения концентрации ODE и записываются в ODE. Затем ODE интегрируется в течение 1 секунды, а результирующие концентрации преобразуются обратно в частицы и записываются в CME. Эта последовательность повторяется вплоть до завершения симуляции.
Результаты
Ниже приведены результаты одной из симуляции. Сравнив их с результатами из статьи, отмечаем сходство значений в теории и на практике, тем самым подтверждается корректность симуляции. После множества прогонов достигается максимальный объём в 6,7e-17 л за срок 60 — 70 минут.
Ниже показаны некоторые пулы доступных метаболитов.
А на этой последней схеме показан расход АТФ при некоторых видах деятельности. Как видите, репликация ДНК начинается примерно на 36-й минуте и завершается по истечении 85 минут в ходе симуляции. От раза к разу эти показатели варьируются, репликация иногда может начаться даже на 6-й минуте. Иногда даже успеваем увидеть два полных цикла репликации генома.
После нескольких шагов оптимизации производительности результат также весьма хорош. На пятилетнем макбуке мне удалось завершить весь цикл за 20 минут. Всё происходит в браузере и всего в одном потоке, поэтому, если задействовать веб-воркеры, то можно будет распараллелить сразу несколько симуляций.
Что далее
В описанном виде эта симуляция не учитывает какую-либо клеточную механику. Поэтому деление клеток моделируется на бумаге по правилам. Думаю, именно клеточное деление – наиболее актуальная тема для ближайших исследований. Ещё одна цель – масштабировать этот эксперимент, чтобы можно было охватить в его рамках более сложные клетки, например, человеческие. В них было бы очень удобно и интересно изучать эффекты от приёма лекарств.
Следующий фронтир – полная молекулярная симуляция целой клетки. В симуляции такого рода учитывается каждый атом (или небольшая группа атомов), и такой опыт гораздо ближе к реальности. Первичные исследования на эту тему уже были проведены и, в частности, описаны в статье Molecular dynamics simulation of an entire cell. На данный момент для такого проекта недостаёт как научных данных, так и программных функций, поэтому до реализации подобного проекта ещё далеко.
Демонстрация и код
Симуляцию можно самостоятельно попробовать здесь.
Весь исходный код к этому проекту выложен на гитхабе автора.
