Стажировка проходила под непосредственным руководством д-ра Омара Валссона (руководитель группы "Методы эффективного сэмплирования для мягкой материи", отдел теории полимеров). Он является одним из разработчиков библиотеки с открытым кодом Plumed (https://www.plumed.org/), посредством которой методы эффективного сэмплирования интегрируются с известными пакетами программ молекулярного моделирования (Gromacs, Amber, LAMPPS и т. д.).
Для работы над проектом по созданию противоопухолевых пептидов был предложен следующий план исследований:
1) выбор силового поля для проведения моделирования
2) расчет поверхности свободной энергии выбранного лиганда к интегрину αVβ6 в физиологическом растворе: определение вероятностей всех конформаций, которые лиганд может принимать, определение устойчивости наиболее заселенной конформации, а также той конформации, которая связывается с интегриновым рецептором.
3) расчет динамических характеристик выбранного лиганда к интегрину αVβ6 в физиологическом растворе (значений 15N{1H} гетероядерного NOE, 15N T1 и T2 времен релаксации), сравнение этих характеристик с имеющимися в литературе экспериментальными данными и построение модели динамики лиганда в растворе. Этот шаг необходим, так как на связывание лиганда с белком влияет не только термодиамическая стабильность отдельных конформаций лиганда, но и его динамическая подвижность. Жесткий лиганд может оказаться менее селективным, чем гибкий, или наоборот. А кроме того только сравнение термодинамических и динамических расчетных данных с экспериментальными позволяет окончательно сказать, насколько правильно было выбрано силовое поле.
4) расчет поверхности свободной энергии и динамических характеристик предложенного химерного пептида, состоящего из лиганда, линкера и «боеголовки». Из сравнения расчетов отдельного лиганда и химерного пептида можно определить, каково влияние линкера и «боеголовки» на стабильность отдельных конформаций лиганда и его динамику. Так, например, если у химерного пептида вероятность конформации лиганда, необходимой для его взаимодействия с рецептором, значительно ниже, чем у отдельного лиганда, то следует оптимизировать последовательность химерного пептида.
5) расчет аффинности отдельного лиганда и химерного пептида к интегрину αVβ6. Если аффинность химерного пептида значительно ниже аффинности отдельного лиганда, на этом шаге также требуется оптимизация последовательности химерного пептида.
6) моделирование химерного пептида в связывающем сайте рецептора с целью определения доступности линкера для катепсина B.

