Метадинамика¶

Создайте рабочую директорию.

• Координаты белка возьмите из записи PDB:1L2Y .
• файл настроек для минимизации энергии em.mdp
• файл настроек для "утряски" воды pr.mdp
• файл настроек для молекулярной динамики md.mdp

• Теперь построим файл топологии системы в силовом поле amber99sb и файл с координатами в формате Gromacs.
  

  gmx pdb2gmx -f p.pdb -o p -p  -ff amber99sb -water tip3p
  

• Делаем небольшой отступ в ячейке от белка.
  

  gmx editconf -f p.gro -o p_ec -d 1.5
  

• Проведём оптимизацию геометрии системы, что бы удалить "плохие" контакты в молекуле.
  

  gmx grompp -f em -c p_ec -p  -o p_em -maxwarn 1
  gmx mdrun -deffnm dna_p -v
  

  
  *Отметье в отчёте изменение максимальной силы в ходе оптимизации геометрии. Знаесите начальное и конечное значение максимальной силы.

• Добавим в ячейку молекулы воды
  

  gmx genbox -cp p_em -p  -cs  -o p_s
  

• Нейтрализуем заряд системы. Это делаем в два шага: строим tpr и запускаем genion. В выводе grompp обратите внимание на информацию о заряде системы.
  

  gmx grompp -f em -p  -c p_s -o p_s
  gmx genion -s p_s -o p_si -p  -neutral -conc 0.15
  

• Проведём "утряску" воды:
  

  gmx grompp -f pr -c p_si -p  -o p_pr -maxwarn 1
  gmx mdrun -deffnm p_pr -v
  

• Cравните визуально изменения в системах p_pr.gro и p_si.gro . Занесите наблюдения в отчёт.
• Запускаем тестовое моделирование молекулярной динамики.
  

  gmx grompp -f md -c p_pr -p  -o p_md -maxwarn 1
  gmx mdrun  -deffnm p_md -v -nsteps 1000000 -ntomp 4
  

  
  Просмотреть ход счёта можно в файле p_md.log

Метадинамика¶

• Попробуем странную штуку https://github.com/spiwokv/af2cv/tree/main
• Pickle для белка можно сделать в ColabFold с выбором AlphaFold2_mmseqs2 и обязательно поставить галочку Save All. Pickle будет в архиве.
• Вероятно, что при утряске воды система будет взрываться, надо в em.mdp поменять integrator на steep. И потом:

gmx grompp -f em -c p_si -p  -o p_em2 -maxwarn 1
gmx mdrun -deffnm p_em2 -v
gmx grompp -f pr -c p_em2 -p  -o p_pr -maxwarn 1
gmx mdrun -deffnm p_pr -v

• Скопируйте себе библиотеку для этой CV и файл plumed.dat

• В начале файла надо внести изменения в путь к библиотеке с CV

Для варианта дизайна этого пептида надо догадаться как поменять номера атомов в двух местах: WHOLEMOLECULES и ATOMS=

• Запускаем

gmx  mdrun -deffnm p_md -ntomp 4 -v  -plumed plumed.dat

Анализ метадинамики¶

• В ходе расчета у вас получились файлы со структурами : trr и xtc
• Файл HILLS это информация об наложенных потенциалах, это тот самый "песок"
• Файл COLVAR это файл в который сохраняли значения CV и сумарное наложенный bias.
• Посомотрите на результат в PyMol

        gmx trjconv -f p_md -s p_md -o view1.pdb -skip 20 -pbc mol
        или
        trjconv -f p_md.xtc -s p_md.tpr -o view2.pdb -skip 20 -fit rot+trans

• Воспользуйтесь plumed sum_hills для востановления поверхности свободной энергии
• Используя matplotlib и подобное постройте график
• Попробуйте придумать способ для оценки сходимости
• Давайте попробуем построить график FES на основе перевзвешенных данных.

         plumed driver --plumed rew.dat --noatoms

• Постройте график и сравните с sum_hills
• Возможно это не самая лучшая переменная, давайте попробуем превзвесить на радиус инерции (предположительное распрямление белка)
• Расчитаем RG c помошью этих директив в отдельном dat файле

        rg: GYRATION TYPE=RADIUS ATOMS=1-284
        PRINT ARG=rg STRIDE=1 FILE=rg

• Расчёт новой переменной будет сделан в новый файл

        plumed driver --plumed rg.dat --mf_trr p_md.trr

• Теперь надо скомбинировать данные meta.rbias из COLVAR и значения радиуса гирации в один новый файл формата COLVAR , обратите внимание на разницу в шаге. Вот некоторые функции в помощь

        def GetPlumedColumns(f):
            with open(f,'r') as dat:
            arr = dat.readline().split()
            return arr[2:]

        def LoadNumpy(f,names):
            a = np.loadtxt(f)
            df = pd.DataFrame()
            for i,n in enumerate(names):
            df[n] = a[:,i]
            return copy.deepcopy(df)

        def combinevars(.......):
            c = LoadNumpy(file, cnames )
            gr = LoadNumpy(file, danames )
            out =  pd.merge(crf, da, on='time')

• Расчитайте график свободной энергии в завимисомти от радиуса инерции, надо сделать новый rev.dat
• Попрбуйте сделать анализ от двух переменных RG и alfascore
• Может есть силы поискать другие переменные которые могут лучше описывать процесс : https://www.plumed.org/doc-v2.9/user-doc/html/_colvar.html
• Для выбранной переменной сделайте анализ мутантного белка и ранее стабилизированного белка, данные лежат выше *

Подсказка¶

Общий пайплан построения FES : Подготовить COLVAR (слить, фильтровать) - Адаптировать rew.dat - Запустить plumed driver --plumed rew.dat --noatoms

• На hub запуск Plumed:

        /home/anur/miniconda3/envs/plumed/bin/plumed driver --plumed rew2d.dat --noatoms

• Для работы дома надо подготовить все возможные COLVAR и положить их в свою директорию на hub