1) Первым силовым полем, которое мы выбрали для моделирования предложенных пептидов, стало поле Amber ff14SB в сочетании с TIP3P водой. Это силовое поле очень часто используется для моделирования белков, и такой выбор является стандартным. Однако в случае лигандов к интегрину αVβ6 использование ff14SB с водой TIP3P привело к появлению очень компактных структур с маленьким радиусом гирации, но в то же время лишенных вторичной структуры. Это явление описано в литературе и является артефактом, связанным с тем, что поле ff14SB в сочетании с водой TIP3P было оптимизировано для моделирования белков, имеющих постоянную вторичную и третичную структуру, а не белков и пептидов с неупорядоченной структурой, к которым, как мы определили впоследствии, относится лиганд к интегрину αVβ6.
Второе силовое поле было выбрано с учетом этого наблюдения. Им стало недавно предложенное поле Amber ff99sb_disp в сочетании со специально оптимизированной под него четырехцентровой водой tp4p-disp (P. Robustelli, S. Piana, D. E. Shaw; PNAS, 2018, 115(21) E4758-E4766). Это поле было оптимизировано в равной мере для упорядоченных белков, для белков с неупорядоченной структурой и для пептидов разной степени упорядоченности. Действительно, моделирование пептидов в этом силовом поле не приводило к появлению компактных конформаций без выраженной вторичной структуры, и мы все последующие моделирования проводили с использованием этого силового поля. Однако следует иметь в виду, что окончательный вывод о применимости этого силового поля как есть, без модификаций, к моделированию исследуемых пептидов может быть сделан только после сравнения окончательных расчетных данных по динамике и термодинамике пептидов с экспериментальными данными.
2) В качестве метода эффективного сэмплирования конформационного пространства для расчета поверхностей свободной энергии пептидов была выбрана метадинамика ( Laio, A.; Parrinello, M.; PNAS, 202, 99(20): 12562–12566.). В этом методе внешний потенциал строится в пространстве нескольких выбранных степеней свободы, называемых коллективными переменными. Коллективными переменными могут быть расстояния между атомами, среднеквадратичное отклонение координат всех атомов от какой-то выбранной конформации, радиус гирации, процент альфа спиралей или бета слоев и т. д. Внешний потенциал обновляется периодически в ходе моделирования, делая энергетически невыгодными те значения коллективных переменных, которые пептид уже принимал. Таким образом, внешний потенциал заставляет пептид исследовать новые области конформационного пространства. Очевидно, что успех моделирования с использованием метадинамики зависит в первую очередь от правильного выбора коллективных переменных. Так, например, если известно, что пептид с некоторой вероятностью образует альфа-спираль, то мера схожести пептида с альфа-спиралью должна быть обязательно среди коллективных переменных, иначе доля альфа-спиральных структур в моделировании будет занижена.
Первый пептид, который был взят для моделирования, - это лиганд, закристаллизованный в комплексе с интегрином αVβ6 (PDB ID 4um9). Его первичная последовательность: HGRGDLGRLKK-NH2. Данный пептид достаточно маленький, и поэтому поверхность его свободной энергии может быть получена, как с использованием методов эффективного сэмплирования конформационного пространства, так и без них, в рамках стандартной молекулярной динамики. Мы решили использовать это обстоятельство и сравнить поверхности свободной энергии, полученные с помощью этих двух подходов. Такое сравнение позволило бы определить, насколько полно выбранные коллективные переменные описывают пептид, и насколько правильно подобраны другие параметры метадинамики. В дальнейшем этот набор коллективных переменных и значения параметров можно будет использовать для моделирования схожих пептидов уже не прибегая к сравнению с классической молекулярной динамикой.
Все расчеты проводились на суперкомпьютерах Cobra и Draco, принадлежащих Институту исследования полимеров.
Стандартная молекулярная динамика: нами было записано 32 траекторий пептида HGRGDLGRLKK-NH2 по 1.5 мкс каждая (суммарное время моделирования 48 мкс; основные параметры моделирования: NPT ансамбль, 300 K, термостат масштабирования скоростей, баростат Парринелло-Рамана). Расчеты были завершены 01.06.2020. На основании полученных данных построена поверхнось свободной энергии конформационных переходов пептида.
Метадинамика: потратив долгое время на опробацию разных коллективных переменных, мы в итоге остановились на семи из них: 1) сходстве с альфа спиралью, определяемом как расстояние между текущими значениями углов φ и ψ всех аминокислотных остатков пептида и набором значений тех же углов в идеальной альфа спирали (φ=-1.0 рад и ψ=-0.82 рад) (см. Plumed коллективная переменная ALPHABETA), 2) сходстве с бета тяжем, определяемом тем же образом (φ=-1.396 рад и ψ=2.618 рад) 3) сходстве с левой альфа спиралью (φ= 1.0 рад и ψ=0.82 рад) 4) сходстве с конформацией, характеризующейся значениями углов: φ=1.396 рад и ψ=2.618 рад. 5) среднеквадратичном отклонении атомов основной цепи пептида от идеального антипараллельного бета листа (см. Plumed коллективная переменная ANTIBETARMSD), 6) среднеквадратичном отклонении атомов основной цепи пептида от идеального параллельного бета листа (см. Plumed коллективная переменная PARABETARMSD), 7) параметре порядка S2, демонстрирующем, в какой мере NH-группы пептида вовлечены в нековалентные внутримолекулярные контакты, а в какой мере сохраняют подвижность. Параметр порядка S2 обычно рассчитывается из скоростей ядерной магнитной релаксации и значений 15N{1H} гетороядерного NOE, но мы использовали эмпирическую формулу, выведенную Брюшвайлером ( Zhang F., Bruschweiler R., J. Am. Chem. Soc., 2002, 124, 12654-12655) и впервые примененную в качестве коллективной переменной Омаром Валссоном (Palazzesi F., Valsson O., Parinello M., J. Phys. Chem. Lett. 2017, 8, 4752−4756).
Тем не менее, семь коллективных переменных - это слишком много для стандартной метадинамики, так как в ней внешний потенциал зависит одновременно от всех коллективных переменных. Мы используем модификацию метадинамики: Parallel biased metadynamics (Pfaendtner J., Bonomi M., J Chem Theory Comput., 2015, 11 (11), 5062-7), в которой к потенциалу, определяемому силовым полем, добавляется ряд одномерных внешних потенциалов, каждый из которых зависит только от одной коллективной переменной. При таком подходе отсутствуют ограничения на количество одновременно используемых коллективных переменных, в результате чего этот метод был успешно применен для моделирования фолдинга минибелка Trp-cage.
На данный момент мы записали 3 мкс траектории с использованием Parallel biased metadynamics, но, судя по предварительному анализу, частоты отдельных конформаций пептида еще не достигли своих предельных значений, и траектория должна быть продолжена.
3) изначально скорости ядерной магнитной релаксации и значения 15N{1H} гетороядерного NOE планировалось рассчитать для пептида HGRGDLGRLKK-NH2. Но оказалось, что для него в литературе отсутствуют соответствующие экспериментальные данные. Такие данные, тем не менее, имеются для пептида FMDV2, имеющего схожую первичную последовательность (NAVPNLRGDLQVLAQKVART) и связывающегося с интегрином αVβ6 с высокой аффинностью и селективностью (Wagstaff J. L., Rowe M. L., Hsieh S-J. et al., RSC Adv., 2012, 2, 11019-11028). Интересно, что FMDV2 не гомологичен первому пептиду: он взят из белка вируса, вызывающего энтеровирусный везикулярный стоматит, а первый пептид является фрагментом естеcтвенного лиганда интегрина αVβ6 – предшественника TGFβ3 (трансформирующего ростового фактора beta 3). Так как FMDV2 является перспективным лигандом интегрина αVβ6 в составе химерного пептида, и для него есть ЯМР данные, мы решили записать с ним ряд траекторий, пока в рамках классической молекулярной динамики, но впоследствии и с использованием методов эффективного сэмплирования. На данный момент нами записано 4 траектории пептида FMDV2 по 1 мкс каждая (суммарное время 4 мкс) при температуре 283 K (эта температура поддерживалась в ЯМР экспериментах), и 4 траектории такой же длины при температуре 300 K (знание о стабильности отдельных конформаций и динамике пептида при этой температуре важно для решения биологических задач). Предварительный анализ траекторий показал, что их недостаточно для того, чтобы частоты всех конформаций FMDV2 достигли своих предельных значений, но ЯМР параметры (15N, 1HN, 1HA химические сдвиги, 15N T1 и T2 времена релаксации, значения 15N{1H} гетероядерного NOE) рассчитанные в первом приближении на основании этих траекторий относительно неплохо согласуются с экспериментальными данными, что говорит о правильном выборе силового поля. Более точный расчет ЯМР параметров возможен после определения предельных частот всех конформаций FMDV2 в растворе с использованием методов эффективного сэмплирования конформационного пространства.
На данный момент полученные результаты готовятся к публикации и ведется работа над остальными пунктами плана